Análise de Expressão Diferencial

A análise de expressão diferencial responde à pergunta mais comum na transcriptômica: quais genes mudam sua expressão entre a condição A e a condição B? Parece simples, mas os desafios estatísticos — dados de contagem, alta dimensionalidade, tamanhos de amostra pequenos, testes múltiplos — requerem métodos especializados desenvolvidos especificamente para dados de RNA-seq.

O Problema Estatístico

Dado:

• n₁ amostras na condição A (por exemplo, tumor)
• n₂ amostras na condição B (por exemplo, normal)
• Medições de expressão para G ~20.000 genes em cada amostra

Encontrar: quais genes diferem significativamente entre as condições?

Por que você não pode usar um teste t em contagens brutas:

1. As contagens de RNA-seq são inteiros, não contínuos — um gene com 0 contagens em uma amostra, mas 100 em outra, não pode ser normalmente distribuído

2. A variância não é constante nos níveis de contagem. Para contagens de Poisson, variância = média. Para RNA-seq, variância >> média (superdispersão) — esta é a variabilidade biológica além da variação de amostragem

3. A profundidade de sequenciamento varia entre amostras — uma amostra com 2× mais leituras totais mostrará 2× mais contagens por gene, que é um artefato técnico, não biologia

Normalização Primeiro

Antes do teste, as contagens devem ser normalizadas para remover os efeitos de profundidade de sequenciamento:

Normalização por fator de tamanho do DESeq2: para cada amostra, calcule um fator de tamanho como a mediana de (contagem do gene g na amostra j) / (média geométrica do gene g em todas as amostras). Dividir cada contagem pelo fator de tamanho da amostra torna as amostras comparáveis.

Isso é mais robusto do que a escala simples por milhão porque é resistente à influência de genes altamente expressos (um gene expresso em 10% das leituras totais dominaria a escala simples).

Quando usar qual normalização:

• Fatores de tamanho DESeq2 (normalização intra-experimento): corrigem as diferenças no tamanho da biblioteca entre amostras no mesmo experimento
• TMM (edgeR): semelhante ao DESeq2, assume que a maioria dos genes não é diferencialmente expressa
• TPM (para comparação entre experimentos): normaliza tanto para o tamanho da biblioteca QUANTO para o comprimento do gene; apropriado para comparar níveis de expressão em amostras de diferentes estudos
• CPM/RPKM/FPKM: métodos mais antigos; evite RPKM/FPKM geralmente devido a inconsistências matemáticas na comparação entre amostras

O Modelo Binomial Negativo

O modelo apropriado para dados de contagem de RNA-seq é a distribuição binomial negativa (BN):

P(Y = k) = Γ(k + μ/α) / (Γ(μ/α) × k!) × (αμ)^k / (1 + αμ)^(k + μ/α)

Onde:

• μ = contagem média (nível médio de expressão)
• α = parâmetro de dispersão (superdispersão; α=0 reduz para Poisson)

A variância é: Var(Y) = μ + α × μ²

A dispersão α captura a variabilidade extra além do ruído de amostragem. Varia por gene e é tipicamente estimada a partir dos dados.

O Problema de n Pequeno e o Encolhimento

Com apenas n=3–5 replicatas por grupo, a estimativa de variância por gene é não confiável. Um gene com 3 amostras de controle pode mostrar alta variância apenas por acaso — levando a intervalos de confiança inflados e perda de poder.

Encolhimento empírico Bayesiano: tanto o DESeq2 quanto o edgeR "pegam emprestado poder" de outros genes. A ideia: a dispersão de um gene com contagem 50 deve ser semelhante às dispersões de outros genes com contagens semelhantes. Ao reunir informações de milhares de genes, você obtém estimativas de dispersão muito melhores para qualquer gene individual.

O DESeq2 ajusta uma curva: dispersão como função da expressão média (a "tendência de dispersão"). A estimativa de dispersão de cada gene é então encolhida em direção a essa tendência. Genes com poucas replicatas (onde a estimativa não é confiável) são encolhidos mais; genes onde os dados apoiam fortemente uma alta dispersão são encolhidos menos.

