做了多年的生物信息处理，和统计显著已经打过很多交道。然而限于生物出身背景，对于统计知识了解不是很深刻。细细想来，每次对外说出各种检验方式的名称（通常都是以人名命名的那种），自己却鲜有细细了解到底这些检验的原理以及具体操作过程。近日听了一场答辩，评委问道：“为何不做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 <- c(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 <- c(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的范围内，真阳性呈偏向分布，而假阳性呈均一分布：
最终能够多筛除假阳性，多保留真阳性：
当然，这种做法也会导致一部分原来p<0.05的真阳性基因被筛掉，也不能排除所有的假阳性。
更多参考：
1.Statquest: https://www.youtube.com/watch?v=K8LQSvtjcEo
2.Stack Overflow:https://stackoverflow.com/questions/10323817/r-unexpected-results-from-p-adjust-fdr
更多生信分析需求请加微信：13120220117
