该博文更新已整理到新地址：

http://qinqianshan.com/math/statics_topic/basic-concepts-in-statistics/

经常看到一些饼图，描述某些事物的组成，比如说有钱人的学历分布，然后我们可以看到高学历所占比例并不高，根据这个比例下结论通常是错的，这些比例说明不了问题，如果把各种学历在总体人口中的分布做为背景进行考虑的话，你就会发现学历还是有点用的。

 

当我们用组学测定了一大堆分子之后，我们希望站在更高的角度去看这些分子和那些生物学过程相关。那么通常各种注释，对这些基因/蛋白进行分类，那么从分类的比例上，是不能草率下结论，正如上面有钱人学历分布的例子一样。我们需要把总体的分布考虑进去。

在做富集分析的时候，会涉及到这么一个概念。

 

超几何分布是统计学上一种离散概率分布。它描述了由有限个物件中抽出n个物件，成功抽出指定种类的物件的个数（不归还）。

   超几何分布和Fisher's Exact Test是完全一模一样的原理，只是两种不同的称谓。 

例如在有N个样本，其中m个是不及格的。超几何分布描述了在该N个样本中抽出n个，其中k个是不合格的的机率：

 

上式可如此理解：(_n^N)表示所有在N个样本中抽出n个，而抽出的结果不一样的数目。（_k^m）表示在m个样本中，抽出k个的方法数目。剩下来的样本都是及格的，而及格的样本有N-m个，剩下的抽法便有(_n-K^N-m)种。

若n=1，超几何分布还原为伯努利分布。

若N接近∞，超几何分布可视为二项分布。

 p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k) for x = 0, …, k. 其中, m 是桶里面白球的个数, n 是黑球的个数, k 是从桶中随机取出的球数, x 是取出球 中白球的个数.

然后计算得到的p-value通过Bonferroni校正之后，以0.05为阈值（小于0.05），满足此条件的GO term定义为显著富集。

（1）超几何分布的模型是不放回抽样

（2）超几何分布中的参数是M,N,n上述超几何分布记作X~H(n，M，N）。

 以文章Gene Expre ssion in Ovarian Cancer Reflects Both Morphology and Biological Behavior, Distinguishing Clear Cell from Other Poor-Prognosis Ovarian Carcinomas所鉴定的差异基因为例。

73个差异基因的Symbol，我把它转为entrezgene ID得到57个（漏掉的不管它，只是做为一个例子）:

> eg

 [1] "7980"   "3081"   "3162"   "3059"   "1545"   "1917"   "6696"   "5797" 

 [9] "6648"   "10397"  "6781"   "5817"   "1282"   "1284"   "6948"   "7077" 

[17] "5744"   "8566"   "1368"   "1474"   "11015"  "3306"   "728441" "2678" 

[25] "4286"   "3929"   "5095"   "2064"   "1428"   "6590"   "3569"   "2745" 

[33] "3912"   "978"    "5950"   "6539"   "9445"   "5004"   "9971"   "7453" 

[41] "2719"   "1410"   "1311"   "4653"   "4162"   "5358"   "3484"   "3486" 

[49] "2261"   "307"    "1672"   "4837"   "22795"  "486"    "4118"   "3915" 

[57] "10140"

下面测试一下这些基因和化学刺激响应的相关性。

goid <- “GO:0042221” # response to chemical stimulus

那么做为背景，总体基因为N，属于“化学刺激响应”这个分类的基因有M个。

 

allgeneInCategory <- unique(get(goid, org.Hs.egGO2ALLEGS))

M <- length(allgeneInCategory)

N <- length(mappedkeys(org.Hs.egGO))

样本的大小是n，属于“化学刺激响应”这个分类的基因有k个。

n <- length(eg)

k <- sum(eg %in% allgeneInCategory)

白球黑球问题，最简单的莫过于用二项式分布，从总体上看，要拿到一个基因属于“化学刺激响应”这个分类的概率是M/N。那么现在抽了n个基因，里面有k个基于这个分类，p值为：

> 1-sum(sapply(0:k-1, function(i) choose(n,i) * (M/N)^i * (1-M/N)^(n-i)))

[1] 8.561432e-10

 

二项式分布，是有放回的抽样，你可以多次抽到同一基因，这是不符合的。所以这个计算只能说是做为近似的估计值，无放回的抽样，符合超几何分布，通过超几何分布的计算，p值为：

> phyper(k-1,M, N-M, n, lower.tail=FALSE)

[1] 7.879194e-10

 

如果用2x2表做独立性分析，表如下：

> d <- data.frame(gene.not.interest=c(M-k, N-M-n+k), gene.in.interest=c(k, n-k))

> row.names(d) <- c("In_category", "not_in_category")

> d

                gene.not.interest gene.in.interest

In_category                  2613               28

not_in_category             15310               29

这个也有很多方法可以做检验，经典的有卡方检验和fisher's exact test。

 

如果用卡方检验来做

> chisq.test(d, )

         Pearson's Chi-squared test with Yates' continuity correction

data:  d

X-squared = 51.3849, df = 1, p-value = 7.592e-13

对于２x２表来说，卡方检验通常也只能做为近似估计值，特别是当sample size或expected ell count比较小的时候，计算并不准确。

 

fisher's exact test

，名副其实，真的就比较exact，因为它使用的是超几何分布来计算p值。这也是为什么fisher's exact test和超几何模式计算的p-值是一样的，

> fisher.test(d)

         Fisher's Exact Test for Count Data

data:  d

p-value = 7.879e-10

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

 0.1013210 0.3089718

sample estimates:

odds ratio

 0.1767937

 

通常各种软件做GO富集性分析，都是使用超几何分布进行计算。

富集性分析应用范围非常广，从Disease Ontology, Gene Ontology, KEGG, 到Reactome Pathway等等。

富集分析需要用到的工具：

GOstat

跟 david 和 onto——tool比较来说了

onto-tool

GOEAST

wego

用go slim

GSEA貌似是最有名的软件，但是搞了半天，发现是基于芯片的分析软件，各种输入文件的要求

DAVID直接通过gene id或者symbol做GO富集分析

参考资料：

http://zh.wikipedia.org/wiki/超几何分布

YGC   http://ygc.name/2012/04/28/enrichment-analysis/

http://www.chenlianfu.com/?p=1122

 

ps:这不就是我高中的时候最喜欢的概率分布问题吗？这个概念吓人，原来这么简单。