Esse encolhimento é uma forma de regularização — exatamente o mesmo princípio que a regressão Ridge. Reduz a variância na estimativa ao custo de introduzir algum viés, mas o efeito líquido no desempenho do teste é fortemente positivo.

Fluxo de Trabalho do DESeq2

Matriz de contagens brutas (genes × amostras)
    ↓
Criar objeto DESeqDataSet
    ↓
Estimar fatores de tamanho (normalização)
    ↓
Estimar dispersões (por gene, ajustar tendência, encolher)
    ↓
Ajustar GLM binomial negativo por gene
    ↓
Teste de Wald ou LRT para log fold changes
    ↓
Encolher log fold changes (apeglm/ashr)
    ↓
Aplicar correção de testes múltiplos (Benjamini-Hochberg)
    ↓
Tabela de resultados: log2FC, SE, estatística Wald, p-valor, padj

Fórmula de Design

A fórmula de design especifica quais comparações fazer:

# Simple two-group comparison
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = metadata,
                              design = ~condition)

# Controlling for batch effect
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = metadata,
                              design = ~batch + condition)

# Interaction: does treatment effect differ between genotypes?
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = metadata,
                              design = ~genotype + treatment + genotype:treatment)

A fórmula de design é um modelo linear no framework de modelo linear generalizado. Cada gene é testado com a mesma estrutura de modelo; os coeficientes e testes correspondem às colunas da matriz de design.

Encolhimento do Log Fold Change

Mesmo após o encolhimento da dispersão, as estimativas de log fold change podem ser ruidosas para genes de baixa contagem. O lfcShrink() do DESeq2 aplica um prior de encolhimento adicional às estimativas de fold change:

• O método apeglm (prior t de Student adaptativo) é a recomendação padrão
• Encolhe grandes LFCs para genes de baixa contagem em direção a zero enquanto preserva estimativas para genes de alta contagem
• NÃO muda o p-valor ou o status de significância — apenas a estimativa de LFC para visualização e classificação

LFCs encolhidos são apropriados para gráficos MA, gráficos vulcão e classificação de genes para análise de vias. LFCs não encolhidos podem produzir valores extremos para genes de baixa contagem que dominam as classificações.

Interpretando Resultados

A Tabela de Resultados

Limiar de significância padrão: padj < 0,05 (5% de FDR) E |log2FoldChange| > 1 (pelo menos mudança de 2 vezes). O segundo critério é um filtro de significância biológica — mudanças de fold change estatisticamente significativas, mas minúsculas, podem não ser biologicamente relevantes.

O Gráfico MA

O gráfico MA mostra o log fold change (eixo y) vs. a expressão média (eixo x). Com encolhimento do fold change:

• Genes de alta contagem: estimativas de LFC se espalham em torno dos valores verdadeiros; genes significativos claramente visíveis
• Genes de baixa contagem: LFCs encolhidos em direção a zero; estimativas ruidosas não dominam o gráfico

Pontos vermelhos = significativos (padj < 0,05). Um gráfico MA bem comportado mostra simetria em torno de LFC = 0 para genes não significativos, com genes significativos se estendendo para longe do centro.

O Gráfico Vulcão

O gráfico vulcão mostra -log₁₀(p-valor) vs. log₂ fold change:

• Eixo x: magnitude da mudança (mais extremo = maior efeito)
• Eixo y: significância estatística (mais alto = mais significativo)
• Superior-esquerdo: significativamente regulados negativamente
• Superior-direito: significativamente regulados positivamente
• Meio inferior: nem significativo nem grande efeito

Pontos coloridos/rotulados são tipicamente aqueles que passam nos limiares de significância e fold change.

Filtragem Independente

