校正测试相关的多个测试的 p 值(遗传学)

机器算法验证 相关性 多重比较 统计学意义 遗传学
2022-01-30 05:14:15

我从很多测试中获得了 p 值,并且想知道在对多次测试进行校正后是否真的有一些重要的东西。并发症:我的测试不是独立的。我正在考虑的方法(Fisher 乘积方法的一种变体,Zaykin 等人,Genet Epidemiol,2002)需要 p 值之间的相关性。

为了估计这种相关性,我目前正在考虑引导案例,运行分析并将结果向量的 p 值相关联。有没有人有更好的主意?甚至对我的原始问题有更好的想法(纠正相关测试中的多个测试)?

背景:我在逻辑上回归我的受试者是否患有特定疾病,取决于他们的基因型(AA、Aa 或 aa)与协变量之间的相互作用。然而,基因型实际上是很多(30-250)个单核苷酸多态性(SNP),它们肯定不是独立的,而是处于连锁不平衡状态。

4个回答

这实际上是全基因组分析研究(GWAS)中的热门话题!我不确定您正在考虑的方法在这种情况下是否最合适。一些作者描述了 p 值的汇集,但在不同的背景下(复制研究或荟萃分析,参见例如 (1) 以获得最近的评论)。当人们想要为给定基因导出唯一的 p 值时,通常使用通过 Fisher 方法组合 SNP p 值。这允许在基因水平上工作,并减少后续测试的维数,但正如你所说,标记之间的非独立性(由空间共位或连锁不平衡引起,LD)引入了偏差。更强大的替代方案依赖于重新采样程序,

我对引导(带替换)的主要担忧是您正在引入一种人为的关联形式,或者换句话说,您创建了虚拟双胞胎,从而改变了 Hardy-Weinberg 平衡(以及最小等位基因频率和调用率)。对于排列单个标签并保持基因分型数据不变的排列方法,情况并非如此。通常,plink软件可以为您提供原始和置换的 p 值,尽管它(默认情况下)使用带有滑动窗口的自适应测试策略,如果 SNP 看起来低于考虑不是“有趣的”;它还具有计算 maxT 的选项,请参阅在线帮助

但鉴于您正在考虑的 SNP 数量较少,我建议依赖在multtest R 包中实现的基于 FDR 或 maxT 测试(请参阅 参考资料mt.maxT),但基因组应用重采样策略的权威指南是Multiple Testing Procedures with Applications to基因组学,来自 Dudoit & van der Laan (Springer, 2008)。另请参阅 Andrea Foulkes 的R 遗传学书,该书在 JSS 中进行了评论。她在多种测试程序方面拥有丰富的资料。

进一步说明

许多作者指出,简单的多重测试校正方法(例如 Bonferroni 或 Sidak)对于调整单个 SNP 的结果来说过于严格。此外,这些方法都没有考虑到 SNP 之间存在的相关性,因为 LD 标记了跨基因区域的遗传变异。已经提出了其他替代方案,例如用于多重比较的 Holm 方法的导数 (3)、隐马尔可夫模型 (4)、条件或正 FDR (5) 或其导数 (6),仅举几例。所谓的差距统计或滑动窗口在某些情况下已被证明是成功的,但您会在 (7) 和 (8) 中找到一个很好的回顾。

我还听说过有效利用单倍型结构或 LD 的方法,例如 (9),但我从未使用过它们。然而,它们似乎更多地与估计标记之间的相关性有关,而不是您所说的 p 值。但实际上,您可能更好地考虑连续测试统计之间的依赖结构,而不是相关 p 值之间的依赖结构。

