例如,在我们对鉴定到的差异蛋白做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.adjustfunction (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}
其实该函数表达的意思是这样的:
联系客服