Corrigindo valores de p para vários testes em que os testes estão correlacionados (genética)

24

Eu tenho valores de p em muitos testes e gostaria de saber se há realmente algo significativo após a correção para vários testes. A complicação: meus testes não são independentes. O método em que estou pensando (uma variante do Método do Produto de Fisher, Zaykin et al., Genet Epidemiol , 2002) precisa da correlação entre os valores de p.

Para estimar essa correlação, atualmente estou pensando em casos de inicialização, executando as análises e correlacionando os vetores resultantes dos valores de p. Alguém tem uma idéia melhor? Ou ainda uma idéia melhor para o meu problema original (corrigindo vários testes em testes correlatos)?

Antecedentes: Regresso logisticamente se meus pacientes sofrem ou não de uma doença específica na interação entre o genótipo (AA, Aa ou aa) e uma covariável. No entanto, o genótipo é na verdade um monte (30-250) de polimorfismos de nucleotídeo único (SNPs), que certamente não são independentes, mas no desequilíbrio de ligação.

S. Kolassa - Restabelecer Monica
fonte

Respostas:

29

Este é realmente um tópico importante nos estudos de análise Genomewide (GWAS)! Não sei se o método que você está pensando é o mais apropriado nesse contexto. O agrupamento dos valores de p foi descrito por alguns autores, mas em um contexto diferente (estudos de replicação ou metanálise, ver, por exemplo, (1) uma revisão recente). A combinação de valores p de SNP pelo método de Fisher é geralmente usada quando se deseja derivar um valor p único para um determinado gene; isso permite trabalhar no nível do gene e reduzir a quantidade de dimensionalidade dos testes subsequentes, mas, como você disse, a não independência entre os marcadores (decorrente da colocação espacial ou desiquilíbrio de ligação, LD) apresenta um viés. Alternativas mais poderosas dependem de procedimentos de reamostragem,

Minhas principais preocupações com o bootstraping (com substituição) são que você está introduzindo uma forma artificial de parentesco ou, em outras palavras, cria gêmeos virtuais, alterando o equilíbrio de Hardy-Weinberg (mas também a frequência alélica mínima e a taxa de chamadas). Este não seria o caso de uma abordagem de permutação em que você permite rótulos individuais e mantém os dados de genotipagem como estão. Normalmente, o software plink pode fornecer valores p brutos e permutados, embora use (por padrão) uma estratégia de teste adaptável com uma janela deslizante que permite interromper a execução de todas as permutações (por exemplo, 1000 por SNP) se parecer que o SNP está sob consideração não é "interessante"; ele também tem opção para calcular o maxT, consulte a ajuda online .

Mas, dado o baixo número de SNPs que você está considerando, eu sugeriria confiar nos testes baseados em FDR ou maxT, conforme implementados no pacote R multtest (veja mt.maxT), mas o guia definitivo para reamostrar estratégias de aplicação genômica é Multiple Procedures Testing with Applications to Genomics , de Dudoit & van der Laan (Springer, 2008). Veja também o livro de Andrea Foulkes sobre genética com R , que é revisado no JSS. Ela tem ótimo material em vários procedimentos de teste.

Notas adicionais

Muitos autores apontaram o fato de que métodos simples de correção de múltiplos testes, como o Bonferroni ou o Sidak, são muito rigorosos para ajustar os resultados para os SNPs individuais. Além disso, nenhum desses métodos leva em consideração a correlação existente entre os SNPs devido ao LD, que marca a variação genética nas regiões gênicas. Outras alternativas foram propostas, como um derivado do método de Holm para comparação múltipla (3), modelo oculto de Markov (4), FDR condicional ou positivo (5) ou seu derivado (6), para citar alguns. As chamadas estatísticas de gap ou janela deslizante provaram ser bem-sucedidas em alguns casos, mas você encontrará uma boa revisão em (7) e (8).

