从Fisher's exact test, pvalue, FDR 谈起

事实上，我有一个形似这样的热图：

即一组 n 行 m 列的数据，在热图上，每一个小格代表每个数据占这一行总数的比例。

之前做到这里把热图一画，看看哪几个格子颜色深就得了，也想做一些统计学的比较，但似乎无从入手。

这几天看文献发现可以做一个 Fisher’s Exact Test ，其实也就是常说的富集分析，平时做的 GO 富集也就是这个道理，只是我没有想到而已。

什么是 Fisher’s exact test 呢？

Fisher’s exact test is a statistical significance test used in the analysis of contingency tables. Although in practice it is employed when sample sizes are small, it is valid for all sample sizes

— From Wikipedia

即是用于分析列联表 contingency tables 统计显著性检验方法，它用于检验两个分类的关联 association 。虽然实际中常常使用于小数据情况，但同样适用于大样本的情况。

	Men	Women	Row Total
Studying	a	b	a+b
Not-studying	c	d	c+d
Col Total	a+c	b+d	a+b+c+d=n

对一个2乘2列联表进行检验，使用公式：

这个 p 值即是这个表出现的概率，而 pvalue 表示在原假设为真前提下，出现该样本或比该样本更极端的概率之和，即 pvalue 是多个表 p 值的和。

现在转到生物学分析来，比如在做 GO 富集分析时，某样品 GO 注释的总基因数为 N ，属于某个条目的基因数为 M ，差异表达分析后得到的基因数量为 n ，这 n 个基因里有 k 个基因属于这个条目，用这些数据判断此条目是否富集。
关于富集性分析，有几种方法可以执行：

超几何分布（常用）
fisher’s exact test （也使用超几何分布）
卡方检验

这里超几何分布和 fisher’s exact test 是一个意思，fisher’s exact test 是使用超几何分布来计算 P 值，所以它们的结果是相同的，而卡方检验通常只做近似值估计，并不精确，不推荐使用。

试着用 fisher’s exact test 进行富集分析，先建立 2×2 列联表：

	gene.in.interest	gene.not.in.interest
In.category	k	M-k
Not.in.category	n-k	N-M-n+k

在 R 中进行 fisher’s exact test ：

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

得到 pvalue 。

我比较用 Python ，在 Python 中也有相应的函数，如下：

>>> from fisher import pvalue
>>> mat = [[12, 5], [29, 2]]
>>> p = pvalue(12, 5, 29, 2)
>>> p.left_tail, p.right_tail, p.two_tail
# doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS 
# (0.04455473783507..., 0.994525206021..., 0.0802685520741...)

有两侧还有单侧的结果，选择合适的即可。

现在我们得到了这一 GO 条目的 P 值，仅仅是第一步，假设你根据 P 值挑选出了显著富集的 R 个 GO 条目，其中只有 S 个是真正富集的，另外的 V 个事实上是没有富集的，也就是假阳性，因此我们希望错误比例 Q＝V/R 平均而言不能超过某个预先设定的值（比如 0.05 ），在统计学上，这也就等价于控制 FDR 不能超过 5％。

这里又提到一个问题，所谓 Pvalue、Qvalue、FDR、adjusted Pvalue究竟有啥区别，Pvalue之前介绍过了，即衡量假阳性率的指标，FDR如上所说，是衡量错误发现率的指标，即：使用Qvalue的这个参数预估FDR。Qvalue 是利用公式从Pvalue校正计算后得到，所以Qvalue通常又被称为adjusted Pvalue。

所以从某种程度来说，Qvalue、FDR、adjusted Pvalue是等价的。

接下来计算 FDR，R中的实现：

> p<-c(0.0003,0.0001,0.02) 
> p 
# [1] 3e-04 1e-04 2e-02 
> p.adjust(p,method="fdr",length(p)) 
# [1] 0.00045 0.00030 0.02000

在 Python 中使用 statsmodels.sandbox.stats.multicomp.multipletests 实现，这里就不细说了。

最后回到文初的问题，我把基因分成了 n 类，n 类的基因分别又注释到了 m 个条目类别里，现在拥有每类基因注释到每个条目的基因数 a，即是每个小格的意义，对每个小格做fisher’s exact test，得到P值，最后对每列，即每类基因集的所有P值做FDR检验，设定FDR阈值，即可获得富集小格。

写的比较急，可能有一些错误，请指正。

参考:

超几何分析和 GO 富集分析
http://www.bakerwm.com/r/2015/01/23/hypergeometric-analysis—enrichment-analysis/

Fishers Exact Test for Python (Cython)
https://github.com/brentp/fishers_exact_test

Statsmodels is a Python module
http://www.statsmodels.org/devel/generated/statsmodels.sandbox.stats.multicomp.multipletests.html#statsmodels-sandbox-stats-multicomp-multipletests

多重检验中的 FDR 错误控制方法与 p-value 的校正及 Bonferroni 校正
https://wenku.baidu.com/view/c0008226a58da0116d17492e.html