Python下实现Fisher富集分析方法
Fisher's精确概率方法利用超几何分布的原理推断每个基因集中的目标基因(ADR基因)的比例是否与整个基因组中目标基因(ADR基因)的比例相同。包括2个原假设:①基因是否目标基因(ADR基因);②基因是否属于的功能基因集(如GO Term),如表:
功能基因集 非功能基因集
合计
注示:N表示全基因组中基因总数;M表示功能基因集中的基因个数,K表示ADR基因数目.。Fisher's得分表示k个ADR基因中至少有x个被功能基因集注释的概率:
ADR基因
非ADR基因 合计 M N-M N
n11 n21
K
n12 n22
N-K
?M??N?M????x?1?i??K?i??p?1???N?i?0???K?
例如我们要研究药物相关基因和功能term之间的相关性,可以利用下面两种情况
功能基因集 非功能基因集
合计
功能基因集 非功能基因集
合计
接下来我们介绍下如何通过Python语言利用FISHER精确检验实现功能富集分析。
miRNA靶基因 Drug显著基因
非drug显著基因 合计 M N-M
N(全基因组个数)
n11 n21
K
n12 n22
N-K
非miRNA靶基因 合计 M N-M
N(所有miRNA靶基
因个数) n11 n21
K
n12 n22
N-K
在Python语言中,Fisher精确检验主要依赖于scipy.stats包以及numpy包,安装和调用命令如下:
假设我们有一个待查基因集,命名为interest,我们想知道这个基因集和通路path之间的相关性,即这个interest基因集是否可以显著富集到通路path中,通路path中包含基因个数为num_sp,另外还有一个背景基因集,即全部的基因个数,命名为num_allgene。接下来就是如何计算四格表,
path nonpath
interest a b num_interestgene
non interest c d num_allgene-num_interestgene num_sp num_allgene-num_sp
这里a为待查基因集和通路基因集的交集,b为不在通路path的待查基因个数,c为不属于interest待查基因的通路基因个数。计算过程如下
a,b,c,d分别为四个表的四个得分值,红圈内的函数即为利用四格表计算的显著性得分。函数fisher_exact()返回的是一个列表,列表第二个值就是我们要的P值了。 最后我们以P小于0.05作为显著性阈值,通过下面命令进行筛选 if pvalue<0.05:
print (sp[0]+'\\t'+`pvalue`)
以上过程是对于一个通路path而言的,如果对于所有的KEGG通路,或者是GO term,通过一个迭代循环,就可以实现批处理的富集分析。
最后一步就是FDR校正,由于经过对所有KEGG通路或GOterm进行富集分析,fisher精确检验会输出一列P值,往往存在假阳性,因此需要经过FDR校正。这里采用的FDR校正方法为BH校正法。命令如下