O DESeq2 executa automaticamente a filtragem independente: remove genes com contagens médias muito baixas da correção de testes múltiplos. Por quê? Esses genes essencialmente não têm poder para ser detectados como significativos, então consomem o orçamento de correção de testes múltiplos sem contribuir com descobertas reais.

O filtro é "independente" porque usa um critério (contagem média) independente da estatística de teste sob a nula — portanto, não envolve viés nos resultados, apenas aumenta o poder reduzindo o número de hipóteses testadas.

Teste LRT vs. Teste de Wald

Teste de Wald (padrão): testa se um coeficiente específico (por exemplo, log2FC para condição) é significativamente diferente de zero. Rápido; testa diretamente um coeficiente.

Teste de razão de verossimilhança (LRT): compara o modelo completo com um modelo reduzido (sem o termo de interesse). Mais apropriado para:

• Testar se a condição explica variância significativa no geral (em um conjunto de condições)
• Comparar modelos com complexidade diferente
• Testar efeitos significativos em qualquer direção simultaneamente

Use LRT quando você tem múltiplas condições e quer identificar genes que mudam em qualquer delas: DESeq2::DESeq(dds, test="LRT", reduced=~1).

edgeR: A Alternativa

O edgeR usa um framework binomial negativo semelhante, mas implementa detalhes diferentes:

• Normalização TMM em vez dos fatores de tamanho do DESeq2
• Abordagem de quasi-verossimilhança (teste F QL) para melhor calibração com amostras pequenas — agora o modo recomendado
• Estimativa de dispersão por tag (por gene) com pooling empírico Bayesiano

Para a maioria dos propósitos, DESeq2 e edgeR fornecem resultados muito semelhantes. DESeq2 tem mais recursos integrados (séries temporais, designs complexos, encolhimento). A abordagem QL do edgeR às vezes é preferida para tamanhos de amostra muito pequenos.

limma-voom: O Padrão para Designs Complexos

Para designs multi-fatoriais complexos, limma-voom é frequentemente a ferramenta mais flexível:

1. Voom: transforma valores log-CPM com pesos de precisão que levam em conta a relação média-variância nos dados de contagem
2. limma: aplica modelagem linear com encolhimento de variância empírico Bayesiano

O ponto forte do limma: ele lida com qualquer design de modelo linear, incluindo covariáveis contínuas, termos de interação, efeitos de lote e comparações de séries temporais. O framework completo de modelo linear do limma é mais flexível do que DESeq2/edgeR para designs experimentais complexos.

Análise de Séries Temporais

Para experimentos com múltiplos pontos de tempo, o LRT do DESeq2 identifica genes com qualquer mudança ao longo do tempo, enquanto maSigPro e ImpulseDE2 podem identificar genes com padrões temporais específicos (aumento monotônico, pico transitório, etc.).

Para padrões periódicos (ritmos circadianos), RAIN e JTK_CYCLE usam testes não paramétricos especificamente projetados para detecção de ritmo.

Considerações Práticas

Tamanho mínimo de amostra: n ≥ 3 por grupo é um mínimo absoluto; n ≥ 5–6 é fortemente recomendado. Com n=2, a única opção estatística é tratar as amostras como replicatas técnicas — o que raramente se justifica.

Verifique amostras outlier antes da análise: PCA ou agrupamento hierárquico de amostras deve mostrar agrupamento por condição, não por amostra. Uma amostra que agrupa com o grupo errado precisa de investigação (possível erro de rotulagem, contaminação ou baixa qualidade).

Filtragem de baixa contagem manual: a filtragem independente do DESeq2 é automática, mas filtrar manualmente antes da análise (por exemplo, manter apenas genes com CPM > 1 em pelo menos 3 amostras) garante que você analise genes com evidências de expressão em pelo menos algumas amostras.

A interpretação biológica requer contexto: 500 genes diferencialmente expressos é o início, não a resposta. Próximos passos: enriquecimento de vias, análise de rede, validação dos principais candidatos, integração com outros tipos de dados (ATAC-seq para elementos regulatórios, expressão de proteínas para confirmação funcional).