######################################################## 
#-------------------------------------------------------
# Topic:超几何分布检验推导
# Author:Wang Haiquan
# Date:Sat Oct 31 08:53:20 2020
# Mail:mg1835020@smail.nju.edu.cn
#-------------------------------------------------------
########################################################


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

#1、该情况的概率
choose(20,2)*choose(80,8)/choose(100,10)
#[1] 0.3181706
#2、该情况是否为随机的
phyper(q=2,
       m=20,
       n=80,
       10)
#[1] 0.6812201
因为p=0.6812201>0.05,因此认为该过程为随机情况（这里计算的左侧累积概率）

另附其它的集中情况的可能性：
choose(20,0)*choose(80,10)/choose(100,10)#
#[1] 0.09511627
choose(20,1)*choose(80,9)/choose(100,10)#
#[1] 0.2679332
choose(20,2)*choose(80,8)/choose(100,10)#
#[1] 0.3181706
choose(20,3)*choose(80,7)/choose(100,10)#
#[1] 0.2092081
choose(20,4)*choose(80,6)/choose(100,10)
#[1] 0.0841073
choose(20,5)*choose(80,5)/choose(100,10)
#[1] 0.02153147
choose(20,6)*choose(80,4)/choose(100,10)
#[1] 0.00354136
choose(20,7)*choose(80,3)/choose(100,10)
#[1] 0.0003679335
choose(20,8)*choose(80,2)/choose(100,10)
#[1] 2.299585e-05
choose(20,9)*choose(80,1)/choose(100,10)
#[1] 7.762311e-07
choose(20,10)*choose(80,0)/choose(100,10)
#[1] 1.067318e-08

结论 ：

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)

重复clusterprofiler包的结果！

我们先使用clusterprofiler的示例数据

data(geneList, package = "DOSE")
de <- names(geneList)[1:100]
yy <- enrichGO(de, 'org.Hs.eg.db', ont="BP", pvalueCutoff=0.01)
head(yy)
                   ID                          Description GeneRatio   BgRatio       pvalue
GO:0140014 GO:0140014             mitotic nuclear division     28/99 264/18670 5.669344e-29
GO:0000280 GO:0000280                     nuclear division     30/99 407/18670 2.557918e-26
GO:0048285 GO:0048285                    organelle fission     30/99 449/18670 4.637016e-25
GO:0000070 GO:0000070 mitotic sister chromatid segregation     20/99 151/18670 9.796542e-23
GO:0007059 GO:0007059               chromosome segregation     25/99 321/18670 1.735044e-22
GO:0000819 GO:0000819         sister chromatid segregation     21/99 189/18670 3.390547e-22
               p.adjust       qvalue
GO:0140014 9.972377e-26 7.954985e-26
GO:0000280 2.249689e-23 1.794581e-23
GO:0048285 2.718837e-22 2.168822e-22
GO:0000070 4.308029e-20 3.436524e-20
GO:0007059 6.103885e-20 4.869081e-20
GO:0000819 9.939952e-20 7.929120e-20

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)
#   a  b   c     d            p        p.adj
#1 28 99 264 18670 5.669344e-29 9.972377e-26
#2 30 99 407 18670 2.557918e-26 2.249689e-23
#3 30 99 449 18670 4.637016e-25 2.718837e-22
#4 20 99 151 18670 9.796542e-23 4.308029e-20
#5 25 99 321 18670 1.735044e-22 6.103885e-20
#6 21 99 189 18670 3.390547e-22 9.939952e-20
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")
#   a  b   c     d            p        p.adj        p.diy    p.adj.diy
#1 28 99 264 18670 5.669344e-29 9.972377e-26 5.669344e-29 9.972377e-26
#2 30 99 407 18670 2.557918e-26 2.249689e-23 2.557918e-26 2.249689e-23
#3 30 99 449 18670 4.637016e-25 2.718837e-22 4.637016e-25 2.718837e-22
#4 20 99 151 18670 9.796542e-23 4.308029e-20 9.796542e-23 4.308029e-20
#5 25 99 321 18670 1.735044e-22 6.103885e-20 1.735044e-22 6.103885e-20
#6 21 99 189 18670 3.390547e-22 9.939952e-20 3.390547e-22 9.939952e-20

有两点疑惑：
1、q的取值中为何要是28-1 ？
2、lower.tail 参数的定义为 logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
也就是说F代表我们计算的是P[X>x]，也就是计算的右侧累计概率，这么做的目的是什么呢？

在医学统计学课本中，对于确定最终p值得描述为：先求出所有组合四格表得概率，然后将小于等于原样本四格表概率得所有四格表概率值相加，得到双侧检验的P值，原样本的四格表以左（包括原样本）的所有四格表概率之和为左侧概率，以右（包括原样本）的所有四格表概率之和为右侧概率，左侧和右侧概率中较小者为单侧检验的P值

一个原函数的链接：https://www.cnblogs.com/leezx/p/12510385.html