生信编程实战第7题（python）
2018-08-20 本文已影响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