参考

  1. Cantor, RM, Lange, K 和 Sinsheimer, JS。优先考虑 GWAS 结果:统计方法及其应用建议的回顾Am J Hum Genet。2010 86(1):6-22。
  2. Corley, RP, Zeiger, JS, Crowley, T 等人。候选基因与青少年反社会药物依赖的关联药物和酒精依赖 2008 96:90–98。
  3. Dalmasso、C、Génin、E 和 Trégouet DA。全基因组关联研究中等位基因频率的加权霍尔姆程序遗传学 2008 180(1):697–702。
  4. Wei, Z, Sun, W, Wang, K 和 Hakonarson, H.通过隐马尔可夫模型进行全基因组关联研究中的多重检验生物信息学 2009 25(21): 2802-2808。
  5. Broberg, P.对未改变基因比例和错误发现率估计值的比较回顾BMC 生物信息学 2005 6:199。
  6. 需要、AC、Ge、D、Weale、ME 等。精神分裂症中 SNP 和 CNV 的全基因组研究公共科学图书馆基因。2009 年 5(2):e1000373。
  7. Han, B, Kang, HM 和 Eskin, E.数百万相关标记的快速准确的多重测试校正和功率估计PLoS 遗传学 2009
  8. Liang, Y 和 Kelemen, A.分析复杂疾病基因组研究中相关高维 snp 数据的统计进展和挑战2008 年统计调查 2:43–60。- 有史以来最好的近期评论
  9. 尼霍尔特博士 相互连锁不平衡中单核苷酸多态性多重检测的简单修正Am J Hum Genet。2004 74(4):765–769。
  10. Nicodemus,KK,Liu,W,Chase,GA,Tsai,YY 和 Fallin,MD。使用主成分与单倍型阻断算法比较大型单核苷酸多态性研究中多次测试校正的 I 型错误BMC 遗传学 2005;6(增补 1):S78。
  11. Peng, Q, Zhao, J 和 Xue, F.基于 PCA 的自举置信区间检验涉及多个 SNP 的基因-疾病关联BMC 遗传学 2010, 11:6
  12. Li, M, Romero, R, Fu, WJ 和 Cui, Y (2010)。用自适应套索映射单倍型-单倍型相互作用BMC Genetics 2010, 11:79 -- 虽然与问题没有直接关系,但它涵盖了基于单倍型的分析/上位效应

使用像 bonferroni 这样的方法很好,但问题是如果你有很多测试,你不可能找到很多“发现”。

您可以使用 FDR 方法进行相关测试(有关详细信息,请参见此处),问题是我不确定您是否可以提前说出您的相关性是否都是正相关。

在 R 中,您可以使用 p.adjust 进行简单的 FDR。对于更复杂的事情,我会看一下multcomp,但我没有通过它来查看依赖关系的解决方案。

祝你好运。

我认为多元正态模型正在用于对相关的 p 值进行建模并获得正确类型的多重测试校正。 对数百万个相关标记进行快速准确的多重测试校正和功率估计。PLoS Genet 2009讨论了它们并提供了其他参考资料。这听起来与您所说的相似,但我认为除了获得更准确的全局 p 值校正之外,还应该使用 LD 结构知识来消除与因果标记相关的标记产生的虚假阳性。

我正在为完全相同的问题寻找可行的解决方案。我发现的最好的是Foulkes Andrea 在他的书Applied Statistical Genetics with R(2009)中介绍的Null Unrestricted Bootstrap与他专门考虑回归的所有其他文章和书籍相反。除了其他方法外,他还建议使用 Null Unrestricted Bootstrap,该方法适用于无法轻松计算残差的情况(如在我的情况下,我对许多独立回归(基本上是简单的相关性)进行建模,每个回归都具有相同的响应变量和不同的片段)。我发现这个方法也被称为maxT方法。

> attach(fms)
> Actn3Bin <- > data.frame(actn3_r577x!="TT",actn3_rs540874!="AA",actn3_rs1815739!="TT",actn3_1671064!="GG")
> Mod <- summary(lm(NDRM.CH~.,data=Actn3Bin))
> CoefObs <- as.vector(Mod$coefficients[-1,1]) 
> B <-1000
> TestStatBoot <- matrix(nrow=B,ncol=NSnps)
> for (i in 1:B){
+    SampID <- sample(1:Nobs,size=Nobs, replace=T)
+    Ynew <- NDRM.CH[!MissDat][SampID]
+    Xnew <- Actn3BinC[SampID,]
+    CoefBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,1]
+    SEBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,2]
+    if (length(CoefBoot)==length(CoefObs)){
+       TestStatBoot[i,] <- (CoefBoot-CoefObs)/SEBoot
+    }
+ }

一旦我们有了所有的TestStatBoot矩阵(在行中我们有引导复制,在列中我们有引导T^统计数据)我们找到了Tcrit.我们准确地观察α=0.05更显着的百分比T^统计数据(更显着意味着绝对值大于Tcrit.)。

我们报告i-th 模型组件显着,如果它Ti^>Tcrit.

最后一步可以用这段代码完成

p.value<-0.05 # The target alpha threshold
digits<-1000000
library(gtools) # for binsearch

pValueFun<-function(cj)
{
   mean(apply(abs(TestStatBoot)>cj/digits,1,sum)>=1,na.rm=T)
}
ans<-binsearch(pValueFun,c(0.5*digits,100*digits),target=p.value)
p.level<-(1-pnorm(q=ans$where[[1]]/digits))*2 #two-sided.