【统计】对P值进行FDR校正的原理（附代码演示）

做了多年的生物信息处理，和统计显著已经打过很多交道。然而限于生物出身背景，对于统计知识了解不是很深刻。细细想来，每次对外说出各种检验方式的名称（通常都是以人名命名的那种），自己却鲜有细细了解到底这些检验的原理以及具体操作过程。近日听了一场答辩，评委问道：“为何不做p值校正？”。闻后心里一惊，因为自己并不清楚校正的原理，如果自己遇到这种问题似乎会很难堪了。

结束后，查阅了多方资料。笔者认为对FDR的数学意义及计算实现已经有了正确的认识。总结一下相关知识，在此以飨读者。

---

利用Benjamini–Hochberg方法计算FDR的计算及R语言实现

FDR的计算相当简单，包括以下几步：

  1.对p值进行从小到大的排序，标记上序号1~n；

  2.其中，最大的FDR(不考虑重复则为第n位)等于最大的p值；

  3.对于n-1位的FDR，取下面两者的较小值：

• 上一步(第n位)计算得出的FDR值；

• p-value*n/(n-1)

  4.不断迭代第三步(n-2,n-3....)，直至计算到最小p值对应的FDR。

下面直接在R里实现：

###例如这10个p值进行校正

temp <->0.01,0.11,0.21,0.31,0.41,0.51,0.61,0.71,0.81,0.91)p.adjust(temp, method='fdr')###[1] 0.1000000 0.5500000 0.7000000 0.7750000 0.8200000 0.8500000 0.8714286 0.8875000 0.9000000 0.9100000

根据上述定义，我们知道以此计算出的FDR是有可能出现相同值的---即使其原始p值不同。在上面的例子中，第9位的p值为0.81，根据公式计算出的待选FDR值为0.81*10/9=0.9，比第10位的0.91要小，因此成为真正的FDR。如果我们更改一下条件，将0.81改成0.90，那么结果就会有所不同：

###将0.81改为0.9

temp <->0.01,0.11,0.21,0.31,0.41,0.51,0.61,0.71,0.9,0.91)p.adjust(temp, method='fdr')###[1] 0.1000000 0.5500000 0.7000000 0.7750000 0.8200000 0.8500000 0.8714286 0.8875000 0.9100000 0.9100000###出现了相同的FDR值

因此这个现象是很正常的。笔者看到Stack Overflow上还有类似的提问，感兴趣的童鞋可进一步参看：

Benjamini–Hochberg方法的几何解释

我们已经知道，p值实际上是拒绝0假设的置信度。这也就意味着，对于n个样本，我们就会期望产生n*p个错误判断。换而言之，对于实际上没有区别的分布进行随机抽样，我们计算出来的p值应该是均一分布的。

注：也就是说有5%的可能性抽出假阳性出来。

如果我们从两个有差别的分布中抽样，那么p值的分布将富集在左侧，而这种趋势与差异程度显然是正相关的：

而对于生物样本分析而言(例如，一个转录谱分析)，我们所看到的常常是以上两种效果的叠加。并且值得注意的是，差异项通常远少于非差异项(例如，对于RNA-Seq而言，通常只有总基因数10%以下的差异基因)：

因此这种叠加的结果是更偏向无差异的分布，可能长得是这样的：

注：也就是说，如图所示的蓝色和红色部分，如果以p=0.2为threshold(当然这只是个例子)，红色的部分都是假阳性。

对于p=0.05的threhold而言，在这里我们有近1:1的真假阳性比例，也就是说这种情况下我们找到一半的差异基因都是假的：

Benjamini–Hochberg方法想解决这一问题的方法的本质，无非就是提高p-value筛选的阈值，因为在p-value<0.05的范围内，真阳性呈偏向分布，而假阳性呈均一分布：>

最终能够多筛除假阳性，多保留真阳性：