Objetivos de aprendizagem

1. Compreender os passos necessários para gerar resultados de comparação (teste de Wald)
2. Resumir os diferentes níveis de filtragem de genes
3. Compreender a retração do log2 fold change

Exploração de resultados

Por padrão, o DESeq2 utiliza o teste de Wald para identificar genes expressos de forma diferencial entre duas amostras. Dado os fatores utilizados na fórmula de design e a quantidade de níveis existentes, podemos extrair resultados para diversas comparações distintas. A seguir, apresentamos como obter resultados do objeto dds e fornecemos explicações sobre sua interpretação.

Nota: O teste de Wald também pode ser utilizado para variáveis contínuas. Se a variável de interesse na fórmula de design for um valor contínuo, o log2FoldChange reportado representa a mudança por unidade dessa variável.

1. Especificando a comparação

---

No nosso conjunto de dados, temos três categorias de amostras, permitindo três comparações possíveis aos pares:

1. Controle vs. Superexpressão de Mov10
2. Controle vs. Knockout de Mov10
3. Superexpressão de Mov10 vs. Knockout de Mov10

Estamos interessados apenas nas comparações 1 e 2. Quando criamos nosso objeto dds, fornecemos ~sampletype como fórmula de design, indicando que sampletype é o fator principal de interesse.

Para indicar quais duas amostras queremos comparar, precisamos especificar a contraste. Utilize a função results() do DESeq2 para extrair os resultados desejados.

A contraste pode ser especificada de duas formas diferentes (a primeira é mais comum):

1. A contraste pode ser fornecida como um vetor de caracteres com três elementos: o nome do fator (de interesse) na fórmula de design, o nome do nível do fator a ser comparado, e o nível base do fator. O nível apresentado por último serve como referência. A sintaxe é:

contrast <- c("fator", "nivel_comparar", "nivel_base")
results(dds, contrast = contrast)

2. A contraste pode ser fornecida como uma lista de dois vetores de caracteres: o nome da dobra que varia com o nível de interesse, e o nome da dobra que varia com o nível base. Estes nomes devem corresponder exatamente aos elementos de resultsNames(object).

resultsNames(dds) 
contrast <- list(resultsNames(dds)[1], resultsNames(dds)[2])
results(dds, contrast = contrast)

Alternativamente, se você tiver apenas dois níveis de fator, não precisará especificar a contraste (ou seja, results(dds)). Neste caso, o DESeq2 selecionará seu nível base fator de acordo com a ordem alfabética.

Primeiro, avaliaremos as mudanças de expressão entre as amostras de superexpressão MOV10 e o controle. Portanto, utilizaremos o primeiro método para especificar a contraste e criar um vetor de caracteres:

contrast_oe <- c("sampletype", "MOV10_overexpression", "control")

2. Resultados

---

Agora que criamos a contraste, podemos utilizá-la como entrada para a função results().

res_oe <- results(dds, contrast=contrast_oe, alpha = 0.05)

Nota: Para nossa análise, além do parâmetro de contraste, fornecemos o valor 0.05 para o parâmetro alpha. Discutiremos isso com mais detalhes quando falarmos sobre filtragem a nível de gene.

O resultado retornado é um objeto DESeqResults, que é uma subclasse simples de DataFrame. Em muitos aspectos, pode ser tratado como um data frame (ou seja, ao acessar/subconjugar dados), mas é importante reconhecer que existem diferenças nas etapas posteriores (como visualização).

Agora, vamos verificar quais informações estão armazenadas nos resultados:

res_oe %>% 
data.frame() %>% 
View()

Podemos utilizar a função mcols() para extrair informações sobre o que cada coluna representa:

mcols(res_oe, use.names=T)

• baseMean: Média das contagens normalizadas de todas as amostras
• log2FoldChange: Mudança em dobramento em log2
• lfcSE: Erro padrão
• stat: Estatística de Wald
• pvalue: Valor-p do teste de Wald
• padj: Valores-p ajustados por Benjamini-Hochberg

3. Valores-p

---

O valor-p é a probabilidade utilizada para determinar se há evidência para rejeitar a hipótese nula. Valores-p menores significam que há evidências mais fortes para支持 a hipótese alternativa. No entanto, como estamos testando cada gene individualmente, precisamos corrigir esses valores-p para múltiplos testes.

