总被审稿人提起的多重假设检验校正是什么？

xubzhu 2021-02-08

展开全文

生物信息学习的正确姿势

NGS系列文章包括NGS基础、在线绘图、转录组分析（Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析（ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析 (step-by-step)）、批次效应处理等内容。

单次检验的I类错误

假设检验是用于检验统计假设的一种方法，其基本思想是“小概率事件”原理，即小概率事件在一次试验中基本上不会发生。

假设检验的基本方法是提出一个空假设（null hypothesis），也叫做原假设或无效假设，符号是H0。一次检验有四种可能的结果，用下面的表格表示：

Type I error，I类错误，也叫做α错误，假阳性。
Type II error，II类错误，也叫做β错误，假阴性。

可以通过下面这张图形象的看到差异。

多次检验使得犯I类错误概率增大

在传统的假设检验中，单个检验的显著性水平或I型错误率 (错误拒绝原假设的概率)为计算出的P-value。但随着检验次数的增加，错误拒绝原假设的概率即I型错误率大大增加。

例如：如果我们进行了m次假设检验，至少有1个假阳性的概率是多少？

错误拒绝原假设的概率 P(Reject H0|H0=True) = α
决策正确的概率 P(No Reject H0|H0=True) = 1-α
P(在m次检验全部决策正确)=(1-α)^m
P(在m次检验中至少一次决策错误) = 1-(1-α)^m

随着检验次数的增多，出现至少一次决策错误的概率快速提高。当说起“根据假设检验的次数校正p值”时，意思是控制整体的I型错误率。

例如：当做差异基因检测时，每个基因分别进行检测生成一个p值。如果p值设置为0.05，每个差异基因识别出错的概率为5%。如果同时分析100个基因，按照p<0.05筛选的差异基因中有5个可能是差异不显著的。如果对一组10000个基因进行检测，按照p<0.05筛选的差异基因中有500个可能是差异不显著的。因此，同时进行多次统计检验时，校正每个基因的p值是很重要的。多重检验校正调整每个基因的p值，以使总体错误率小于或等于用户指定的p-cutoff value。

如何进行多重假设检验校正？

Family Wise Error Rate校正法控制假阳性率为0

Family Wise Error Rate是控制全部比较中至少出现一次Type I error的概率，也就是控制假阳性率为0。这是很严格的方式。

通常有两种计算方法：

Bonferroni correction方法

如果要维持整个检测 (做了m次检测)的Type I error rate < 0.05，则需要设定p-value为0.05/m作为筛选标准。反过来，如果我们做了10000次统计检测，采用Bonferroni correction方法校正后的p值就是原始P-value * 10000。

当然，我们也只是借这个方法理解校正的计算方式，实际却不用这个方法。

这对其中任何一个检测是否差异统计显著是不公平的，因为它取决于检测的总数目。一个检测放在有100次检测的操作集合中可能统计显著，而放在有1000次检测的操作集合中可能统计就不显著了，这是不合适的。

Bonferroni adjustments are, at best, unnecessary and, at worst, deleterious to sound statistical inference.
Perneger (1998)

Holm 校正方法

Holm 校正方法相对没有那么严苛。假设针对10000个基因进行了统计检验，对所有的原始P-value进行由小到大的排序分别为p1, p2, ..., p10000，校正后的p为：p1*10000, p2*9999, ..., p10000*1。

FDR校正法：允许一定的假阳性率

在实际应用中，我们希望减少Type I Error出现的可能，但也可以容许一定的假阳性率的存在。

Benjamini and Hochberg FDR (BH)

这是我们最常用的校正P-value控制假阳性率的方式。假设针对10000个基因进行了统计检验，对所有的原始P-value进行由小到大的排序分别为p1, p2, ..., p10000，校正后的FDR为：p1*10000/1, p2*10000/2, ..., p10000*10000/10000。与Bonferroni correction一致的地方是都乘以了检测总数，不一致的地方是BH算法在此基础上除去了各个原始p-value的排序值。

具体计算方式见下表（总检测次数为10次；控制FDR小于0.1）

Rank	P-value	FDR	FDR_formula	Reject H0	Reject_formula
1	0.0008	0.008	=B2*10/A2	TRUE	=C2<0.1
2	0.009	0.045	=B3*10/A3	TRUE	=C3<0.1
3	0.165	0.55	=B4*10/A4	FALSE	=C4<0.1
4	0.205	0.5125	=B5*10/A5	FALSE	=C5<0.1
5	0.396	0.792	=B6*10/A6	FALSE	=C6<0.1
6	0.45	0.75	=B7*10/A7	FALSE	=C7<0.1
7	0.641	0.915714286	=B8*10/A8	FALSE	=C8<0.1
8	0.781	0.97625	=B9*10/A9	FALSE	=C9<0.1
9	0.9	1	=B10*10/A10	FALSE	=C10<0.1
10	0.993	0.993	=B11*10/A11	FALSE	=C11<0.1

BH法有时也称fdr法，是我们最常用的多重假设检验校正方法，可以很好的控制假阳性率和维持统计检出力。R函数p.adjust可用来计算一组p-value校正后的fdr值。(DESeq2中返回的padj也是用BH方法控制的FDR)

q-value是什么？

q-value是Storey和Tibshirani提出的基于p-value分布的FDR计量方法，具体见什么，你算出的P-value看上去像齐天大圣变的庙？。

如何尽量减少统计检验次数

我们看到上面的校正方法多于统计检测次数有关，统计检测次数越多，校正也会越强烈。有没有合适的办法来规避一些无意义的统计检验呢？

WGCNA方法通过把基因聚类为模块再进行统计分析，大大降低了统计检验次数，具体见WGCNA分析，简单全面的最新教程
GSEA、GO等富集分析时合并相似的GO/KEGG通路再进行富集分析，如一文掌握GSEA，超详细教程中提到的合并共有基因数目超过70%的通路。
差异基因分析时过滤掉极低表达的基因 (低表达基因通常生物意义小或检测噪声大，即便有差异也难分清是生物差异还是技术差异)，如高通量数据中批次效应的鉴定和处理 - 系列总结和更新提到的方法。
DESeq2中还额外进行了independent filtering进行进一步过滤提高统计检出率。
没有通过过滤标准的基因校正后的padj赋值为NA (这也是之前总被问起的DESeq2结果中NA的来源)。