生信编程实战第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 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

可以看到这个在线工具最终的结果是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

可以查到这些基因确实在我下载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

最后编辑于：2018.08.26 20:11:54

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 202,529评论 5赞 475
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,015评论 2赞 379
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 149,409评论 0赞 335
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,385评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,387评论 5赞 364
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,466评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,880评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,528评论 0赞 256
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,727评论 1赞 295
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,528评论 2赞 319
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,602评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,302评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,873评论 3赞 306
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,890评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,132评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,777评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,310评论 2赞 342

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

推荐阅读更多精彩内容