Modelos Lineares Generalizados e Plots com edgeR – Análise de Expressão Diferencial Avançada

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[This article was first published on R – Myscape, and kindly contributed to R-bloggers]. (Você pode relatar problemas sobre o conteúdo desta página aqui)


Quer compartilhar seu conteúdo em R-bloggers? clique aqui se você tiver um blog, ou aqui se não tiver.

Modelos lineares generalizados (GLM) são um método clássico para analisar dados de expressão de RNA-seq. Em contraste com os testes exatos, os GLMs permitem comparações mais gerais.

Gráficos de exemplo das etapas de preparação de dados (linha superior) e resultados filtrados de testes F de quase verossimilhança (linha inferior) para a amostra Daphnia conjunto de dados.

Os tipos de comparações que você pode fazer dependerão do projeto do seu estudo. No exemplo a seguir, usaremos as contagens brutas de genes expressos diferencialmente (DE) para comparar o seguinte Daphnia genótipos.

O projeto de nosso experimento é descrito por três repetições de tratamento com radiação ultravioleta (UV) e três repetições de controle de luz visível (VIS) para cada um dos Daphnia genótipos.

Após a normalização da contagem de genes brutos, realizaremos testes F de quase-verossimilhança (QL) usando GLMs em edgeR. Esses testes nos permitem realizar testes do tipo ANOVA para diferenças nas médias dentro e entre nossos fatores de agrupamento.

As hipóteses que iremos testar:

  1. Os meios de tratamento fator são iguais
  2. Os meios da tolerância fator são iguais
  3. Não há interação entre os dois fatores

Observe que, optamos por usar testes QL F em vez de testes de razão de verossimilhança (LRT), uma vez que mantém a incerteza na estimativa da dispersão para cada gene.

D. melanica das montanhas olímpicas em WA (à esquerda) e das montanhas de Sierra Nevada na Califórnia (à direita).

QL F-Tests e script de plotagem – glmQLFTest_edgeR.r

O seguinte script R será usado para preparar contagens de genes brutos para testes de QL F em edgeR. Observe que, ao trabalhar com leituras de RNA-seq, você primeiro precisa realizar:

  1. Controle de qualidade – FastQC
  2. Corte – Trimmomatic
  3. Alinhamento – Hisat2
  4. Classificação – Samtools
  5. Quantificação –HTSeq

Como alternativa para alinhar as leituras a um genoma de referência, você pode desejar fazer uma montagem de novo de suas leituras usando um programa como o Trinity. Observe também que existem muitos programas disponíveis para realizar as etapas anteriores. Os programas que você escolher dependerão de seus dados e organismo de estudo.

Depois que os dados foram processados, podemos realizar testes QL F usando GLMs em edgeR. Usaremos as seguintes funções de plotagem para gerar gráficos informativos de nossos conjuntos de dados de entrada e ANOVA como resultados de teste com as bibliotecas edgeR e statmod.

  • trama de barra – Gráficos de barras de tamanhos de biblioteca de leitura de sequência
  • plotMDS – Gráfico de dispersão de escala multidimensional (MDS) de distâncias entre as amostras
  • plotBCV – Gráficos do coeficiente de variação biológica genômica (BCV) contra a abundância genética
  • plotMD – Gráficos de diferença média (MD) de todas as alterações logarítmicas (logFC) em relação ao tamanho médio da contagem
  • plotQLDisp – Gráfico da dispersão de QL genewise contra a abundância do gene (log2 CPM)

Entradas de script

Um arquivo csv formatado de contagens brutas de genes é esperado como a primeira entrada para o script. Abaixo está um subconjunto de exemplo de uma tabela de contagem de genes bruta, onde colunas representam amostras e linhas são genes.

A segunda entrada para o script deve ser um arquivo csv que descreve o projeto experimental e os fatores de agrupamento para cada amostra. Abaixo está um arquivo csv de exemplo que descreve o projeto experimental de nosso estudo.


R Script

Primeiro, devemos adicionar alguns comentários no início do script R descrevendo como executar o script. Devemos também adicionar um comentário com um lembrete do propósito do script.

#!/usr/bin/env Rscript
#Usage: Rscript glmQLFTest_edgeR.r rawGeneCountsFile experimentalDesign
#Usage Ex: Rscript glmQLFTest_edgeR.r daphnia_rawGeneCounts_htseq.csv daphnia_experimentalDesign.csv
#R script to perform QL F-tests using generalized linear models in edgeR

