GWAS分析中的多重校正方法

根据一些资料,整理一下GWAS分析中可能用到的多重校正方法。

多重比较问题

首先我们看什么是多重比较问题,下面的内容来自于维基百科。

当统计分析涉及多个同时进行的统计测试时,就会出现多重比较,每个测试都有可能产生相同数据集或相关数据集的“发现”。规定的置信水平通常仅适用于单独考虑的每个测试,但通常希望对整个系列的同时测试都有一个置信水平。未能对多重比较进行补偿可能会产生重要的现实后果,如以下示例所示:

  • 假设治疗是一种新的写作教学方式,而控制是写作教学的标准方式。两组学生可以在语法、拼写、组织、内容等方面进行比较。随着更多的属性被比较,处理组和控制组由于随机抽样误差而在至少一个属性上出现差异的可能性越来越大。
  • 假设我们从减轻多种疾病症状中的任何一种来考虑药物的功效。随着更多症状被考虑,就至少一种症状而言,该药物似乎比现有药物有所改善的可能性越来越大。

在这两个示例中,随着比较次数的增加,被比较的组更有可能在至少一个属性方面出现差异。如果将结果作为涉及多重比较的分析的一部分而不是仅涉及单一比较的分析进行观察,我们对结果将推广到独立数据的信心通常会减弱。 例如,如果在 5% 的水平上执行一项测试并且相应的原假设为真,则错误拒绝原假设的可能性只有 5%。但是,如果进行了 100 次测试并且所有相应的零假设都为真,则错误拒绝(也称为误报或 I 类错误)的预期数量为 5。如果测试在统计上彼此独立,则至少一个错误的拒绝是 \(1-0.95^{100} = 99.4\%\) (当且仅当 \(P(AB)=P(A)P(B)\) ,事件 \(A\)\(B\) 独立)。 请注意,多重比较问题当然不会出现在经验检验多个假设的每种情况下,无论是顺序还是并行(发);粗略地说,每当在同一数据集(或不独立的数据集)上测试多个假设或在多个数据集中测试同一个假设时,就会出现多重比较问题。 多重比较问题也适用于置信区间。具有 95% 覆盖概率水平的单个置信区间将包含 95% 实验中的总体参数。但是,如果同时考虑 100 个置信区间,每个区间的覆盖概率为 95%,则非覆盖区间的预期数量为 5。如果区间在统计上彼此独立,则至少一个区间不包含总体的概率参数为 99.4%。 已经开发了一些技术来防止多重统计测试中出现的假阳性率和未覆盖率的膨胀。 - https://cn.wikiterhal.com/558787-multiple-testing-correction-IUCOPW

简单地理解一下,假设你想比较的两个分组实际上没有区别(即来自于同一个总体),但是只要你比较的属性/水平足够多,你总能在某个属性/水平找到所谓的显著性差别。但是,此时实际上造成这种“显著性差别”的原因是抽样误差

多重校正标准

FWER

族错误率 (Family-wise error rate, FWER) 指这多次检验中至少出现一次假阳性的概率,设 \(f\) 为出现假阳性的次数,则有 \[ FWER = P(f \geq 1) \]

FDR

错误发现率 (false discovery rate, FDR) 指在所有显著的检验中假阳性所占比例的期望,设 \(f\) 为出现假阳性的次数,\(d\) 为检验方法得到的显著性位点数目,则有 \[ FDR = E(f/d) \] 不同的多重校正方法都是通过控制整理的 FWER 或 FDR 低于给定水平,从而减少假阳性。可以证明,只要控制 \(\text{FWER} \leq \alpha\) ,就能保证 \(\text{FDR} \leq \alpha\) 。也就是说,基于 FWER 标准的方法比基于 FDR 标准的方法更加严格。

GWAS的多重校正

多重检验问题

GWAS 分析中,我们是逐个SNP进行检验,假设SNP数目为 \(n\) ,那么我们就做了 \(n\) 次检验。假设其中真正与性状有关联的 SNP 数目为 \(a \ll n\) , 每次检验显著性水平我们定为 \(\alpha\) ,那么在 \(n-a\) 个无关联位点的检验中 (原假设为真) 出现假阳性的次数为 \(\alpha \times (n-a)\) ,约等于 \(\alpha n\) 。比如假设总位点数目为 100 万,\(\alpha = 0.05\) ,那么出现假阳性的位点数目约为 5 万,很显然这是不可接受的。这就是GWAS研究中的多重检验问题。为了解决这一问题,提出了以下几种方法。

Bonferroni 校正法

首先,Bonferroni 校正是一种基于 FWER 标准的方法。假设真正关联的 SNP 的集合为 \(A\) ,所有 SNP 的集合为 \(B\) ,其中某个位点 \(i\) 出现假阳性的事件为 \(B_{i}\) ,基于事件集合,FWER 可以写为 \[ F W E R=P\left(\bigcup_{B_{i} \notin A} B_{i}\right) \] 根据 Bonferroni 不等式,我们有 \[ P\left(\bigcup_{B_{i} \notin A} B_{i}\right) \leqslant \sum_{B_{i} \notin A} P(B_{i}) \] 这个不等式可以根据 \(P(A \cup B)=P(A)+P(B)-P(A \cap B)\) 来证明。