Também ouvi falar de métodos que fazem uso efetivo da estrutura do haplótipo ou LD, por exemplo (9), mas nunca os usei. Eles parecem, no entanto, mais relacionados à estimativa da correlação entre marcadores, e não ao valor de p, como você quis dizer. Mas, na verdade, é melhor pensar em termos da estrutura de dependência entre estatísticas sucessivas de teste, do que entre valores de p correlacionados.

Referências

  1. Cantor, RM, Lange, K. e Sinsheimer, JS. Priorizando os resultados do GWAS: uma revisão dos métodos estatísticos e recomendações para sua aplicação . Sou J Hum Genet. 2010 86 (1): 6–22.
  2. Corley, RP, Zeiger, JS, Crowley, T et al. Associação de genes candidatos à dependência de drogas antissociais em adolescentes . Dependência de Drogas e Álcool 2008 96: 90–98.
  3. Dalmasso, C, Génin, E e Trégouet DA. Um Procedimento de Holm Ponderado Contabilizando Frequências Alélicas em Estudos de Associação Genomewide . Genetics 2008 180 (1): 697–702.
  4. Wei, Z, Sun, W, Wang, K e Hakonarson, H. Testes Múltiplos em Estudos de Associação em Todo o Genoma através de Modelos Markov Ocultos . Bioinformatics 2009 25 (21): 2802-2808.
  5. Broberg, P. Uma revisão comparativa das estimativas da proporção de genes inalterados e a taxa de falsas descobertas . BMC Bioinformatics 2005 6: 199.
  6. Necessidade, CA, Ge, D, Arma, ME, et a. Uma investigação em todo o genoma de SNPs e CNVs na esquizofrenia . PLoS Genet. 2009 5 (2): e1000373.
  7. Han, B, Kang, HM e Eskin, E. Correção rápida e precisa de testes múltiplos e estimativa de energia para milhões de marcadores correlacionados . PLoS Genetics 2009
  8. Liang, Y e Kelemen, A. Avanços e desafios estatísticos para analisar dados snp de alta dimensão correlacionados em estudo genômico para doenças complexas . Pesquisas Estatísticas 2008 2: 43–60. - a melhor revisão recente de sempre
  9. Nyholt, DR. Uma correção simples para testes múltiplos de polimorfismos de um único nucleotídeo no desequilíbrio de ligação entre si . Sou J Hum Genet. 2004 74 (4): 765-769.
  10. Nicodemos, KK, Liu, W, Chase, GA, Tsai, AA e Fallin, MD. Comparação do erro do tipo I para correções de teste múltiplas em grandes estudos de polimorfismo de nucleotídeo único usando componentes principais versus algoritmos de bloqueio de haplótipos . BMC Genetics 2005; 6 (suplemento 1): S78.
  11. Peng, Q, Zhao, J e Xue, F. Testes de intervalo de confiança de autoinicialização baseados em PCA para associação gene-doença envolvendo múltiplos SNPs . BMC Genetics 2010, 11: 6
  12. Li, M, Romero, R., Fu, WJ e Cui, Y (2010). Mapeando interações haplótipo-haplótipo com LASSO adaptável . BMC Genetics 2010, 11:79 - embora não esteja diretamente relacionado à questão, abrange análise baseada em haplótipos / efeito epistático
chl
fonte
11
Uau, obrigado por enfrentar todo esse problema! Entendo suas dúvidas sobre o bootstrap e estou quase convencido. Penso que minha principal complicação é a covariável numérica que tenho que certamente será necessária (por si só ou em interação com o genótipo), e isso parece descartar o mt.maxT e o plink, embora eu possa precisar procurar novamente o plink. Mas certamente vou procurar nas referências que você forneceu!
S. Kolassa - Restabelece Monica
Você sempre pode trabalhar com os resíduos do seu GLM para obter informações sobre suas covariáveis, embora tenha perdido algum Df que será difícil de contabilizar ou reintroduzir posteriormente (por exemplo, para calcular o valor p).
chl
Hm, resíduos da minha regressão logística? Isso seria legítimo?
S. Kolassa - Restabelece Monica
Sim, porque não? Não é incomum remover a variação contabilizada por outras covariáveis ​​e depois passar para a análise de segundo nível com seus dados residualizados. Geralmente é mais rápido (por exemplo, o plink é bem lento com as covariáveis ​​categóricas, enquanto as com contínuas são aceitáveis; snpMatrixou simplesmente glm()funciona muito bem nesse ponto, mas você não pode incorporar muitos SNPs dentro de glm()...); o problema é que obter o valor p corrigido no final de sua 2ª análise é bastante complicado (porque você precisa contabilizar os parâmetros já estimados).
chl
Para uma ilustração de como as pessoas estão trabalhando com resíduos, veja, por exemplo, p. 466 de Heck et al. A investigação de 17 genes candidatos para traços de personalidade confirma os efeitos do gene HTR2A na busca de novidades. Genes, cérebro e comportamento (2009) vol. 8 (4) pp. 464-72
chl
2