A coluna padj nos resultados representa os valores-p ajustados para múltiplos testes, sendo a coluna mais importante nos resultados. Tipicamente, limiares como padj < 0.05 são bons pontos de partida para identificar genes importantes. O método padrão de correção para múltiplos testes no DESeq2 é a implementação de Benjamini-Hochberg para taxa de falsos descobrimentos (FDR). Outros métodos de correção estão disponíveis e podem ser alterados adicionando o parâmetro padjustMethod à função results().

4. Filtragem

---

Observe atentamente nossos resultados. Ao navegar, você notará que para certos genes, existem valores NA nas colunas pvalue e padj. O que isso significa?

Valores ausentes indicam genes que foram filtrados como parte da função DESeq(). Antes de realizar a análise de expressão diferencial, é benéfico ignorar genes que têm pouca ou nenhuma chance de serem detectados como diferencialmente expressos. Isso aumenta o poder de detecção de genes diferencialmente expressos. O DESeq2 não remove nenhum gene da matriz de contagem bruta; portanto, todos os genes aparecerão na sua tabela de resultados. Os genes omitidos pelo DESeq2 satisfazem um dos três critérios de filtragem:

1. Genes com contagem zero em todas as amostras

Se em uma linha, todas as amostras tiverem contagem zero, não há informação de expressão, portanto esses genes não são testados.

res_oe[which(res_oe$baseMean == 0),] %>% 
data.frame() %>% 
View()

Nestes genes, a coluna baseMean será zero, e as estimativas de log2 fold change, valores-p e valores-p ajustados serão definidos como NA.

2. Genes com valores discrepantes extremos de contagem

A função DESeq() calcula um teste de diagnóstico de valores discrepantes para cada gene e cada amostra, chamado distância de Cook. A distância de Cook mede a influência de uma amostra individual no coeficiente de ajuste do gene, e valores grandes indicam contagens discrepantes. Genes com distância de Cook acima de um limiar são marcados, mas o标记 requer pelo menos 3 repetições, pois é difícil determinar qual amostra pode ser um valor discrepante com apenas 2 repetições. Podemos desativar essa filtragem usando o parâmetro cooksCutoff na função results().

res_oe[which(is.na(res_oe$pvalue) & 
             is.na(res_oe$padj) &
             res_oe$baseMean > 0),] %>% 
data.frame() %>% 
View()

3. Genes com média baixa de contagens normalizadas

O DESeq2 define um limiar de contagem baixa empiricamente determinado a partir dos seus dados, onde a proporção de genes significativos pode ser aumentada reduzindo o número de genes considerados para múltiplos testes. Isso se baseia na ideia de que genes com contagens muito baixas geralmente não apresentam diferenças significativas devido à alta dispersão.

No valor especificado pelo usuário (alpha = 0.05), o DESeq2 avalia a mudança no número de genes significativos à medida que filta porçõesprogressivamente maiores de genes com base em sua contagem média, como mostrado na figura acima. O ponto onde o número de genes significativos atinge o pico é o limiar de média baixa usado para filtrar genes testados múltiplas vezes. Há também um parâmetro para desativar a filtragem definindo independentFiltering = F.

res_oe[which(!is.na(res_oe$pvalue) & 
             is.na(res_oe$padj) & 
             res_oe$baseMean > 0),] %>% 
data.frame() %>% 
View()

Nota: Por padrão, o DESeq2 executa a filtragem descrita acima; no entanto, outras ferramentas de DE (como EdgeR) não o fazem. A filtragem é uma etapa necessária mesmo se você estiver usando o método de quasi-verossimilhança do limma-voom e/ou edgeR. Ao usar outras ferramentas, siga as etapas de pré-filtragem conforme descrito no guia do usuário no Bioconductor, pois geralmente apresentam melhor desempenho.

5. Mudança de dobramento

---

Outra coluna importante nos resultados é o log2FoldChange. Para listas grandes de genes, é difícil extrair relevância biológica significativa. Para ajudar a aumentar o rigor, um limiar de fold change pode ser adicionado. Lembre-se de que ao definir esse valor, estamos lidando com mudanças de dobramento em log2, então um corte de log2FoldChange < 1 se traduz em dobramento real de 2.

