浅谈多重检验校正FDR

例如，在我们对鉴定到的差异蛋白做GO功能注释后，通常会计算一个p值。当某个蛋白的p值小于0.05（5%）时，我们通常认为这个蛋白在两个样本中的表达是有差异的。但是仍旧有5%的概率，这个蛋白并不是差异蛋白。那么我们就错误地否认了原假设（在两个样本中没有差异表达），导致了假阳性的产生（犯错的概率为5%）。如果检验一次，犯错的概率是5%；检测10000次，犯错的次数就是500次，即额外多出了500次差异的结论（即使实际没有差异）。为了控制假阳性的次数，于是我们需要对p值进行多重检验校正，提高阈值。

第一种方法Bonferroni，最简单严厉的方法。
例如，如果检验1000次，我们就讲阈值设定为5% / 1000 = 0.00005；即使检验1000次，犯错误的概率还是保持在N×1000 = 5%。最终使得预期犯错误的次数不到1次，抹杀了一切假阳性的概率。但是该方法虽然简单，但是检验过于严格，导致最后找不到显著表达的蛋白（假阴性）。

第二种方法FDR（False Discovery Rate）
相对Bonferroni来说，FDR用比较温和的方法对p值进行了校正。其试图在假阳性和假阴性间达到平衡，将假/真阳性比例控制到一定范围之内。例如，如果检验1000次，我们设定的阈值为0.05（5%），那么无论我们得到多少个差异蛋白，这些差异蛋白出现假阳性的概率保持在5%之内，这就叫FDR＜5%。

那么我们怎么从p value 来估算FDR呢，人们设计了几种不同的估算模型。其中使用最多的是Benjamini and Hochberg方法，简称BH法。虽然这个估算公式并不够完美，但是也能解决大部分的问题，主要还是简单好用！

FDR的计算方法
除了可以使用excel的BH计算方法外，对于较大的数据，我们推荐使用R命令p.adjust。 p.adjust(p, method = p.adjust.methods, n = length(p))

p.adjust.methods
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
#   "fdr", "none")

我们还可以从R命令p.adjust的源代码，了解其运行的机制是什么。

> p.adjust
function (p, method = p.adjust.methods, n = length(p)){
    method <- match.arg(method)
    if (method == "fdr") 
    method <- "BH"
    nm <- names(p)
    p <- as.numeric(p)
    ……
    BH = {
        i <- lp:1L
        o <- order(p, decreasing = TRUE)
        ro <- order(o)
        pmin(1, cummin(n/i * p[o]))[ro]
    }
    ……
     p0
}

其实该函数表达的意思是这样的：

我们将一系列p值、校正方法（BH）以及所有p值的个数（length(p)）输入到p.adjust函数中。
将一系列的p值按照从大到小排序，然后利用下述公式计算每个p值所对应的FDR值。公式：p * (n/i)， p是这一次检验的p value，n是检验的次数，i是排序后的位置ID（如最大的P值的i值肯定为1，第二大则是2，依次至最小为n）。
将计算出来的FDR值赋予给排序后的p值，如果某一个p值所对应的FDR值大于前一位p值（排序的前一位）所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值。因此会产生连续相同FDR值的现象；反之则保留计算的FDR值。
将FDR值按照最初始的p值的顺序进行重新排序，返回结果。

参考：http://www.omicshare.com/forum/thread-173-1-1.html

本文出自于http://www.bioinfo-scrounger.com转载请注明出处