生信编程实战第7题（python）

image.png

做这个题目之间必须要了解一些背景知识

1.超几何分布
超几何分布是统计学上一种离散概率分布。它描述了由有限个物件中抽出n个物件，成功抽出指定种类的物件的次数（不归还），称为超几何分布。

2.富集分析的原理
基于筛选的差异基因，或其他自己定义的一组基因，采用超几何检验，判断上调或下调基因在哪些GO或KEGG或其他定义的通路富集。
假设背景基因的数目为m
背景基因中某一通路的pathway中的基因有n个
上调基因有k个
上调基因中落于通路的基因有1个

简单来说，就是比较1/k是否显著高于（也可能低于）n/m，也就是上调基因落在通路中的比例是否显著高于背景基因在这一个通路中的比例

3.例子
举个例子
假设：
一个人类的背景基因有30000(m)个，其中有40(n)个落在p53 signaling pathway上。
在一个给定的基因集中，这个基因集共有300(k)个基因是和背景基因overlap的，
其中3个基因落在p53 signaling pathway上。
这时候就有一个问题：
与背景的40/30000相比，3/300是否显著高于随机概率
这时候:
Fisher Exact P-Value=0.008
因为pvalue <0.01
所以3/300是否显著高于随机概率

但是由于很多时候有多个pathway，这时候就涉及到一个多重假设检验计算FDR

4.背景基因
关于背景基因
收集一
凡是富集分析，都要有背景和选择集
有参的，那就找参考对应的注释信息，作为背景
无参的，那就自己注释，得到背景

收集二
其实pathway富集分析本身也只是提供一些参考，并非非要富集不可。因为某些pathway的调控，基因直接并非相互调控，而是共同参与某个产物合成过程中的不同步骤。例如，某代谢性物X的合成，需要合成酶 A、B、C、D 四个合成步骤。那么A表达的变化，并不会直接影响B、C、D基因的表达，只是影响代谢物X的合成量。如果没有富集到，你就当这个是基因注释了，讨论这些落在你感兴趣的pathway中的基因，也是一种策略。

5.问题解析

所以这道题的关键在于得到4个值
人的kegg pathway中注释的所有基因的数目作为背景基因m
具体到某个通路上的所有基因数目n

基因集中与m重合的基因数目 k
具体到某个通路，k中落在该通路上的基因数目a

得到这个4个值，然后分别计算每个通路的pValue，将得到的pValue进行多重假设检验得到FDR，再筛选即可。

6.代码
基因集来自我自己整理的差异基因集

import osimport reimport pandas as pdimport numpy as npfrom scipy.stats import hypergeomfrom collections import OrderedDict

kegg = OrderedDict()#计算m和kpop=[]for line in open("hsa.kegg.txt"):
    lineL=line.strip().split("\t")    if len(line)>3:     #只有len（line）>3才表示这条是有基因的pathway
        gene_id=lineL[-2]
        pathway=lineL[2]
        gene_idL=gene_id.strip().split(";")
        kegg[pathway]=gene_idL        for i in gene_idL:            if i not in pop:
                pop.append(i)
DEG=[]for line in open("ehbio.DESeq2.all.DE.entrez.txt") :
    line=line.strip()    if line in pop:
        DEG.append(line)
pop_number=len(pop)
DEG_number=len(DEG)

print("the number of all gene in pathway is :%d" %pop_number )#注意%d而不是d%print("the number of diff gene in pathway is :%d" %DEG_number)#超几何分布检验keggEnrichment=OrderedDict()for k,v in kegg.items():
    cnt=[]
    pathway_gene_cnt=len(v)    for i in DEG:        if i in v:
            cnt.append(i)
    dif_in_pathway=len(cnt)    if dif_in_pathway == 0:
        print(k)    else:
        pValue=hypergeom.sf(dif_in_pathway-1,pop_number,pathway_gene_cnt,DEG_number)
        keggEnrichment[k]=[dif_in_pathway-1,pop_number,pathway_gene_cnt,DEG_number,pValue,";".join(cnt)]#用pandas的dataframe便于操作keggOUT=pd.DataFrame.from_dict(keggEnrichment,orient="columns",dtype=None)
