校正

p.adjust()

library("fdrtool")

https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html

http://www.360doc.com/content/17/1228/11/50153987_717073620.shtml

http://www.360doc.com/content/18/0914/23/52043738_786760347.shtml

在生信分析中，隔三差五地就需要和p值探讨是否显著差异，还要搬出FDR对p值进行校正。让每个基因根据p值大小从小到大排个队，拿个号牌，然后把自己的p值乘上总基因数，再除以自己号牌上的数就是FDR校正后的p值啦。这个过程用数学语言表示为：

 

其中，q-valuei是校正后的p值，length(p)是总基因数，rank(p)是每个基因排队的号牌。

那么，为什么要用这个公式完成校正呢？且听小锐从原理说起。

 

寻找差异基因的漫漫长路，从假设开始……

我们在处理宏基因组差异基因的选择时，需要对两个样本的每个基因进行一次假设检验。如果我们有m个基因，那么我们就要做m次假设检验。每一次的假设检验的零假设H0为：两个样本的这个基因没有显著性差异。其中有m0个零假设是正确的，即这个基因在两个样本中确实没有显著性差异；但有m1=m-m0个零假设是错误的，即两个样本的这个基因是有显著性差异。m次检验之后，被拒绝的零假设的个数记为R。为了方便记忆，可用一张表格来表示假设检验的结果：

m次假设检验的结果

 

此时，我们需要多考虑一点：找到的“差异基因”里面，是否有找错了的呢？

为了解决这一问题，统计学家提出了控制错误率的标准。其中，Benjamini和Hochberg(1995)提出的FDR准则最为流行，它的标准相对来说更宽泛，在实际应用中能够获得更大的功效。也就是说，通过这种方式，我们可以在错误率很低的情况下找到真正有差异的基因。

根据“m次假设验证的结果”（记性不好的童鞋可以翻看前面的表格），我们可以得到FDR准则，即要求控制错误拒绝率。令Q=V / (V+S)，它表示被错误拒绝的零假设数目占所有被拒绝的零假设数目的比例。Q也是一个不可观测的随机变量。定义FDR如下：

 

当R=0时，定义Q=0

Benjamini和Hochberg提出了以FDR作为多重检验的准则，但是其检验的方法采用的是Simes(1986)提出的算法。

怎么修正？

重回开始的例子：在做差异基因的筛选时，有两个样本，m个基因。根据上述，我们需要做m次假设检验，会得到m个p值，分别为p1, p2,…, pm。将这m个p值由小到大排序，得到另外一组p(1)≤p(2)≤…≤p(m)，其中下标表示p值的排序位置。做这m次假设检验的控制错误率为

 

上式表示期望，也就是均值，即平均控制错误率。

假设我们按照每个P值小于0.05这个标准，我们一共可以挑出V+S=R个差异基因，这R个差异基因中有V个其实是没有差异表达的，S个是真正有差异表达的基因。如果有10,000个基因接受显著性检验，其中有90%，也就是9000个（即m0）满足原假设，那我们依然会从这些满足原假设的基因中找到450（即V）个有差异的，即使剩下1,000个（S）全部都被认为是差异基因，如果不利用FDR值进行修正，我们找到的差异基因中450/(450+1000)=31%（即Q=V/(V+S)）的差异基因是错误的。

为了使得这个错误率很低，我们要求平均错误控制率要小于一个预先设定的q值（比如0.05）即E(Q)≤q*

在统计学上，这也就等价于控制FDR不能超过5%。将E(Q)≤q*进行拆分：

 

因为V代表实际的差异基因，是固定不变的；要想使得控制错误率达到最小，就必须要求分母（V+S）最大化。

我们做了次假设检验，得到了V+S个差异基因（V+S个基因的P值≤0.05），V+S也表示排序后的p值小于0.05中的最大的那个p值对应的排序号（秩），即最大的那个p值为上述的pi，排序号为上述i。而此时V等于pi*m，所以E(Q)≤q*可以写成

 

修正的q*值取最小值，即为我们平时工作中用到的修正的FDR值q-valuei。这里的公式即为Bonferroni型多重检验过程中的公式。也是开始FDR的计算公式：

 

最后是FDR校正后的p值计算的一个小例子。

大家可以移步该网页查看

http://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html