Antes de prosseguir, precisaremos instalar as bibliotecas edgeR e statmod. Essas bibliotecas precisarão ser instaladas apenas uma vez.

#Install edgeR and statmod, this should only need to be done once
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("edgeR")
install.packages("statmod")

As bibliotecas que estamos usando precisarão ser carregadas cada vez que executarmos o script usando o biblioteca função.

#Load the edgeR and statmod libraries
library("edgeR")
library("statmod")

Em seguida, recuperamos os caminhos de entrada dos arquivos csv para as contagens brutas de genes e o projeto experimental.

#Retrieve inputs to the script
args = commandArgs(trailingOnly=TRUE)

Agora podemos importar nossos dados dos caminhos do arquivo de entrada usando o read.csv função wrapper. Se os nomes das linhas para as contagens de genes tiverem um cabeçalho, ele precisará ser especificado com o row.names argumento.

Leia Também  Usando Python para trapacear no Scrabble
cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
#Import gene count data and view the first few rows
countsTable 

The grouping factors need to be added as a column to our experimental design before we can create a DGE list object. Here we use paste and factor to combine the factors (treatment and genotype) into one string (group) separated by a period for each sample. These strings are stored as a vector in the group data frame.

The relationship of the grouping factors to our experimental design.

Using the imported raw gene counts and grouping factors, we can create a DGE list object that has the grouping factors assigned to the gene count samples using the DGEList function.

#Setup and view the grouping factors
group 

As a first step in the analysis, we will plot the library sizes of our sequencing reads before normalization using the barplot function. The resulting plot is saved to the glmQLF_plotLibrarySizes jpg file.

#Plot the library sizes before normalization
jpeg("glmQLF_plotLibrarySizes.jpg")
barplot(list$samples$lib.size*1e-6, names=1:ncol(list), ylab="Library size (millions)")
dev.off()

Em seguida, precisamos filtrar as contagens de genes brutos por níveis de expressão e remover contagens de genes expressos de forma inferior.

#Retain genes only if it is expressed at a minimum level
keep 

The filtered raw counts are then normalized with calcNormFactors according to the weighted trimmed mean of M-values (TMM) to eliminate composition biases between libraries. The normalized gene counts in counts per million (CPM) are output to the glmQLF_normalizedCounts csv file using the write.table function.

#Use TMM normalization to eliminate composition biases
list 

Now that we have normalized gene counts for our samples we should generate some informative plots of our data.

First, we can verify the TMM normalization with a mean difference (MD) plot of all log fold change (logFC) against average count size. The resulting plot is saved to the glmQLF_plotNormalizedMD jpg file.

#Verify TMM normalization using a MD plot
jpeg("glmQLF_plotNormalizedMD.jpg")
plotMD(cpm(list, log=TRUE), column=1)
abline(h=0, col="red", lty=2, lwd=2)
dev.off()

Em seguida, vamos usar plotMDS para exibir as semelhanças relativas das amostras e ver as diferenças entre os perfis de expressão de diferentes amostras. O gráfico resultante é salvo no glmQLF_plotMDS arquivo jpg.

Uma legenda é adicionada ao canto superior esquerdo do gráfico, especificando "topleft" no lenda comando. A localização da legenda pode ser alterada atualizando o primeiro argumento do comando.

Observe que os quadros de dados do pontos e cores precisará ser atualizado se você tiver um número diferente de conjuntos de amostra do que em nosso exemplo. Esses quadros de dados são usados ​​para definir o pch e col argumentos do plotMDS função.

Aqui temos 2 Daphnia conjuntos de amostras de genótipos, então o pontos para cada amostra são especificados como os seguintes pares:

  1. 0 para quadrado vazio (tratamento Y023 UV), 15 para quadrado preenchido (tratamento Y023 VIS)
  2. 1 para círculo vazio (tratamento UV Y05), 16 para círculo preenchido (tratamento VIS Y05)

o cores para cada par de pontos de amostra é especificado da seguinte forma:

  1. azul (genótipo Y023)
  2. vermelho (genótipo Y05)
#Use a MDS plot to visualizes the differences
 between samples
#Set the point shapes and colors
points 

The design matrix for our data also needs to be specified before we can perform the F-tests. The experimental design is parametrized with a one-way layout and one coefficient is assigned to each group.

The design matrix for our experiment, which we can verify by comparing to our input gene counts.
#Setup the design matrix
design 

With the normalized gene counts and design matrix we can now generate the negative binomial (NB) dispersion estimates using the estimateDisp function. The NB dispersion estimates reflect the overall biological variability under the QL framework in edgeR.