keggOUT=pd.DataFrame.transpose(keggOUT) #行列的转置keggOUT.columns=["a","m","n","k","pValue","gene"]
keggOUT=keggOUT.sort_values(by="pValue",axis=0)#FDRdef p_adjust_bh(p):
    """Benjamini-Hochberg p-value correction for multiple hypothesis testing."""
    p = np.asfarray(p)
    by_descend = p.argsort()[::-1]
    by_orig = by_descend.argsort()
    steps = float(len(p)) / np.arange(len(p), 0, -1)
    q = np.minimum(1, np.minimum.accumulate(steps * p[by_descend]))    return q[by_orig]

pValue=keggOUT["pValue"]
fdr=p_adjust_bh(pValue)
fdr=pd.DataFrame(fdr,index=keggOUT.index)
keggOUT.insert(5,"FDR",fdr)

keggOUT.loc[keggOUT["FDR"]<0.05] #筛选FDR<0.05

结果如下

    a   m   n   k   pValue  FDR genehsa05200:Pathways in cancer 47  13712   526 571 2.59037e-07 0.000086    9252;6387;624;196883;815;7042;1030;3162;623;11...hsa04360:Axon guidance  22  13712   175 571 1.00737e-06 0.000167    6387;6091;815;10371;80031;2049;2048;8503;21969...hsa04923:Regulation of lipolysis in adipocyte   11  13712   54  571 1.67365e-06 0.000185    196883;114;5593;8503;134;5743;364;57104;8660;5...hsa04052:Cytokines and growth factors   25  13712   237 571 6.35655e-06 0.000528    85480;3976;6387;9518;10673;146433;7042;9966;82...hsa00536:Glycosaminoglycan binding proteins 22  13712   200 571 9.98448e-06 0.000663    1277;6387;7042;7130;7422;3592;2263;7472;255631...hsa04933:AGE-RAGE signaling pathway in diabetic complications   14  13712   99  571 1.3803e-05  0.000764    1277;7412;5581;7042;8503;7422;3725;2308;7048;5...hsa04750:Inflammatory mediator regulation of TRP channels   13  13712   99  571 5.83834e-05 0.002769    5321;624;196883;5581;815;623;114;8503;4803;40;...hsa04068:FoxO signaling pathway 15  13712   132 571 0.000118444 0.004638    1901;7042;1030;8503;1026;10769;8743;10912;8660...hsa05224:Breast cancer  16  13712   147 571 0.000131748 0.004638    672;8503;1026;7472;3725;3714;2099;2255;10912;8...hsa04512:ECM-receptor interaction   11  13712   82  571 0.000139684 0.004638    1277;3673;3680;3908;7057;1282;3161;8515;1286;1...hsa04350:TGF-beta signaling pathway 11  13712   84  571 0.000176636 0.005331    10468;7042;1030;8200;83729;130399;7048;7057;36...hsa04213:Longevity regulating pathway - multiple species    9   13712   62  571 0.000220627 0.006104    196883;114;8503;3310;8660;2309;2308;51422;5295...hsa04060:Cytokine-cytokine receptor interaction 25  13712   294 571 0.000248419 0.006344    85480;3976;6387;9518;10673;146433;7042;9966;82...hsa04926:Relaxin signaling pathway  14  13712   130 571 0.000330597 0.007462    1277;9586;196883;114;8503;7422;59350;3725;5433...hsa05414:Dilated cardiomyopathy (DCM)   11  13712   90  571 0.000341367 0.007462    196883;7042;114;3673;6444;3680;783;3908;488;85...hsa04977:Vitamin digestion and absorption   5   13712   24  571 0.000359612 0.007462    8029;8884;151056;80704;10560;338hsa05222:Small cell lung cancer 11  13712   93  571 0.000463659 0.009055    1030;8503;3673;1026;5743;10912;4616;3908;5295;...hsa00350:Tyrosine metabolism    6   13712   36  571 0.000609243 0.010704    259307;218;4128;125;316;130;4129hsa04010:MAPK signaling pathway 24  13712   295 571 0.000612603 0.010704    9252;5321;1649;11221;7042;7422;3310;2263;3725;...hsa04510:Focal adhesion 18  13712   199 571 0.000655508 0.010881    1277;8503;7422;3673;3725;5228;5063;3680;3908;7...hsa05412:Arrhythmogenic right ventricular cardiomyopathy (ARVC) 9   13712   72  571 0.000758135 0.011649    3673;6444;3680;783;3908;488;8515;8516;1756;5318hsa05410:Hypertrophic cardiomyopathy (HCM)  10  13712   85  571 0.00077191  0.011649    7042;3673;6444;3680;783;3908;488;51422;8515;85...hsa04211:Longevity regulating pathway - mammal  10  13712   89  571 0.00113772  0.016423    9586;814;196883;114;8503;8660;2309;2308;51422;...hsa00750:Vitamin B6 metabolism  2   13712   6   571 0.00130739  0.017691    29968;316;493911hsa04666:Fc gamma R-mediated phagocytosis   10  13712   91  571 0.00136822  0.017691    4082;65108;5321;5581;8613;8503;3635;10810;273;...hsa05226:Gastric cancer 14  13712   149 571 0.00138548  0.017691    7042;1030;8503;1026;2263;7472;2255;10912;8325;...hsa99995:Signaling proteins 15  13712   165 571 0.00144771  0.017801    26045;4974;23767;9628;5999;349667;323;7079;481...hsa00410:beta-Alanine metabolism    5   13712   31  571 0.00153551  0.018207    218;339896;18;4329;219;51733hsa04390:Hippo signaling pathway    14  13712   154 571 0.00192623  0.021526    154796;7042;8200;7472;84552;23418;8325;7048;85...hsa05031:Amphetamine addiction  8   13712   68  571 0.00194509  0.021526    9586;814;815;2903;84152;3725;4128;5500;4129hsa04713:Circadian entrainment  10  13712   96  571 0.00211575  0.022659    8863;9252;196883;815;114;5593;2903;54331;8864;...hsa04015:Rap1 signaling pathway 17  13712   206 571 0.00244447  0.025361    196883;57568;114;8503;7422;2903;11069;2263;480...hsa04380:Osteoclast differentiation 12  13712   128 571 0.00260978  0.026256    814;7042;8503;4982;8651;29760;3725;9103;54;704...hsa05210:Colorectal cancer  9   13712   86  571 0.00297969  0.029096    7042;8503;1026;3725;10912;7048;4616;5881;5295;...hsa04022:cGMP - PKG signaling pathway   14  13712   163 571 0.00334068  0.031689    9586;624;196883;5581;8654;114;5593;150;134;866...hsa04024:cAMP signaling pathway 16  13712   198 571 0.00381188  0.033934    9586;814;196883;815;114;8503;2903;11069;6662;1...hsa05212:Pancreatic cancer  8   13712   75  571 0.00383541  0.033934    7042;8503;7422;1026;10912;7048;4616;5881;5295hsa01001:Protein kinases    34  13712   523 571 0.00388403  0.033934    152110;9252;11113;79858;814;23043;5581;815;244...hsa04020:Calcium signaling pathway  15  13712   183 571 0.00412696  0.034471    5027;814;624;196883;57620;815;623;114;2903;891...hsa04974:Protein digestion and absorption   9   13712   90  571 0.00415309  0.034471    1277;10008;5222;255631;54407;1294;1282;1286;12...hsa04210:Apoptosis  12  13712   136 571 0.00441474  0.035088    7277;1649;7846;8503;8743;3725;4803;10912;4616;...hsa05225:Hepatocellular carcinoma   14  13712   168 571 0.00443879  0.035088    7042;3162;8503;1026;7472;10912;8325;7048;4616;...hsa00360:Phenylalanine metabolism   3   13712   17  571 0.00459179  0.035453    259307;218;4128;4129hsa04978:Mineral absorption 6   13712   51  571 0.00493439  0.037232    3162;4502;26872;261729;4493;4501;55503hsa04031:GTP-binding proteins   15  13712   192 571 0.006532    0.046618    338382;5912;8153;6236;23551;51285;57381;54331;...hsa05146:Amoebiasis 9   13712   96  571 0.00656538  0.046618    338382;1277;7042;8503;3592;3908;5295;1282;1286...hsa05219:Bladder cancer 5   13712   41  571 0.00659953  0.046618    9252;7422;1026;7057;1890;23604hsa04151:PI3K-Akt signaling pathway 24  13712   354 571 0.00711077  0.048196    1277;9586;672;8503;7422;3673;1026;2263;54331;4...hsa04261:Adrenergic signaling in cardiomyocytes 12  13712   144 571 0.00711328  0.048196    9586;9252;196883;815;6324;114;11069;783;147;48...hsa01522:Endocrine resistance   9   13712   98  571 0.00757401  0.049299    196883;114;8503;1026;3725;3714;2099;8202;5295;...hsa05211:Renal cell carcinoma   7   13712   69  571 0.00770521  0.049299    7042;8503;7422;1026;3725;5063;112399;5295hsa04810:Regulation of actin cytoskeleton   16  13712   213 571 0.00784549  0.049299    6387;624;623;8503;3673;26999;2263;5063;2255;36...hsa04072:Phospholipase D signaling pathway  12  13712   146 571 0.00795893  0.049299    1759;5321;196883;8613;114;8503;11069;1607;9265...hsa05165:Human papillomavirus infection 23  13712   339 571 0.0080185   0.049299    1277;9586;8503;7422;3673;3955;1026;7472;84552;...

