JS散度评估特征距离

直接上代码:

import pandas as pd
import numpy as np
import os,sys
from sklearn.preprocessing import StandardScaler
from scipy.stats import shapiro, ttest_ind, wilcoxon,levene,kruskal,mannwhitneyu
from statistics import median
from statistics import mean
from scipy.stats import normaltest
import itertools
import random
from scipy.spatial.distance import jensenshannon
from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
import pickle
os.makedirs('./Js_plot', exist_ok=True)
scaler = StandardScaler()
work="./"
tsv_files = [f for f in os.listdir(work) if f.endswith('LE5X.csv')]
tsv_files = [i  for i in tsv_files if i.find('Frag2023_FragArm2023')==-1]
def js_divergence_scipy(p, q,p_SampleID,q_SampleID,name):
    data1 =np.array(np.nan_to_num(p), dtype=np.float64)
    data2=np.array(np.nan_to_num(q), dtype=np.float64)
    kde1 = gaussian_kde(data1)
    kde2 = gaussian_kde(data2)
    x_range = np.linspace(min(data1.min(), data2.min()), max(data1.max(), data2.max()), 1000)
    prob_dist1 = kde1.evaluate(x_range)
    prob_dist2 = kde2.evaluate(x_range)
    out=jensenshannon(prob_dist1, prob_dist2)
    plt.clf()
    plt.figure(figsize=(10, 6))
    plt.plot(x_range, prob_dist1, 'r', label=p_SampleID)
    plt.plot(x_range, prob_dist2, 'b', label=q_SampleID)
    plt.xlabel('Feature Value')
    plt.ylabel('Probability Density')
    plt.title(name.split('.')[0]+' Probability Distribution')
    plt.legend()
    F=p_SampleID+"_"+q_SampleID+"_"+name.split('.')[0]
    plt.savefig('./Js_plot/JSD_' + F + '.png')
    if(out >=0):
        F=p_SampleID+"_"+q_SampleID+"_"+name.split('.')[0]
        with open('./Js_plot/'+F,'w+') as f:
            f.write(str(out)+"\n")
        return out
    return  None
def paired_delta(df1, df2):
    allsum=0
    # 删除元素
    list1,list2=def_na(df1, df2)
    m=len(list1)
    n=len(list2)
    if(m*m==0 or m==0 or n==0 or m!=n):
        return []
    for i in range(0,m):
        allsum=computer(list1[i],list2[i])+allsum
    d=allsum/(m)
    return abs(d)
def def_na(list1,list2):
    nan_idx1 = [i for i, v in enumerate(list1) if np.isnan(v)]
    nan_idx2 = [i for i, v in enumerate(list2) if np.isnan(v)]
    # 合并两个列表的缺失索引
    nan_idx = list(set(nan_idx1 + nan_idx2))
    # 删除元素
    list1 = [v for i, v in enumerate(list1) if i not in nan_idx]
    list2 = [v for i, v in enumerate(list2) if i not in nan_idx]
    return list1,list2
def tmpdffun(tmpdf):
    scaler = StandardScaler()
    grouptmp = tmpdf[['SampleID', 'group']].copy()
    features = tmpdf.drop(['SampleID', 'group'], axis=1)
    scaled_features = scaler.fit_transform(features)
    scaled_df = pd.concat([grouptmp, pd.DataFrame(scaled_features, columns=features.columns)], axis=1)
    return scaled_df
groupname_list=sys.argv[1]
groupFlag=sys.argv[2] #'TypeOfCollection'
#group1=pd.read_table(work+"group1.TypeOfCollection.info.list")
group1=pd.read_table(work+groupname_list)
if 'LabID' not in group1.columns.to_list():
    group1['LabID']=group1['SampleID'].map(lambda x:x.split("-")[0])

group1=group1[['SampleID',groupFlag,'LabID']]  ###LabID 是具体配对编号
print(set(group1[groupFlag]))
group1name=groupFlag
group1.columns=['SampleID','group','LabID']
mygroup_name=group1
mygroup_name_list=list(set(mygroup_name['group']))[0]
mygroup_name_dict_a=[]
mygroup_name_dict_b=[]
def paired_dict(tmpdf):
    a=tmpdf.query("group==@mygroup_name_list")['SampleID'].to_list()[0]
    b=tmpdf.query("group!=@mygroup_name_list")['SampleID'].to_list()[0]
    mygroup_name_dict_a.append(a)
    mygroup_name_dict_b.append(b)
mygroup_name.groupby("LabID").apply(paired_dict)

mygroup_name=mygroup_name[['SampleID','group']]
mydict={}
for name in tsv_files:
    print(name)
    mydict[name]=[]
    datafile=work+name
    frature=pd.read_csv(datafile,header=0)
    df=pd.merge(mygroup_name,frature,on=['SampleID'])
    df=df.dropna(axis=1)
    #df=df.dropna(axis=1,thresh=round( df.shape[0] / 10) )  ##列NA限制10%
    df1=tmpdffun(df)  ## scale
   # df1=df
    groupA=df1.query("SampleID ==@mygroup_name_dict_a")
    groupA['SampleID'] = pd.Categorical(groupA['SampleID'], categories=mygroup_name_dict_a, ordered=True)
    groupA = groupA.sort_values('SampleID')
    groupB=df1.query("SampleID ==@mygroup_name_dict_b")
    groupB['SampleID'] = pd.Categorical(groupB['SampleID'], categories=mygroup_name_dict_b, ordered=True)
    groupB = groupB.sort_values('SampleID')
    for lengroup in range(groupB.shape[0]):
        data_b=np.array(groupA.iloc[lengroup,2:])
        b_SampleID=groupA.iloc[lengroup,0]
        data_p=np.array(groupB.iloc[lengroup,2:])
        p_SampleID=groupB.iloc[lengroup,0]
        JSdis=js_divergence_scipy(data_b,data_p,b_SampleID,p_SampleID,name)
        if(JSdis is None):
            pass
            #mydict[name].append(0)
        else:
            mydict[name].append(JSdis)
name=[]
value=[]
for k,v in mydict.items():
    for j in v:
        name.append(k.replace(".LE5X.csv","").split(".")[-1])
        value.append(j)
newdf=pd.DataFrame({'Categories': name, 'Values': value})
newdf['Categories']=newdf['Categories'].map(lambda x:x.replace("854TF_OCF_Tcell","griffin_ocf"))
newdf['group']=groupFlag
newdf.to_csv(work+"JS_"+groupFlag,index=False,sep="\t")
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,496评论 6 501
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,407评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,632评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,180评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,198评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,165评论 1 299
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,052评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,910评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,324评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,542评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,711评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,424评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,017评论 3 326
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,668评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,823评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,722评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,611评论 2 353

推荐阅读更多精彩内容