This allows us to use the plotBCV function to generate a genewise biological coefficient of variation (BCV) plot of dispersion estimates. The resulting plot is saved to the glmQLF_plotBCV jpg file.

#Generate the NB dispersion estimates
list 

Next, we estimate the QL dispersions for all genes using the glmQLFit function. This detects the gene-specific variability above and below the overall level. The dispersion are then plotted with plotQLDisp, and the resulting plot is saved to the glmQLF_plotQLDisp jpg file.

#Estimate and view the QL dispersions
fit 

Now we are ready to begin defining and testing contrasts of our experimental design. The first comparison we will make is used to test our hypothesis that the means of the treatment factor are equal.

So, we will make a contrast that tests whether the average across all UV groups is equal to the average across all VIS groups using makeContrasts. Note that UV and VIS are the two levels in our treatment group.

#Design a contrast to test overall effect of the first factor
con.UVvsVIS 

With the contrast defined we can use the glmQLFTest function to test the effect of our first factor. The DE genes in our test are plotted using plotMD and written to the glmQLF_UVvsVIS_plotMD jpg.

The top tags for the DE genes that are above a 0.05 false discovery rate (FDR) cutoff are written to the glmQLF_UVvsVIS_topTags csv. A FDR of 0.05 is used to filter the results to the most biologically relevant genes.

#Look at genes expressed across all UV groups
test.anov.one 

We can use the glmTreat function to filter out genes that have a low logFC and are therefore less significant to our analysis. In this example we will use a log2 FC cutoff of 1.2, since we see a relatively large number of DE genes from our test.

The DE genes from our filtered test results are plotted using plotMD and written to the glmQLF_UVvsVIS_plotMD_filtered jpg. The top tags for the DE genes are selected using the topTags function, which default ranks the genes by p-value. The top tags that are above a 0.05 false discovery rate (FDR) cutoff are written to the glmQLF_UVvsVIS_topTags_filtered csv.

#Look at genes with significant expression across all UV groups
treat.anov.one 

The second comparison we will make is used to test our hypothesis that the means of the tolerance factor are equal.

So, we will need to make a second contrast that tests whether the average across the Tolerant group (Y05 genotype) is equal to the average across the Not Tolerant group (Y023 genotype).

#Design a contrast to test overall effect of the first factor
con.TvsN 

With the second contrast defined we again use glmQLFTest to test the effect of our second factor. The DE genes in our second test are plotted using plotMD and written to the glmQLF_TvsN_plotMD jpg. The top tags for the DE genes that are above a 0.05 false discovery rate (FDR) cutoff are written to the glmQLF_TvsN_topTags csv.

#Look at genes expressed across all tolerant groups
test.anov.two 

Again we use glmTreat to filter out genes that have a low logFC. Then the DE genes from our filtered test results are plotted using plotMD and written to the glmQLF_TvsN_plotMD_filtered jpg. The top tags that are above 0.05 FDR are written to the glmQLF_TvsN_topTags_filtered csv.

#Look at genes with significant expression across all tolerant groups
treat.anov.two 

The final comparison we will make is to test our hypothesis that there is no interaction between the two factors (treatment and tolerance). To make our last contrast we will test whether the mean effects of the two factors are equal.

#Test whether there is an interaction effect
con.Interaction 

With the last contrast defined we use glmQLFTest to test the interaction effect. The DE genes in our final test are saved to the glmQLF_Interaction_plotMD jpg. The top tags above 0.05 FDR are written to the glmQLF_Interaction_topTags csv.

#Look at genes expressed across all tolerant groups
test.anov.Interaction 

Finally, the DE genes from our last filtered test results are saved to the glmQLF_Interaction_plotMD_filtered jpg. The top tags that are above 0.05 FDR are written to the glmQLF_Interaction_topTags_filtered csv.

#Look at genes with significant expression across all tolerant groups
treat.anov.Interaction 


Example Outputs – glmQLFTest_edgeR.r

The following plots are example outputs from the above R script.

Example plots from data preparation steps, from left to right: bar plot of library sizes, MD plot of normalized counts, MDS plot of samples, BCV plot of dispersion estimates, and QL dispersion plot of dispersions for all genes.
Example plots of quasi-likelihood F-tests for the sample Daphnia data set. The top row of plots show the differentially expressed genes in each of our hypotheses before filtering. The bottom row of plots show the differentially expressed genes in each of our hypotheses after filtering.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br