Normalização de dados em redes de co expressão - lmigueel/Bioinformatica GitHub Wiki

1. Sobre

Quem trabalha com Transcriptômica já deve ter pensando em realizar uma análise de co-expressão utilizando o pacote WGCNA. Com ele podemos identificar módulos de co-expressão, e linkar com meta dados também. Um pacote completo. Logo mais teremos um post sobre redes de co-expressão.

Mas logo no começo da análise, você deve normalizar os dados. O mais clássico de todos é utilizar a normalização pelo log. Mas será que sempre ela é a correta? Aqui vamos discutir isso.

Este post é focado na normalização de read counts para análises de co-expressão. Logo, caso você tenha os dados abaixo, segue as dias:

    1. Dados de TPM (Transcritos por milhão), RPKM ou FPKM:

Neste caso você deve aplicar a normalização inicial pelo log().

    1. Dados do GEO (Gene Expression Omnibus):

Neste caso você deve verificar se o dado já está normalizado ou não pelo plot da distribuição. Caso já esteja, sucesso, caso não, aplique o log(). Cuidado é sempre necessário para não ficar uma distribuição log-log.

2. Dados de RNA-Seq

Quando estou analisando dados de RNA-Seq, há duas propriedades específicas de RNA-Seq que sempre mantenho em mente:

  • A presença de valores extremos
  • A dependência média-variância (também conhecida como heterocedasticidade)

Essas duas propriedades são importantes a serem consideradas, dependendo da análise. Por exemplo, se você estiver fazendo uma análise diferencial de expressão gênica com DESeq ou edgeR, isso não é um problema. Nesse caso, você usará os read counts, sem qualquer transformação. O motivo pelo qual você não precisará transformar seus dados é porque DESeq e edgeR foram desenvolvidos especificamente para dados RNA-Seq e, portanto, cuidam desse tipo de propriedades internamente. No entanto, existem muitas outras aplicações de RNA-Seq que foram historicamente desenvolvidas com dados de Microarray e, portanto, são mais susceptíveis de serem influenciadas pela natureza altamente distorcida dos dados de RNA-Seq. Por exemplo, você pode querer considerar a transformação de dados ao construir assinaturas de genes prognósticos (leia Zwiener et al., 2014 para mais detalhes) ou ao construir redes co-regulatórias de genes, especialmente se usar a correlação de Pearson (porque os resultados de PC são fortemente afetados por outliers).

3. Valores extremos e distribuição enviesada

Na literatura, não é incomum ver valores de RNA-Seq representados na escala logarítmica. Na verdade, a transformação do logaritmo eliminará alguns valores extremos. Outro método frequentemente usado é a transformação de estabilização de variância (variance-stabilizing transformation - VST). O VST para dados RNA-Seq foi proposto por Anders e Huber, 2010 e está implementado no pacote DESeq.

Para verificar como o seu dado de expressão de comporta após as normalizações, você pode executar os comandos abaixo. Vou supor que você carregou a sua matriz de expressão em read counts com o nome de data. Você tem os read counts e deseja realizar análises posteriores de co-expressão. Logo, você deve normalizar os dados pelo log() ou pela função vst()? ** eu prefiro a vst() **.

A função vst() criada tem como padrão genérico o uso do DESeq, mas você deve alterar a variável 'condition' segundo a sua amostra. Um exemplo de dados para usar abaixo você pode baixar AQUI e apenas rodar o comando load("rnaseq_lusc_example_SeqQC.Rda") que ele já irá aplicar na variável data.

load("rnaseq_lusc_example_SeqQC.Rda")

vst <- function(countdata){
  library(DESeq)
  condition <- factor(rep("Tumour", ncol(countdata)))
  countdata <- newCountDataSet(countdata,condition )
  countdata <- estimateSizeFactors( countdata )
  cdsBlind <- DESeq::estimateDispersions( countdata, method="blind")
  vstdata <- varianceStabilizingTransformation( cdsBlind )
  return(exprs(vstdata))
}

data.log2 <- log2(data+1)
data.vst <- vst(data)

par(mfrow=c(1,3))
plot(density(as.numeric(data[1,])), main="counts", cex.main=2)
plot(density(as.numeric(data.log2[1,])), main="log2", cex.main=2)
plot(density(as.numeric(data.vst[1,])), main="VST", cex.main=2)

A figura terá três gráficos, o tipo de distribuição vai mudar completamente. Em quanto os dados brutos (variável data) possui uma binomial negativa, as outras se aproximam de uma distribuição normal.

4. Heteroscedasticidade

Em dados de RNA-Seq, genes com expressão média maior têm, em média, maiores variações observadas entre as amostras, isto é, eles variam em expressão de amostra para amostra mais do que outros genes com expressão média mais baixa. Esse fenômeno é conhecido como heterocedasticidade de dados. Heterocedasticidade é o oposto de homocedasticidade e significa literalmente "ter a dispersão diferente". Isso significa que os valores dos dados estão espalhados, ou espalhados, em diferentes extensões para diferentes genes.

Este efeito é claramente visível quando traçamos o desvio padrão por gene (obtido em N amostras), contra a classificação da expressão média. Você pode usar a função meanSdPlot do pacote vsn e demonstra a dependência da média-variância claramente. Após a transformação de log, ainda há variâncias desiguais entre as amostras, especialmente para as contagens mais baixas (mostrando um desvio padrão elevado). Após o VST, o dp torna-se mais constante ao longo de toda a faixa dinâmica, mas ainda são desiguais para todos os genes.

par(mfrow=c(1,3))
library(vsn)
meanSdPlot(as.matrix(data),ylim=c(0,8e3),main="counts")
meanSdPlot(as.matrix(data.log2), ylim = c(0,2.5),main="log2")
meanSdPlot(as.matrix(data.vst), ylim = c(0,2.5),main="VST")

5. Log vs. VST: qual você escolheria?

Depende! Uma coisa que normalmente faço é usar a transformação de log (quero dizer, o logaritmo deslocado log2(n + 1)) quando quero ter uma inspeção visual rápida de meus dados. Também faço muitas inferências de rede a partir de dados RNA-Seq e, para isso, uso o VST. Ao decidir entre o log e o VST, há algumas coisas a serem consideradas:

  • Após a transformação do log, há menos valores extremos quando comparados aos dados não transformados, mas ainda há variâncias desiguais

  • Após VST, o desvio padrão por gene torna-se mais constante ao longo de todo o intervalo dinâmico, mas observe que as variâncias ainda são desiguais para todos os genes

  • Um problema adicional da transformação log2 é que log2 de zero é infinito! Para evitar obter o logaritmo de zero, é comum adicionar um pseudo-valor de 1 antes de obter o log. E, é claro, temos que assumir que adicionar 1 não influencia muito as contagens baixas diferentes de zero

  • Para contagens mais baixas, a transformação do log é muito dramática. VST leva a menos valores variáveis ​​em toda a faixa dinâmica.

Em resumo, ambas são técnicas de transformação populares e a decisão depende de sua aplicação. O log é simples e fácil de interpretar, mas não permite que você tire o log do zero e pode não corrigir a heterocedasticidade completamente. VST pode ser a solução preferida para sua aplicação.

Referências

Zwiener, I., Frisch, B., & Binder, H. (2014). Transforming RNA-Seq data to improve the performance of prognostic gene signatures. PloS one, 9(1), e85150.

Anders, S., & Huber, W. (2010). Differential expression analysis for sequence count data. Nature Precedings, 1-1.