可以得出,只要把单次检验的显著性水平设为 \(\alpha/n\) ,就可以使得整体检验的 \(\mathrm{FWER} \leq(\mathrm{n}-\mathrm{a}) \alpha / \mathrm{n} < \alpha\) 。因此该方法只需要将单次检验的显著性水平设为 \(\alpha/n\) ,就能控制整体检验出现假阳性的概率低于 \(\alpha\)

Bonferroni 校正是最严格的多重校正方法,由于不同的位点间存在连锁不平衡 (LD),即很多位点之间不是独立的,因此Bonferroni 校正可能过于保守,导致出现假阴性的概率增加。

置换检验法

置换检验法 (Permutation test) 也采用 FWER 标准,不同的是,它可以很好的处理不同检验之间存在的相关。置换检验法是Fisher于20世纪30年代提出的一种基于大量计算(computationally intensive),利用样本数据的全(或随机)排列,进行统计推断的方法,因其对总体分布自由,应用较为广泛,特别适用于总体分布未知的小样本资料,以及某些难以用常规方法分析资料的假设检验问题。

在GWAS分析中,由于不同的位点间存在连锁不平衡 (LD),使得等效的独立SNP数目低于位点总数 \(n\) ,置换检验法可以给出一个控制假阳性的合理的标准。其做法是将基因型和表型的对应关系随机打乱,我们认为随机打乱之后的基因型和表型之间不存在关联,也就是说此时所有的SNP均满足零假设。然后我们用打乱后的数据做GWAS,每次记录下最小的 P 值。重复此过程多次,比如1000次,我们就得到了原假设成立下最小 P 值的一个经验分布。如果我们将这些最小P值从小到大排序,如果我们选择显著性水平α=0.05,则将排序P值的0.05分位数作为校正后的显著性水平。这里用的最小 P 值也可换成最大统计量。

我们来直观解释一下其中的原理,FWER 标准是至少出现一次假阳性的概率,等价于零假设成立下至少有一个位点的 P 值小于校正后的显著性水平的概率,进而等价于零假设成立下最小 P 值小于校正后的显著性水平的概率。当我们将上面例子中 1000 次抽样排序P值的0.05分位数作为校正后的显著性水平时,此时 FWER ,即零假设成立下最小 P 值小于校正后的显著性水平的概率就是 0.05 ,满足要求。

置换检验法的每次随机打乱的过程中考虑到了 SNP 之间的连锁不平衡关系,因此相比于Bonferroni 校正可能更符合实际情况,其最主要的问题是计算量太大,短时间内几乎不可能实现。

FDR 的 BH 法

FDR 标准下最常用的方法就是 BH 法,其原理如下。

我们先将所有 SNPs 检验的 \(p\) 值从小到大排列, 假设我们找到的显著位点数目为 \(d\)\(p_{\mathrm{d}}\) 表示最后一个显著的 SNP 的 \(p\) 值 。设 FDR 下的显著水平为 \(\alpha\), 因为报告显著的 SNPs 有 \(\mathrm{d}\) 个, 所以按 FDR 标准, 允许错误发现的 SNPs 个数为 \(\alpha \times \mathrm{d}\) ,而此时理论上假阳性的位点数目约等于 \(\mathrm{n} \times p_{\mathrm{d}}\) (实际为 \(\mathrm{(n-a)} \times p_{\mathrm{d}}\))。为了把假阳性数目控制在允许错误发现数以内, 即 \(\mathrm{n} \times p_{\mathrm{d}} \leq \alpha \times \mathrm{~d}\), 等价于把 \(p\) 值校正为 \((\mathrm{n} / \mathrm{d}) \times p_{\mathrm{d}}\), 再和 \(\alpha\) 比较。 \(\mathrm{BH}\) 法就是确定 \(\mathrm{d}\) 的过程,具体步骤如下:

  1. \(\mathrm{n}\) 个 SNPs 的 \(p\) 值按从小到大排列;
  2. 把排在第 1 的 \(p\) 值乘以 \(\mathrm{n}\), 排在第 \(\mathrm{i}\)\(p\) 值 乘以 \(\mathrm{n} / \mathrm{i}\), 依次类推, 排在最后的 \(p\) 值乘以 \(\mathrm{n} / \mathrm{n}\), 得到校正后的 \(p\) 值 (有人也称为 \(q\) 值)
  3. 从排在最后的 SNP 开始, 按从后往前的顺序, 首次发现校正后 \(p\) 值小于 \(\alpha\) 的 SNP (设序号为 d)后, 报告序号为 1 到 \(\mathrm{d}\) 的 SNPs 显著。

我们可以简单证明一下,BH 法比Bonferroni 校正法要宽松,也就是上面说的基于 FDR 标准的方法比基于 FWER 标准的方法更加宽松。Bonferroni 校正法的判定标准是 \(p < \alpha/n\) ,即 \(np < \alpha\) ,我们可以理解为其校正后的 \(p\) 值为 \(np\) ,而 BH 法的校正后的 \(p\) 值为 \(np/d \leq \alpha\) ,因此 BH 法得到的校正后的 \(p\) 值更小,得到的显著位点数目只会更多,即该方法比Bonferroni 校正法宽松。

参考文献

  1. BENJAMINI, 1995, Controlling the False Discovery Rate: a Practical and Powerful Approach to Multiple Testing

  2. 黄杨岳,2013,全基因组关联研究中的多重校正方法比较

  3. [Permutation Test 置换检验(转)](https://www.cnblogs.com/Dzhouqi/p/3440902.html)

  4. GWAS中的多重假设检验

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2026 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信