美文网首页生物信息学与算法
生信编程实战第7题(python)

生信编程实战第7题(python)

作者: 天秤座的机器狗 | 来源:发表于2018-08-20 22:12 被阅读35次

    题目来自生信技能树论坛

    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 os
    import re
    import pandas as pd
    import numpy as np
    from scipy.stats import hypergeom
    from collections import OrderedDict
    
    kegg = OrderedDict()
    
    #计算m和k
    
    pop=[]
    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)
    
    
    #FDR
    def 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 gene
    hsa05200: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;338
    hsa05222: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;4129
    hsa04010: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;5318
    hsa05410: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;493911
    hsa04666: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;51733
    hsa04390: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;4129
    hsa04713: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;5295
    hsa01001: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;4129
    hsa04978:Mineral absorption 6   13712   51  571 0.00493439  0.037232    3162;4502;26872;261729;4493;4501;55503
    hsa04031: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;23604
    hsa04151: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;5295
    hsa04810: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

    (3).dataframe的行列转置的方法
    pd.DataFrame.transpose()
    (4).给dataframe添加列名
    keggOUT.columns=["a","m","n","k","pValue","gene"]
    (5).dataframe按某一列的值排序
    keggOUT=keggOUT.sort_values(by="pValue",axis=0)
    (6).dataframe删除某一列
    keggOUT=keggOUT.drop("FDR",1)
    (7).将array转成dataframe
    fdr=pd.DataFrame(fdr,index=keggOUT.index)
    (8).dataframe按条件筛选 行
    keggOUT.loc[keggOUT["FDR"]<0.05]
    (9).富集分析一般是p<0.01
    FDR<0.05

    相关文章

      网友评论

        本文标题:生信编程实战第7题(python)

        本文链接:https://www.haomeiwen.com/subject/oelriftx.html