Usar um método como o bonferroni é bom, o problema é que, se você tiver muitos testes, provavelmente não encontrará muitas "descobertas".

Você pode usar a abordagem FDR para testes dependentes (consulte aqui para obter detalhes ). O problema é que não tenho certeza se você pode dizer com antecedência se suas correlações são todas positivas.

Em R, você pode executar um FDR simples com p.adjust. Para coisas mais complexas, eu daria uma olhada multcomp , mas não procurei por soluções em casos de dependências.

Boa sorte.

Tal Galili
fonte
11
Oi Tal, obrigado! Bonferroni não parece apropriado para mim - se um dos meus SNPs é causal e outros estão associados a ele, deve haver um sinal, e Bonferroni sempre pareceu muito conservador para mim (eu geralmente prefiro a correção gradual de Holm). O FDR ao qual você vincula e o p.adjust não considera evidências combinadas (e o FDR ainda exige que eu compreenda a correlação dos meus testes, a pergunta original). o multcomp pode ajudar, embora à primeira vista pareça que lide mais com vários testes em um único modelo, enquanto eu tenho vários modelos. Eu vou cavar mais fundo ...
S. Kolassa - Restabelece Monica
Olá Stephan. Eu entendo, desculpe por não ajudar mais. Boa sorte! Tal
Tal Galili
Olá Stephan, eu ainda acho que você ainda pode usar o método = BY (para o procedimento Benjamini Hochberg Yekuteli) em p.adjust em R, como apontado por Tal. Definitivamente, usar Bonferroni pode ser conservador.
suncoolsu
suncoolsu, acho que esse método só funciona quando a correlação é positiva (não negativa) entre as variáveis. Felicidades.
Tal Galili
2

Acho que os modelos normais multivariados estão sendo usados ​​para modelar os valores p correlacionados e obter o tipo certo de várias correções de teste. Correção rápida e precisa de testes múltiplos e estimativa de energia para milhões de marcadores correlacionados. O PLoS Genet 2009 fala sobre eles e também fornece outras referências. Parece semelhante ao que você estava falando, mas acho que, além de obter uma correção global mais precisa do valor de p, o conhecimento da estrutura do LD também deve ser usado para remover positivos falsos decorrentes de marcadores correlacionados com marcadores causais.

alta largura de banda
fonte
2

Estou procurando uma solução funcional para exatamente o mesmo problema. O melhor que encontrei é o Bootstrap nulo e irrestrito, introduzido por Foulkes Andrea em seu livro Applied Statistical Genetics with R (2009) . Ao contrário de todos os outros artigos e livros, ele considera especificamente as regressões. Além de outros métodos, ele aconselha o Null Unrestricted Bootstrap, que é adequado onde não é possível calcular facilmente resíduos (como no meu caso, onde modelo muitas regressões independentes (basicamente correlações simples), cada uma com a mesma variável de resposta e snip diferente). Eu achei esse método também chamado de método 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
+    }
+ }

TestStatBootT^Tcrit.α=0,05T^Tcrit.

EuTEu^>Tcrit.

A última etapa pode ser realizada com este código

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.
Adam Ryczkowski
fonte