近年来,微阵列芯片技术(microarray technology)和高通量测序(high-throughput sequencing)技术的蓬勃发展为组学(omics)研究提供了有力的技术支撑,科学家们试图在全基因组、蛋白组、转录组、代谢组以及表观遗传学范围内探寻目标群体之间的差异[1-2]。这些高通量组学研究往往可以通过少量样本获得海量数据,而对大量高维数据的处理要求在给定的显著水平上进行假设检验,选择感兴趣的“显著”标记,并在后续研究中进行验证。由于验证通常只能在低通量下进行,因此如何准确的筛选出具有显著差异的候选标记至关重要。目前,单个假设检验理论现已成熟,即事先规定可能犯Ⅰ类错误(假阳性)的最大概率为α,在Ⅰ类错误不超过此概率的前提下减小犯Ⅱ类错误的概率[3]。然而,组学数据中存在数以千计的变量,从而面临同时进行多个假设检验的挑战。此时,传统的假设检验理论已不能满足此类数据分析处理的要求,增加的假阳性率会掩盖数据间的真实关联。其实,早在20世纪八九十年代便有学者提出常用的多组比较方法和多重检验的理论及过程;进入21世纪后,组学研究的发展和海量数据处理的需求使多重检验问题得到广泛关注,研究者开始在既往理论基础上寻求更好的方法来解决这一问题。本文主要讨论了多重假设检验问题的产生、处理方法及其在现实中应用和相关建议。
一、多重假设检验问题的提出在单个假设检验中,是否做出拒绝零假设(H0)主要取决于P值和预先给定的检验水准α,当P≤α时拒绝H0,反之则不能拒绝H0,其基本原则是在保证犯Ⅰ型错误的概率不超过α的情况下,寻找一种假设检验方法,使该检验犯Ⅱ型错误的概率β尽可能小。然而,当同时对多个假设进行检验时,情况将变得复杂。这时每个单独的假设检验都有其本身的Ⅰ型错误。如果不对其进行任何控制措施,犯Ⅰ型错误的概率会随着假设检验的次数的增加而迅速增加。一般地,若同时进行n次假设检验,每次检验水准为α,则不犯Ⅰ型错误的概率为1-α,n次均不犯Ⅰ型错误的概率为(1-α)n,此时总的检验水准为1- (1-α) n,其大小远远超过原本的α。
再者,单个假设检验规定的4种结果,即真阴性(H0成立不拒绝H0),假阳性(I型错误,H0成立拒绝H0)、假阴性(Ⅱ型错误,H0不成立不拒绝H0)和真阳性(H0不成立拒绝H0)的发生概率分别为U、V、T、S。而对于多重假设检验来说,若总共n个假设检验Hi(i=1, 2, 3…n),当每单个检验结果都落于单个检验规定的4种结果之一时,H0成立的个数为nt,H0不成立的个数为nf,n=nt+nf。接受H0的所有假设检验的个数为W,拒绝H0的所有假设检验的个数为R。此时,上述的U、V、T、S分别代表这n次假设检验中得到真阴性、假阳性(Ⅰ型错误)、假阴性(Ⅱ型错误)和真阳性结果的个数,而不是他们各自发生的概率(1-α、α、β和1-β)[4]。在当前的组学研究中,每次实验都要进行几十个或数百个甚至更多的假设检验,实际发现的假阳性结果数目远远超过预期的真实数量,这给后续的验证性研究增加了极大的困难和不确定性。因此,多重假设检验的核心问题是如何使用一种合理的错误测度来控制错误拒绝H0(Ⅰ型错误)的个数V或比率V/R,进而寻找恰当的检验方法尽可能提高检验功效。
二、多重检验的错误测度及检验方法目前,常用的错误测度有以下几种:平均比较错误率(per-comparison error rate,PCER)、平均族错误率(per-family error rate,PFER)、总体/族错误率(family-wise error rate,FWER)、错误发现率(false discovery rate,FDR)和阳性错误发现率(positive false discovery rate,pFDR)。各种与之相应的多重假设检验的控制过程,即用于控制多重假设检验错误率的新的检验方法,也应运而生。这些不同的过程主要是对分析过程中得到的P值进行调整。这里主要介绍传统的错误测度FWER和当前受到广泛关注的错误测度FDR及其控制过程。
1、错误测度FWERFWER指在多重假设检验中至少有一个检验犯Ⅰ型错误的概率,FWER=Pr(V≥1)[5]。FWER错误测度应用于多重假设检验问题的初期,是最为传统的方法。FWER的控制过程主要有传统Bonferroni程序,以及在此基础上的一些改进程序(Holm、Hochberg、Hommel程序等),研究者必须根据自己的数据类型选择适合的控制程序。此外,值得注意的是,该测度在实际应用中的关键问题是如何定义总体/族(family)。例如研究8个剂量点的辐射对10个指标的影响,那么单个剂量点的10次检验是否为一个family或是将在8个剂量点进行的共80次检验作为一个family,若后期又补充了2个剂量点,那么是否要与前期8个剂量点合并后重新进行检验合计100次检验作为一个family。实际中对此没有严格的规定,研究者需要根据假阳性的出现对该研究的不利程度,在检验开始之前做出决定,以避免不合理的调整family中检验数量来获得预期结果。
(1) 传统Bonferroni程序:传统的Bonferroni校正比较简单保守,它既不考虑检验统计量的顺序,也不考虑原始P值的顺序,平等的对待所有假设。具体方法如下:
① 假定同时检验n个独立假设,给定总体检验水平α,则设置截断点
② 令每一个假设对应的原始P值Pi与S相比较,即对于每一个独立假设而言,它们各自的临界值
③ 若第i个检验的P值Pi≤S,则拒绝原假设H0i(i=1, 2, 3…n)。
举例说明,原本使用P=0.05为标准判断是否存在差异,那么对于1 000个实际没有差异的基因也会得出“50个基因差异表达”的结果,而将拒绝阈提高到0.05/1 000后,同样1 000次比较得到的假阳性次数依然被控制在0.05。
传统Bonferroni过程能够强控制FWER在α水平,且计算简单方便,但这种方法过度调整了P值,可能导致过高的假阴性率,特别是在数据高度相关的情况下。
(2) 控制FWER的自上而下(Step-down)和自下而上(Step-up)过程:为了克服传统Bonferroni程序过于保守的缺点,Holm[6]对其进行改进,提出了一种基于FWER测度的Step-down过程,称为Holm程序。该程序考虑了原始P值的顺序,是一种Step-up过程。具体如下:
① 给定总体检验水平α,依次将n个假设检验所对应的原始P值Pi(i=1, 2, 3…n)按从小到大(显著性逐渐减弱)的顺序排列,得到P1≤P2≤P3≤…≤Pn。
② 取截断点
③ 按从小到大的顺序逐个将Pi与截断点Si相比较,若Pi≤Si,则拒绝Pi所对应的原假设H0i。当发现第1个不被拒绝的原假设时,程序终止,后续所有假设失效。
从以上过程可以看出当i≥ 2时,
Hochberg[7]也提出了另一种用于校正传统Bonferroni程序的Step-up过程,即Hochberg程序。它与Holm程序类似,但检验顺序恰好相反,从最不显著的P值开始倒序检验,将Pi与截断点Si相比较,当出现第一个被拒绝的原假设时,程序停止,后续所有检验的零假设均被拒绝。除此之外,用于FWER控制的Step-up过程还有Hommel程序[8]和Rom[9]过程,这两个程序的功效比Hochberg程序高,但与它们定义式和计算的复杂程度相比,提高的功效似乎并没有那么有吸引力。
尽管目前已经提出了多种用于控制FWER的程序,但整体而言,FWER作为多重假设检验错误测度仍然过于保守,检验功效偏低,会导致更高的假阴性率。而且,传统的FWER并不特别适用于当前的大规模组学研究:一是FWER控制方法通常假定检验统计量服从多元正态或t分布,但现实数据往往是各种分布混合的情况;二是组学研究是一种探索性的分析策略,以期为后续的验证性研究提供候选基因,这类数据的n往往成千上万,此时需要显著性P足够小才能得到一个差异具有统计学意义的指标,此时控制FWER显然没有多大意义。
2、FDR错误测度Benjamini和Hochberg[10]在1995年首次提出了FDR的概念,并给出了在多重假设检验中对它的控制方法。尽管这在当时并未受到重视,但随着微阵列检测和高通量测序技术的应用,以及组学研究的开展,高维、海量数据的出现使FDR得到日益广泛的应用。组学研究的探索性决定了它需要在假阳性个数和拒绝原假设的总个数之前寻求平衡,即在检出更多候选基因的同时将假阳性率控制在可接受的范围,而FDR能有效地满足这一需求。与FWER相比,FDR具有以下优点:①取值灵活,作为假设检验错误率的控制指标,其取值可以根据需要进行调整,而传统FWER的取值则较为固定,通常定为0. 05。②FDR可以作为筛选出的差异变量的评价指标,而FWER主要用于控制I型错误。③保守度低,检验功效更高,当nt=n时,FDR=FWER;但当nt<n时,控制FWER必然控制FDR,但反之则不然。目前,研究者提出了许多FDR的控制过程,其中最经典的本杰明-霍克伯格(Benjamini-Hochberg,BH)过程。
(1) BH过程:Benjamini和Hochberg[10]在提出FDR概念的同时提出此程序,它实际上是对Simes(1986)[5]提出的FWER控制程序进行了改进,使其在控制FDR上具有较高的功效。该程序是一个Step-up过程,原理如下:
① 给定FDR值Q,n个原假设H0i(i=1, 2, 3,…n)所对应的原始P值Pi按从小到大的顺序排列,得到P1≤P2≤P3≤…≤Pn。
② 取截断点
③ 从i=n开始,倒序将Pi与Si比较。当出现第一个被拒绝的原假设时,程序停止,后续所有检验的零假设均被拒绝,即使那些假设检验得到的原始P值并不小于截断点Si。
一般而言,在文章中最好能够提供原始P值,并使用基于FDR控制的BH过程发现有意义的指标,但若学科领域内的文章普遍使用经BH法调整后的P值时,则需注意调整后的P值
(2)其他FDR控制过程:1999年Yekutieli和Benjamini[11]对BH过程进行了改进,提出BY过程,原理是利用重复抽样的方法计算P值,可以在变量相关条件下控制FDR在给定的Q水平,但相对保守;同年Benjamini和Liu[12]则提出了一种Step-down的控制方法,即BL过程,该过程与经典的BH法类似,只是进行P值比较的顺序不同。Benjamini和Hochberg[13]将对真实原假设nt个数或比例的估计糅合进原来的BH过程。随后,Benjamini和Yekutieli[14]对进一步提出了可用于不同变量间独立和相关条件下的FDR的控制过程。Benjamini等[15]又提出了一种两阶段线性Step-up过程,可以在不同的显著水平下两次使用上述介绍的基础过程,特别是在变量相关条件下得到的FDR估计较为稳健。
3、pFDR错误测度Storey[16]认为在现实的组学或微阵列分析数据中,通常可以检出有差别的“显著”标记,即几乎肯定有Pr(R>0)→1,此时FDR可表达为
目前,对pFDR控制过程的研究大多在Bayesian框架下进行[17]。Storey[18]提出q-值法来估计pFDR,其基本思想是在最初的FDR控制过程基础上,估计了无差异变量占总变量个数的比例,以提高原有方法的检验效能。多重检验统计量在独立同分布条件下,pFDR可用一个后验概率的形式表示,pFDR随着Ⅰ型错误的增大而增大,随着检验效能的增大而减小。与单个假设检验中控制Ⅰ型错误α的P值类似,q值是为了在多重检验中控制pFDR,实质上就是一个后验Bayesian P值,具体定义为拒绝检验量的最小的pFDR值的大小,当q≤α时可保证pFDR≤α。然而,q-值法要求无效假设下的P值服从均匀分布,但现实中数据的复杂性往往无法满足这一条件,为此,有些方法直接使用样本统计量的经验分布对P值进行估计,如SAM法直接使用d统计量,通过重复抽样得到无效假设的统计量分布[19]。另外,一些经验Bayesian法,例如Locfdr法、Fdrtool法和PRML法等也在一定程度上解决了这一问题[20-22]。值得注意的是,这些pFDR估计方法并非适用于所有类型的数据,如何合理选择无效假设,构建更合适的统计量或得到准确的P值,需要进一步的研究。
三、60Co γ射线事故受照者外周血白细胞的基因表达谱实例分析为说明多重检验在高维海量数据分析中的必要性,采用Chi等[23]公开发布的基因表达数据,该研究报道了60Co γ射线意外辐照者外周血白细胞的基因表达谱,纳入3名受照者和3名健康对照,3名受照者采血时间分别是受照后1、2和3年。现以全部受照者与健康对照的部分基因表达差异比较为例,统计分析使用R4.0.2软件,借助Bioconductor的GEOquery包、Limma包、fdrtool包等R语言包和p.adjusted命令来实现。Stata软件中的Michael Andersons代码也可实现FDR控制。按照α=0.05,若不考虑多重假设检验的问题,不做任何调整地计算原始P值(Pi)可得到115个阳性结果。经4种FWER控制过程(Bonferroni法、Holm法、Hochberg法、Hommel法)校正后的adj-Pi均>0.05,没有发现阳性结果,而经基于控制FDR的BH法矫正后可得到5个阳性结果,见表 1。在此,简述以FWER为错误测度的Bonferroni校正和以FDR为错误测度的BH校正过程。在检验水准α=0.05情况下,Bonferroni过程中单个检验的临界值为0.05/358=0.000 13,因此只有单个检验的P < 0.000 13时才具显著性,而本例中最小的原始P值为0.000 15,因此,在该标准下不能发现阳性结果。
![]() |
表 1 FWER和FDR控制过程校正后的组学数据差异[23] Table 1 Difference comparison of omics data using the correction process based on FWER and FDR error measures |
BH法是一种Step-up过程,如表 1中第3列所示,将原始P值从小到大排序,最小P值的秩次为i=1。此时不再控制单个检验水准,而是将错误发现率Q控制在0.05。从i=n开始,倒序将每个原始P值Pi与其所对应的BH法临界值(i/n)Q(i为原始P值的秩次,n为总检验次数,Q为规定的错误发现率)相比较,直到出现第一个被拒绝的原假设。本例中经过一系列计算,可发现Gene108所对应的原始P值P6(0.002 2)超过其BH法临界值0.000 8(6/358×0.05),而Gene206所对应的原始P值P5(0.000 51)小于其BH法临界值0.000 6(5/358×0.05),因此当i=5时,第1次拒绝原假设,BH程序终止,认为前5个检验是显著的。由BH法得到的调整后P值称为q值,可以由公式
事实上,多重假设检验最初的研究方向是的多组比较,特指通过方差分析(计量资料且满足独立、正态性和方差齐性)、χ2检验(计数资料)或秩和检验(计量资料且不满足独立、正态性和方差齐性),得到多组总体参数或总体分布存在显著差异的结论后,进一步进行组间两两比较来回答具体哪两个组间有差别,包括样本均数、样本率和秩均值的多重比较[3]。而广义的多重假设检验一般指通过对多变量的反复检验来回答同一问题,如多元回归中各自变量的假设检验,微阵列数据或组学中对各组间差异表达蛋白/基因的筛选[24]。在实际应用中,研究者需根据具体的研究目的、研究设计和数据性质决定是否进行多重检验,并选择恰当的校正方法。
尽管对多重假设检验进行校正可以降低犯Ⅰ型错误的概率,但会相应提高犯Ⅱ型错误的概率,因此研究者根据研究具体情况进行代价利益评估以确定是否要进行校正。一般来说,如若不想对多重假设检验进行任何校正,可采用Bayes方法分析数据,查看均数差异的等效区间是否超过95%最高后验概率密度可信区间;也可在报告中给出研究所进行的所有假设检验的次数,例如共进行40个检验,仅报告10个检验,其中7个差异有统计学意义;还有一些学者认为在计划比较情况下不应对多重假设检验进行任何校正,但计划比较这一术语的定义当前并不明确。因此,建议通过控制FWER或FDR对多重假设检验进行校正,尤其当面对微阵列数据或组学数据这类变量高度相关并涉及大规模假设检的资料时,应采用控制FDR或pFDR的相关过程,无论是采用经典的本杰明-霍克伯格(BH)法,还是改进后的本杰明-耶库蒂利(BY)法、BL法或经验Bayes方法,使发现差异同时能将假阳性率控制在可接受的范围。然而,在验证性分析时建议采用FWER控制过程(通常采用Bonferroni校正)以尽可能降低假阳性率,发现真实差异。
利益冲突 无
作者贡献声明 高宇负责调研文献,撰写文章;苏垠平负责修改论文;孙全富提供写作思路,指导写作
[1] |
Lappalainen T, Scott AJ, Brandt M, et al. Genomic analysis in the age of human genome sequencing[J]. Cell, 2019, 177(1): 70-84. DOI:10.1016/j.cell.2019.02.032 |
[2] |
Uhlén M, Fagerberg L, Hallström BM, et al. Tissue-based map of the human proteome[J]. Science, 2015, 347(6220): 1260419. DOI:10.1126/science.1260419 |
[3] |
孙振球. 医学统计学[M]. 北京: 人民卫生出版社, 2018. Sun ZQ. Medical statistics[M]. Peking: People's Medical Publishing House, 2018. |
[4] |
刘晋, 张涛, 李康. 多重假设检验中FDR的控制与估计方法[J]. 中国卫生统计, 2012, 29(2): 305-308. Liu J, Zhang T, Li K. The control and estimation method of FDR in multiple hypothesis testing[J]. Chin J Health Stat, 2012, 29(2): 305-308. DOI:10.3969/j.issn.1002-3674.2012.02.054 |
[5] |
Simes J. An improved Bonferroni procedure for multiple tests of significance[J]. Biometrika, 1986, 73: 751-754. DOI:10.1093/biomet/73.3.751 |
[6] |
Holm S. A simple sequentially rejective Bonferroni test procedure[J]. Scand J Stat, 1979, 6(2): 65-70. DOI:10.1093/biomet/75.4.800 |
[7] |
Hochberg Y. A sharper Bonferroni procedure for multiple tests of significance[J]. Biometrika, 1988, 75(4): 800-802. |
[8] |
Homeel G. A stagewise rejective multiple test procedure based on a modified Bonferronitest[J]. Biometrika, 1988, 75(2): 383-338. DOI:10.1093/biomet/75.2.383 |
[9] |
Rom DM. A sequentially rejective procedure based on a modified Bonferroni inequality[J]. Biometrika, 1990, 77(3): 663-665. DOI:10.1093/biomet/77.3.663 |
[10] |
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing[J]. J R Stat Socser B, 1995, 57(1): 289-300. DOI:10.1111/i.2517-6161.1995.tb02031.x |
[11] |
Yekutieli D, Benjamini Y. Resampling-based false discovery rate controlling multiple test procedure for correlated test statistics[J]. J Stat Plann Inference, 1999, 82(1-2): 171-196. DOI:10.1016/S0378-3758(99)00041-5 |
[12] |
Benjamini Y, Liu W. A step-down multiple testing procedure that controls the false discovery rate under independence[J]. J Stat Plann Inference, 1999, 82(1-2): 163-170. DOI:10.1016/S0378-3758(99)00040-3 |
[13] |
Benjamini Y, Hochberg Y. On the adaptive control of the false discovery rate in multiple testing with independent statistics[J]. J Educ Behav Stat, 2000, 25(1): 60-83. DOI:10.3102/10769986025001060 |
[14] |
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency[J]. Ann Stat, 2001, 29(4): 1165-1188. DOI:10.1214/aos/1013699998 |
[15] |
BenjaminiY, Krieger AM, Yekutieli D. Adaptive linear step-up procedures that control the false discovery rate[J]. Biometrika, 2006, 93(3): 491-507. DOI:10.1093/biomet/93.3.491 |
[16] |
Storey JD. A direct approach to false discovery rates[J]. J R Stat Soc B, 2002, 64(3): 479-498. DOI:10.1111/1467-9868.00346 |
[17] |
Farcomeni A. A review of modern multiple hypothesis testing, with particular attention to the false discovery proportion[J]. Stat Methods Med Res, 2008, 17(4): 347-388. DOI:10.1177/0962280206079046 |
[18] |
Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value[J]. Ann Stat, 2003, 31(6): 2013-2035. DOI:10.1214/aos/1074290335 |
[19] |
Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response[J]. Proc Natl Acad Sci USA, 2001, 98(9): 5116-5121. DOI:10.1073/pnas.091062498 |
[20] |
Bradley E. Large-scale simultaneous hypothesis testing[J]. J Amer Stat Assoc, 2004, 99(465): 96-104. DOI:10.1198/016214504000000089 |
[21] |
Strimmer K. A unified approach to false discovery rate estimation[J]. BMC Bioinformatics, 2008, 9: 303. DOI:10.1186/1471-2105-9-303 |
[22] |
Martin R, Tokdar ST. A nonparametric empirical Bayes framework for large-scale multiple testing[J]. Biostatistics, 2012, 13(3): 427-439. DOI:10.1093/biostatistics/kxr039 |
[23] |
Chi C, Tian R, Liu H, et al. Follow-up study of abnormal biological indicators and gene expression in the peripheral blood of three accidentally exposed persons[J]. J Radiat Res, 2013, 54(5): 840-851. DOI:10.1093/jrr/rrt022 |
[24] |
荀鹏程, 赵杨, 柏建岭, 等. 微阵列数据的多重比较[J]. 中国卫生统计, 2006, 23(1): 5-8. Xun PC, Zhao Y, Bai JL, et al. Multiple comparisons of microarray data[J]. Chin J Health Stat, 2006, 23(1): 5-8. DOI:10.3969/j.issn.1002-3674.2006.01.002 |