A mudança de dobramento nos resultados é calculada como:

log2 (contagens_normalizadas_grupo1 / contagens_normalizadas_grupo2)

O problema é que essas estimativas de mudança de dobramento não são totalmente precisas, pois não consideram a dispersão que observamos em contagens baixas de leituras. Para resolver isso, é necessário ajustar a mudança de dobramento em log2.

Retração de LFC

• Estimativas mais precisas de LFC

Para gerar estimativas mais precisas de log2 fold change (LFC), o DESeq2 permite que as estimativas de LFC sejam retraídas para zero quando a informação do gene é escassa, o que pode incluir:

• Contagens baixas
• Valores de dispersão altos

A retração de LFC usa informações de todos os genes para gerar estimativas mais precisas. Especificamente, a distribuição das estimativas de LFC de todos os genes (como priori) é usada para retrair as estimativas de LFC de genes com pouca informação ou alta dispersão para estimativas de LFC mais prováveis (mais baixas).

Na figura acima, temos um exemplo usando genes verdes e roxos. Para cada gene, os valores de expressão de cada amostra em duas linhagens diferentes de camundongos (C57BL/6J e DBA/2J) são plotados. Ambos os genes têm a mesma média para os dois grupos de amostras, mas o gene verde tem pouca variação dentro do grupo, enquanto o gene roxo tem um nível alto de variação. Para o gene verde com baixa variação dentro do grupo, a estimativa de LFC não retraída (vértice da linha verde sólida) é muito semelhante à estimativa de LFC retraída (vértice da linha verde tracejada). No entanto, devido à alta dispersão, a estimativa de LFC para o gene roxo é muito diferente. Portanto, mesmo que dois genes possam ter valores de contagem normalizados semelhantes, podem ter diferentes graus de retração de LFC. Observe que as estimativas de LFC são retraídas em direção ao priori (linha sólida preta).

A retração da mudança de dobramento em log2 não altera o número total de genes identificados como diferencialmente expressos significativos. A retração de fold change é para ajudar na avaliação posterior dos resultados. Por exemplo, se você deseja criar um subconjunto de genes importantes com base na mudança de dobramento para avaliação adicional, pode deseja usar os valores retraídos. Além disso, para ferramentas de análise funcional como o GSEA, que requerem valores de fold change como entrada, você pode querer fornecer os valores retraídos.

Para gerar as estimativas retraídas de log2 fold change, você deve executar uma etapa adicional usando a função lfcShrink() no objeto de resultados (que criaremos abaixo).

res_oe_original <- res_oe

# Aplicar retração de fold change
res_oe <- lfcShrink(dds, coef="sampletype_MOV10_overexpression_vs_control", type="apeglm")

Dependendo da versão do DESeq2 utilizada, o método padrão para estimativas retraídas varia. Conforme mostrado acima, o padrão pode ser alterado adicionando o parâmetro de tipo à função lfcShrink(). Para as versões mais recentes do DESeq2, type="normal" é o padrão e era o único método em versões anteriores. Demonstrou-se que existem métodos alternativos com menos viés do que o método "normal" na maioria dos casos, portanto optamos por usar o apeglm.

Gráfico MA

Um gráfico disponível para explorar nossos resultados é o gráfico MA. O gráfico MA mostra a média das contagens normalizadas versus a mudança de dobramento em log2 para todos os genes testados. Genes significativamente DE são coloridos para facilitar a identificação. Isso também é uma boa maneira de ilustrar o efeito da retração de LFC. O pacote DESeq2 fornece uma função simples para gerar gráficos MA.

Vamos começar com os resultados não retraídos:

plotMA(res_oe_original, ylim=c(-2,2))

Resultados retraídos

plotMA(res_oe, ylim=c(-2,2))

No lado esquerdo, você tem os valores de fold change não retraídos, e pode ver alta dispersão para genes com baixa expressão. Isto é, muitos com baixa expressão mostram mudanças de dobramento muito altas. Após a retração, vemos que as estimativas de mudança de dobramento são muito menores.

Além da comparação acima, o gráfico também nos permite avaliar a magnitude das mudanças de dobramento e sua distribuição em relação à expressão média. Geralmente, queremos ver genes significativos em toda a faixa de níveis de expressão.