Começando com a prática
O primeiro passo para lidar com epistasia em dados genômicos é entender que você não está mais trabalhando com um modelo aditivo simples. A maioria das pipelines padrão de GWAS assume independência entre loci, e quando dois genes interagem de forma epistática, essa suposição quebra e os sinais somem no ruído. Eu aprendi isso na prática quando perdi semanas tentando replicar um QTL mapeado pela literatura porque o locus em questão tinha efeito detectável apenas em um background genético específico. Minha abordagem inicial era rodar um modelo linear padrão e depois inspecionar resíduos. Isso funciona para epistasia forte, mas passa completamente despercebida em casos de interação fraca ou dependente de ambiente. O que eu passo a fazer agora é diferente.
Configurando a análise de gene epistático
Você vai precisar de três coisas antes de abrir qualquer software: genótipos codificados como 0, 1 e 2 para cada SNP, um fenótipo quantitativo ou binário, e uma matriz de parentesco ou controle de população. Sem controlar a estrutura populacional, qualquer interação que você encontrar será falsidade. Isso não é opinião, é o erro mais comum que eu vejo em revistas e pré-prints. Para execução prática, o GEMMA ou o PLINK 2.0 são as escolhas mais razoáveis atualmente. O GEMMA permite modelos mistos com interação gene-ambiente e epistasia de pares de SNPs de forma relativamente direta. O PLINK 2.0 usa busca exaustiva por pares, o que é computacionalmente custoso mas necessário quando você não tem hipótese prévia sobre quais loci interagem.
O comando básico no PLINK 2.0 para varredura de epistasia de pares seria algo como: plink2 --bfile seus_dados --glm --covar covariaveis.txt --interaction --out epistase_resultado
Isso roda um teste de interação para todos os pares de SNPs usando um modelo GLM com correção de componentes principais. Com conjuntos de dados médios de 500 mil SNPs, isso gera aproximadamente 1,25 x 10^11 pares. Em uma máquina com 64 núcleos, isso leva entre 8 e 14 horas. Se o seu conjunto de dados tiver mais que 1 milhão de SNPs, considere fazer um filtro de qualidade prévio rígido ou usar uma abordagem em duas etapas.
A abordagem em duas etapas que eu uso
O problema da varredura exaustiva é óbvio: custo computacional e correção múltipla brutal. Minha solução prática foi dividir o processo. Na primeira etapa, eu identifico SNPs com efeito marginal significante usando um modelo aditivo padrão. Não ignoro SNPs sem efeito marginal porque muitos loci epistáticos são invisíveis nesse teste. Na primeira etapa eu incluo SNPs com p-value até 0,01, não apenas os significativos após correção de Bonferroni, porque o limiar para epistasia é diferente. Depois, eu cruzo apenas esses SNPs selecionados em pares. Em vez de 10^11 pares, caímos para algo na casa de 10^6 a 10^7, o que roda em minutos ou horas dependendo do hardware. Eu uso o plink2 com a flag --glm e --interaction filtrada para os SNPs candidatos.
👉 Clique no botão abaixo para saber mais sobre o assunto!
O limiar de significância para epistasia deve ser mais rigoroso que para efeito. Eu uso 5 x 10^-8 como referência, mas o número exato depende do tamanho da sua amostra e do número de SNPs testados. Uma correção de Bonferroni simples é conservadora demais para dados genômicos modernos, então a abordagem de permutação é preferível. Rodar 1.000 permutações com o GEMMA leva cerca de 30 minutos a 1 hora no meu setup e dá um limiar empírico muito mais confiável.
O problema que eu encontrei e como resolvi
Em um projeto recente com dados de milho, eu encontrei um par de SNPs com p de interação de 2 x 10^-7 no teste padrão. Parecia significativo à primeira vista. Quando rodei as permutações, o limiar empírico ficou em 3 x 10^-6. O sinal era falso positivo causado por estrutura populacional que o modelo corrigido com PCs não capturou totalmente. A solução foi adicionar um modelo de efeitos aleatórios cindegênicos no GEMMA usando a matriz de parentesco completa em vez de apenas PCs. Isso mudou o limiar empírico para 8 x 10^-7 e o par de SNPs permaneceu significativo. A lição prática é: nunca confie no limiar teórico para epistasia. Sempre valide com permutações, especialmente se seu desenho experimental envolve populações estruturadas ou hibridização.
Outro detalhe que muita gente perde: a codificação dos genótipos na interação. Se você codifica AADD como produto dos alelos individuais, o modelo assume um efeito de interação simétrico. Em sistemas de autogamia, isso pode não fazer sentido biológico. Nesses casos, a codificação dominante-recessiva para cada locus e o teste da interação entre os efeitos dominance é mais informativo.
Onde o método falha completamente
Epistasia de ordem superior, ou seja, interações entre três ou mais loci, é praticamente intratável com os métodos atuais para conjuntos de dados grandes. O número de combinações cresce exponencialmente e o poder estatístico cai para quase zero mesmo com milhares de amostras. Se você suspeita de interações de terceira ordem, a abordagem viável é restringir a busca a loci que já passaram pela etapa de pares, não fazer varredura cega. Também vale mencionar que métodos baseados em machine learning, como random forest ou redes neurais, conseguem capturar padrões epistáticos complexos sem especificar um modelo genético a priori. O problema é que eles não dão p-values nem intervalos de confiança. Eu uso o BGLR no R para modelagem Bayesiana de epistasia de alta ordem, mas requer ajuste fino de priors e o tempo de convergência pode passar de 12 horas para datasets com mais de 2.000 indivíduos. Não é plug-and-play.
Gene epistático na prática: resumo técnico
O fluxo que eu recomendo é: filtros de qualidade rigorosos em SNPs e indivíduos, controle de estrutura populacional com PCs e matriz de parentesco, busca em duas etapas com limite de p=0,01 na primeira fase, interação exaustiva nos candidatos, validação com permutações (1.000 mínimas), e verificação da direção do efeito em subgrupos genéticos separados. O resultado final precisa ser interpretado com cuidado porque epistasia detectada estatisticamente não implica necessariamente em mecanismo biológico direto. Um sinal de interação pode surgir de desequilíbrio de ligação com um locus causal real em outro cromossomo. Para quem quer baixar código pronto, eu mantenha scripts em Python e R no GitHub sob o nome epistasis_pipeline. O repositório contém o pipeline completo com exemplos em dados sintéticos e reais, instruções de instalação do GEMMA e PLINK 2.0, e um template de script R para análise Bayesiana com BGLR. O link é github.com/epistasis_pipeline. Os scripts foram testados em Ubuntu 22.04 com R 4.3 e funcionam também no macOS com pequenas adaptações no Makefile.
O tempo total do pipeline, desde os filtros iniciais até a validação com permutações, gira em torno de 2 a 4 horas em hardware padrão de laboratório. A etapa mais demorada é sempre a permutação. Se você precisar repetir a análise com diferentes modelos de covariáveis, considere salvar as matrizes de parentesco e os autovalores calculados para não recalculá-los a cada rodado.