7.验证一下结果对不对
用一个在线工具

image.png

image.png

image.png

image.png

image.png

可以看到这个在线工具最终的结果是59个通路
而我的脚本得到的是55个通路富集

大部分是一样的
有一些不同
可以看到这个工具算出来的背景基因是6910，而我算出来的是13712
猜测有可能是使用的kegg版本不同

验证一下
就拿pathway in cancer这条通路来看

image.png

在线工具计算的基因有37个
而我计算的47

可以看看这相差的10个在不在我下载的kegg数据库中就知道了

写个脚本找到这几个不同的基因,工具的结果存到1.txt

import sys
args=sys.argv
filename=args[1]
tmp=[9252,6387,624,196883,815,7042,1030,3162,623,114,8503,7422,3673,1026,3592,2263,7472,3725,54331,3714,2099,5228,5743,112399,2255,10912,7704,8325,2308,7048,4616,7296,3908,8202,5881,7855,10235,5295,1282,8313,1286,1285,27006,5336,23604,5468,54567,185]
tmp1=[]for line in open(filename):
   lineL=line.strip().split("\t")
   gene_ID=int(lineL[0])   if gene_ID not in tmp1:
        tmp1.append(gene_ID)for i in tmp:   if i not in tmp1:      print (i)

python3 tmp.py 1.txt

9252
815
3162
3592
3714
2099
10912
4616
7296
8202
54567

可以查到这些基因确实在我下载kegg数据库的cancer in pathway的通路中，所以确实是因为kegg版本的原因

8.知识点总结

(1).python中计算p-value用的是scipy.stats中的hypergeom.sf
(2).pandas中dataframe生成的几种方式：

image.png

作者：天秤座的机器狗
链接：https://www.jianshu.com/p/54088c519017