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

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

R进⾏超⼏何分布计算

超⼏何分布是统计学中的⼀种离散分布，它描述了由限个物件中抽出n个物件，成功抽出指定种类的物件的个数。也是⽣物信息学中常⽤

的⼀种统计分布模型。我们在GO，KEGG的富集分析时就采⽤的是这个模型。其数学表达式如下：

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

其中：

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.tale(/nfs2/igeneook/old/paotongmerip_AJMeR2191122028MeRIP/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],

=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)