详解富集分析超几何分布检验

R进行超几何分布计算

超几何分布是统计学中的一种离散分布，它描述了由有限个物件中抽出n个物件，成功抽出指定种类的物件的个数。 也是生物信息学中常用的一种统计分布模型。我们在GO，KEGG的富集分析时就采用的是这个模型。 其数学表达式如下：

f(k;n,m,N)={{{m \choose k}{{N-m} \choose {n-k}}} \over {N \choose n}}.

其中：
N： 总样本数
m： 总样本中“特定类别”的数量
n: 从总样本中随机抽取的数量
k: 抽取样本中属于“特定类别”的数量
但是我们做统计模型时，统计不是分布概率，而是累积概率，所以计算的p_value 值计算公式为：

对应的R代码实现为：

1-phyper(k-1,m, N-m, n)

富集分析中的超几何分布概率p值等计算

choose(5,2) # 排列组合数
gamma(5) # 阶乘数

#我们通过一次分析，从100个总基因(其中位于A通路的有20个)中筛到了10个差异基因，其中位于A通路的基因有2个，不位于A中的有8个
#问:1、出现该情况的概率是多少？
#2、该抽样过程是否是随机情况？
#1、该情况的概率：

choose(20,2)*choose(80,8)/choose(100,10)
#[1] 0.3181706

#2、该情况是否为随机的：

1 - phyper(q=2,
       m=20,
       n=80,
       10)

1-phyper(k-1,m, N-m, n)

结论：

phyper(q=2,m=20,n=80,k=10)
##### 等价于
choose(20,0)*choose(80,10)/choose(100,10)+choose(20,1)*choose(80,9)/choose(100,10)+choose(20,2)*choose(80,8)/choose(100,10)

实战示例：

yy = read.table('/nfs2/igenebook/old/paotongmerip_AJMeR2191122028_20200506/MeRIP/06.GOKEGG/BB_IP_3_vs_BB_input_3_KEGGEnrich_enrich_filt.txt',header = T,sep = '\t')
yy[1:4,1:6]
xx<-data.frame(a=str_split(yy$GeneRatio,"\\/",simplify = T)[,1],
               b=str_split(yy$GeneRatio,"\\/",simplify = T)[,2],
               c=str_split(yy$BgRatio,"\\/",simplify = T)[,1],
               d=str_split(yy$BgRatio,"\\/",simplify = T)[,2],
               p=yy$pvalue,
               p.adj=yy$p.adjust)

xx<-apply(xx,2,as.numeric)
xx<-as.data.frame(xx)
#注意！！！我这里使用的是q-1 ,同时设置了lower.tail=F
xx$p.diy<-apply(xx,1,function(a){phyper(a[1]-1,a[3],a[4]-a[3],a[2],lower.tail = F)})
xx$p.adj.diy=p.adjust(xx$p.diy,"BH")

head(xx)

注：
phyper(q=131,m=268,n=7662-268,k=1781,lower.tail = F) 
## 等价于  
1-phyper(q=131,m=268,n=7662-268,k=1781)