Três métodos diferenciais de análise de expressão para sequenciamento de RNA: limma, EdgeR, DESeq2

Resumo

Foi fornecido um protocolo detalhado de métodos de análise de expressão diferencial para sequenciamento de RNA: limma, EdgeR, DESeq2.

Resumo

O sequenciamento de RNA (RNA-seq) é uma das tecnologias mais utilizadas em transcriptômica, pois pode revelar a relação entre a alteração genética e os processos biológicos complexos e tem grande valor em diagnósticos, prognósticos e terapêuticas de tumores. A análise diferencial dos dados do RNA-seq é crucial para identificar transcrições aberrantes, e limma, EdgeR e DESeq2 são ferramentas eficientes para análise diferencial. No entanto, a análise diferencial do RNA-seq requer certas habilidades com linguagem R e a capacidade de escolher um método adequado, que está faltando no currículo da educação médica.

Aqui, fornecemos o protocolo detalhado para identificar genes expressos diferencialmente (DEGs) entre o cholangiocarcinoma (CHOL) e os tecidos normais através de limma, DESeq2 e EdgeR, respectivamente, e os resultados são mostrados em parcelas vulcânicas e diagramas de Venn. Os três protocolos de limma, DESeq2 e EdgeR são semelhantes, mas têm etapas diferentes entre os processos da análise. Por exemplo, um modelo linear é usado para estatísticas em limma, enquanto a distribuição binomial negativa é usada em edgeR e DESeq2. Além disso, os dados de contagem de RNA-seq normalizados são necessários para EdgeR e limma, mas não é necessário para o DESeq2.

Aqui, fornecemos um protocolo detalhado para três métodos de análise diferencial: limma, EdgeR e DESeq2. Os resultados dos três métodos são parcialmente sobrepostos. Todos os três métodos têm suas próprias vantagens, e a escolha do método depende apenas dos dados.

Introdução

O RNA-sequenciamento (RNA-seq) é uma das tecnologias mais utilizadas em transcriptômicas com muitas vantagens (por exemplo, reprodutibilidade de dados elevados), e aumentou drasticamente nossa compreensão das funções e dinâmicas dos processos biológicos complexos^1,2. A identificação de transcrições aberratas em diferentes contextos biológicos, também conhecidos como genes expressos diferencialmente (DEGs), é um passo fundamental na análise do RNA-seq. O RNA-seq torna possível obter uma compreensão profunda dos mecanismos moleculares relacionados à patogênese e funções biológicas. Portanto, a análise diferencial tem sido considerada valiosa para diagnósticos, prognósticos e terapêuticas de tumores^3,4,5. Atualmente, mais pacotes de R/Bioconductor de código aberto foram desenvolvidos para análise de expressão diferencial RNA-seq, particularmente limma, DESeq2 e EdgeR^1,6,7. No entanto, a análise diferencial requer certas habilidades com linguagem R e a capacidade de escolher o método adequado, que está faltando no currículo da educação médica.

Neste protocolo, com base nos dados da contagem de RNA-seq de cholangiocarcinoma (CHOL) extraídos do Atlas do Genoma do Câncer (TCGA), três dos métodos mais conhecidos (limma⁸, EdgeR⁹ e DESeq2¹⁰) foram realizados, respectivamente, pelo programa R¹¹ para identificar os DEGs entre CHOL e tecidos normais. Os três protocolos de limma, EdgeR e DESeq2 são semelhantes, mas têm etapas diferentes entre os processos de análise. Por exemplo, os dados de contagem de RNA-seq normalizados são necessários para EdgeR e limma^8,9, enquanto o DESeq2 usa suas próprias discrepâncias de biblioteca para corrigir dados em vez de normalização¹⁰. Além disso, edgeR é especificamente adequado para dados RNA-seq, enquanto o limma é usado para microarrays e RNA-seq. Um modelo linear é adotado pela limma para avaliar os DEGs^12,enquanto as estatísticas em edgeR são baseadas nas distribuições binomiais negativas, incluindo estimativa empírica de Bayes, testes exatos, modelos lineares generalizados e testes de quase-probabilidade⁹.

Em resumo, fornecemos os protocolos detalhados da análise de expressão diferencial RNA-seq utilizando limma, DESeq2 e EdgeR, respectivamente. Ao se referir a este artigo, os usuários podem facilmente realizar a análise diferencial do RNA-seq e escolher os métodos de análise diferencial adequados para seus dados.

Protocolo

NOTA: Abra o programa R-studio e carregue o arquivo R "DEGs.R", o arquivo pode ser adquirido a partir de arquivos suplementares/Scripts.

1. Download e pré-processamento de dados

1. Baixe os dados da contagem de colhidos (HTSeq) de cholangiocarcinoma (CHOL) do Atlas do Genoma do Câncer (TCGA). Esta etapa pode ser facilmente alcançada pelo seguinte código R.
   1. Clique em Executar para instalar pacotes R.
   2. Clique em Executar para carregar pacotes R.
      se(!requeroospace("BiocManager", silenciosamente=TRUE))
      + install.packages("BiocManager")
      BiocManager::install(c("TCGAbiolinks", "SummarizedExperiment"))
   3. Defina o diretório de trabalho.
      biblioteca (TCGAbiolinks)
      biblioteca (SummarizedExperiment)
      setwd ("C:/Users/LIUSHIYI/Desktop")
   4. Escolha o tipo de câncer.
      câncer <- "TCGA-CHOL"
   5. Execute o código R do arquivo "GDCquery.R" para baixar os dados. O arquivo "GDCquery.R" pode ser adquirido a partir de arquivos suplementares/Scripts:
      fonte ("Arquivos suplementares/Scripts/GDCquery.R")
      cabeça(cnt)
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      ##ENSG00000000003 4262
      ##ENSG00000000005 1
      ##ENSG00000000419 1254
      ##ENSG00000000457 699
      ##ENSG00000000460 239
      ##ENSG00000000938 334
      NOTA: Após a execução, os dados da contagem de CHOLHTSeq serão baixados e nomeados "cnt", onde as linhas representam IDs e colunas genéticas do conjunto representam IDs de amostra. Observe os números nas posições 14-15 nos IDs da amostra; números que variam de 01 a 09 indicam tumores e variando de 10 a 19 indicam tecidos normais.
2. Converta os genes do conjunto em símbolos genéticos.
   1. Importe o arquivo de anotação em R de acordo com seu caminho de armazenamento. O arquivo de anotação (gencode.v22.anotação.gtf) pode ser adquirido a partir de arquivos suplementares.
      gtf_v22 <-rtracklayer::import('Arquivos suplementares/gencode.v22.anotação.gtf')
   2. Execute o código R do "gtf_v22. Arquivo R", que pode ser adquirido a partir de arquivos suplementares/Scripts:
      fonte ("Arquivos suplementares/Scripts/gtf_v22. R")
   3. Aplique a função "ann" para converter os IDs genéticos do conjunto em símbolos genéticos.
      cnt=ann(cnt,gtf_v22)
3. Filtrando genes de baixa expressa
   1. Clique em Executar para instalar o pacote R "edgeR".
      BiocManager::install("edgeR")
   2. Clique em Executar para carregar o pacote R "edgeR".
      biblioteca (edgeR)
   3. Execute o seguinte código R para manter os genes com valores de contagem por milhão (CPM) maiores que um em pelo menos duas amostras.
      manter <-rowSums(cpm(cnt)>1)>=2
      cnt <- as.matrix(cnt[keep,])
      NOTA: O valor de contagem por milhão (CPM) é usado em vez da contagem de leitura para eliminar o desvio causado por diferentes profundidades de sequenciamento.

2. Análise de expressão diferencial através de "limma"

1. Clique em Executar para instalar o pacote R "limma".
   BiocManager::install("limma")
2. Clique em Executar para carregar os pacotes R "limma", "edgeR".
   biblioteca (limma)
   biblioteca (edgeR)
3. Execute o seguinte código R para criar a matriz de design.
   grupo <- substring(colnames(cnt),14,15) # Extract group information
   grupo [grupo %em% "01"] <- "Cancer" # set '01' as tumor tissue
   grupo [grupo %em% "11"] <- "Normal" # set '11' as normal tissue
   grupo <- factor (group, levels = c("Normal","Cancer"))
   1. Crie a matriz de design.
      design <-model.matrix (~grupo)
      nomes de linha (design) <-colnames(cnt)
   2. Crie o objeto DGEList.
      dge <- DGEList(contagens = cnt, grupo = grupo)
   3. Normalize os dados.
      dge <-calcNormFactors(dge, método = "TMM")
   4. Execute o seguinte código R para executar a análise de expressão diferencial baseada no método limma-trend.
      dge
      ##An objeto da classe "DGEList"
      ##$counts
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      ##TSPAN6 4262
      ##DPM1 1254
      ##SCYL3 699
      ##C1orf112 239
      ##FGR 334
   5. Calcule o valor do CPM.
      logdge <-cpm(dge, log=TRUE, prior.count=3)
   6. Clique em Executar para encaixar em um modelo linear para prever os dados ou inferir a relação entre variáveis.
      caber <-lmFit (logdge, design)
   7. Calcule o valor T, o valor F e as probabilidades de log com base no Bayesian.
      ajuste <-eBayes (fit, trend=TRUE)
   8. Extrair a tabela de resultados.
      res_limma<-as.data.frame (topTable(fit,n=Inf))
      
      cabeça (res_limma)
      ## logFC AveExpr t P.Value adj. P.Val B
      ##RP11-252E2.2 -4.899493 -2.488589 -20.88052 2.386656e-25 4.931786e-21 47.28823
      ##BX842568.1 -4.347930 -2.595205 -20.14532 1.082759e-24 1.118706e-20 45.83656
      ##CTC-537E7.3 -5.154894 -2.143292 -19.59571 3.452354e-24 2.216114e-20 44.72001
      ##RP11-468N14.3 -6.532259 -2.029714 -19.49409 4.289807e-24 2.216114e-20 44.51056
      ##AP006216.5 -4.507051 -2.670915 -19.25649 7.153356e-24 2.956339e-20 44.01704
      ##RP11-669E14.4 -4.107204 -2.828311 -18.93246 1.448209e-23 4.987633e-20 43.33543
      #The resultado da análise de expressão diferencial é salvo em "res_limma", que inclui o id genético, o valor de alteração do fold log2 (logFC), o nível médio de expressão log2 do gene no experimento (AveExpr), a estatística t modificada (t), o valor de relavent p (P.Value), a taxa de descoberta falsa (FDR) corrigida p valor (adj). P.Val) e as probabilidades de log de genes expressos diferencialmente (B)
      NOTA: A função "calcNormFactors()" do "edgeR" foi usada para normalizar os dados para eliminar a influência causada pela preparação da amostra ou construção e sequenciamento da biblioteca. Na construção da matriz de design, é necessário combinar o design experimental (por exemplo, tipo de tecido: tecidos normais ou tumorais) para amostrar IDs da matriz. limma-trend é adequado para dados cuja profundidade de sequenciamento é a mesma, enquanto limma-voom é adequado: (i) quando o tamanho da biblioteca de amostras é diferente; (ii) dados não normalizados pela TMM; (iii) há muito "ruído" nos dados. Um logFC positivo significa que o gene é regulado no experimento, enquanto o número negativo significa que o gene é regulado.
   9. Identifique os DEGs.
      res_limma$sig <- como.factor(
      ifelse (res_limma$adj. P.Val < 0,05 & abs (res_limma$logFC) > 2,
      ifelse (res_limma$logFC > 2 'up'',''down'),''')) # O valor adj.p < 0,05 e o |log2FC| >= 2 são limiares para identificar os DEGs
      resumo (res_limma$sig)
      ##down não para cima
      ##1880 ​17341 1443
   10. Saída a tabela de resultados para um arquivo.
       escrever.csv(res_limma, arquivo = 'result_limma.csv')
   11. Clique em Executar para instalar o pacote R "ggplot2".
       install.packages("ggplot2")
   12. Clique em Executar para carregar o pacote R "ggplot2".
       biblioteca(ggplot2)
   13. Execute o código R do vulcão. R" para criar o enredo do vulcão. O arquivo "vulcão. R" pode ser adquirido a partir de arquivos complementares.
       fonte ("Arquivos suplementares/Scripts/volcano. R")
       vulcão (res_limma,"logFC","adj. P.Val",2,0,05)
       NOTA: Os genes podem ser mapeados para diferentes posições de acordo com seus valores log2FC e adj-p, os DEGs regulamentados para cima são coloridos em vermelho e os DEGs regulamentados são coloridos em verde.
   14. Clique em Exportar para salvar o enredo do vulcão.
       NOTA: As parcelas do vulcão podem ser geradas e baixadas em diferentes formatos (por exemplo, pdf, TIFF, PNG, formato JPEG). Os genes podem ser mapeados para diferentes posições de acordo com seus valores log2FC e adj p, os DEGs up-regulated (log2FC > 2, adj p < 0,05) são coloridos em vermelho, e os DEGs não regulamentados (log2FC < -2, adj p < 0,05) são coloridos em verde, não-DEGs são coloridos em cor cinza.

3. Análise diferencial da expressão através do "edgeR"

1. Clique em Executar para carregar o pacote R "edgeR".
   biblioteca (edgeR)
2. Execute o seguinte código R para criar matriz de design.
   grupo <-substring (colnames(cnt),14,15)
   grupo [grupo %em% "01"] <- "Câncer"
   grupo [grupo %em% "11"] <- "Normal"
   group=fator (grupo, níveis = c("Normal","Câncer"))
   design <-model.matrix(~grupo)
   nomes de linha (design) = colnames(cnt)
3. Clique em Executar para criar o objeto DGEList.
   dge <-DGEList (counts=cnt)
4. Normalize os dados.
   dge <-calcNormFactors(dge, método = "TMM")
5. Clique em Executar para estimar a dispersão dos valores de expressão genética.
   dge <- estimativaDisp(dge, design, robusto = T)
6. Clique em Executar para ajustar o modelo para contar dados.
   caber <-glmQLFit(dge, design)
7. Faça um teste estatístico.
   caber <-glmQLFTest (fit)
8. Extrair a tabela de resultados. O resultado é salvo em "res_edgeR", que inclui o valor de alteração do fold de log, o cpm de registro, f, valor p e valor p corrigido fdr.
   res_edgeR=as.data.frame(topTags(fit, n=Inf))
   cabeça (res_edgeR)
   ## logFC logCPM F PValue FDR
   ##GCDH -3.299633 5.802700 458.5991 1.441773e-25 2.979280e-21
   ##MSMO1 -3.761400 7.521111 407.0416 1.730539e-24 1.787993e-20R
   ##CL1 -3.829504 5.319641 376.5043 8.652474e-24 5.516791e-20
   ##ADI1 -3,533664 8.211281 372.6671 1.067904e-23 5.516791e-20
   ##KCNN2 -5.583794 3.504017 358.6525 2.342106e-23 9.679455e-20
   ##GLUD1 -3.287447 8.738080 350.0344 3.848408e-23 1.194406e-19
   #The resultado é salvo em "res_edgeR", que inclui o valor de alteração da dobra de log (logFC), cpm de registro, F, valor p e valor p corrigido por FDR
9. Identifique os DEGs.
   res_edgeR$sig = as.factor(
   ifelse (res_edgeR$FDR < 0,05 & abs (res_edgeR$logFC) > 2,
   ifelse (res_edgeR$logFC > 2 'up','down'),'''))
   resumo (res_edgeR$sig)
   ##down não para cima
   ##1578 15965 3121
10. Saída a tabela de resultados para um arquivo.
    escrever.csv(res_edgeR, arquivo = 'res_edgeR.csv')
11. Crie o enredo do vulcão.
    vulcão (res_edgeR,"logFC","FDR",2,0,05)
12. Clique em Exportar para salvar o enredo do vulcão.

4. Análise de expressão diferencial através do "DESeq2"

1. Clique em Executar para instalar pacotes R "DESeq2".
   BiocManager::install("DESeq2")
2. Clique em Executar para carregar pacotes R "DESeq2".
   biblioteca (DESeq2)
3. Execute o seguinte código R para determinar o fator de agrupamento.
   grupo <-substring (colnames(cnt),14,15)
   grupo [grupo %em% "01"] <- "Câncer"
   grupo [grupo %em% "11"] <- "Normal"
   group=fator (grupo, níveis = c("Normal","Câncer"))
4. Crie o objeto DESeqDataSet.
   dds <-DESeqDataSetFromMatrix (cnt, DataFrame(grupo), design = ~grupo)
   Dds
   ##class: DESeqDataSet
   ##dim: 20664 45
   ##metadata(1): versão
   ##assays(1): conta
   ##rownames(20664): TSPAN6 DPM1 ... RP11-274B21.13 LINC01144
   ##rowData nomes(0):
   ##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ...
   ##colData nomes(1): grupo
5. Faça a análise.
   dds <-DESeq(dds)
6. Gere a tabela de resultados.
   res_DESeq2 <-data.frame(resultados(dds))
   
   cabeça (res_DESeq2)
   ## baseMean log2FoldChange lfcSE stat pvalue padj
   ##TSPAN6 4704.9243 -0,8204515 0,3371667 -2.433370 1.495899e-02 2.760180e-02
   ##DPM1 1205.9087 -0,3692497 0.1202418 -3.070894 2.134191e-03 4.838281e-03
   ##SCYL3 954.9772 0.2652530 0.2476441 1.071106 2.841218e-01 3.629059e-01
   ##C1orf112 277.7756 0.7536911 0.2518929 2.992109 2.770575e-03 6.101584e-03
   ##FGR 345.8789 -0,6423198 0.3712729 -1.730047 8.362180e-02 1.266833e-01
   ##CFH 27982.3546 -3.8761382 0.5473363 -7.081823 1.422708e-12 1.673241e-11
   NOTA: O resultado é salvo em "res_DESeq2", que inclui a média da contagem de leitura normalizada (baseMean), valor de alteração do fold de log (log2FoldChange), erro padrão de alteração de dobra de log (lfcSE), a estatística de Wald (stat), valor p original (pvalue) e valor p corrigido (padj)
7. Identifique deGs.
   res_DESeq2$sig = as.factor(
   ifelse (res_DESeq2$padj < 0,05 & abs (res_DESeq2$log2FoldChange) > 2,
   ifelse (res_DESeq2$log2FoldChange > 2 'up','down'),'''))
   resumo (res_DESeq2$sig)
   ##down não para cima
   ##1616 16110 2938
8. Saída a tabela de resultados para um arquivo.
   escrever.csv(res_DESeq2, arquivo = 'res_DESeq2.csv')
9. Crie o enredo do vulcão.
   vulcão (res_DESeq2,"log2FoldChange","padj",2,0,05)
10. Clique em Exportar para salvar o enredo do vulcão.

5. Diagrama de Venn

1. Clique em Executar para instalar o pacote R "VennDiagram".
   install.packages("VennDiagram")
2. Clique em Executar para carregar o pacote R "VennDiagram".
   biblioteca (VennDiagram)
3. Faça um diagrama de Venn de DEGs regulamentados.
   grid.newpage()
   grid.draw (venn.diagram(lista(Limma=rownames(res_
   limma[res_limma$sig=="up",]),
   edgeR=rownames(res_edgeR[res_edgeR$sig=="up",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig==
   "para cima",])),
   NULO, altura = 3,largura = 3,unidades = "in",
   col="preto",lwd=0,3,fill=c ("#FF6666","#FFFF00",
   "#993366"),
   alfa=c (0,5, 0,5, 0,5), principal = "DEGs regulados"))
4. Clique em Exportar para salvar o diagrama de Venn.
5. Faça um diagrama de Venn de DEGs regulamentados.
   grid.newpage()
   grid.draw (venn.diagram(lista(Limma=rownames(res_
   limma[res_limma$sig=="down",]),
   edgeR=rownames(res_edgeR[res_edgeR$sig==
   "para baixo",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig="down",]),
   NULO, altura = 3,largura = 3,unidades = "in",
   col="preto",lwd=0,3,fill=c ("#FF6666","#FFFF00",
   "#993366"),
   alfa=c (0,5, 0,5, 0,5), principal = "DEGs não regulados"))
6. Clique em Exportar para salvar o diagrama de Venn.

Resultados

Existem várias abordagens para visualizar o resultado da análise de expressão diferencial, entre as quais o enredo do vulcão e o diagrama de Venn são particularmente utilizados. limma identificou 3323 DEGs entre o CHOL e tecidos normais com o |logFC|≥2 e adj. P.Val <0,05 como limiares, entre os quais 1880 foram regulados em tecidos CHOL e 1443 foram regulados(Figura 1a). Enquanto isso, o edgeR identificou os 1578 DEGs para baixo regulados e 3121 DEGs up-regulated(...

Discussão

Transcrições abundantes de aberrate em cânceres podem ser facilmente identificadas pela análise diferencial RNA-seq⁵. No entanto, a aplicação da análise de expressão diferencial RNA-seq é muitas vezes restrita, pois requer certas habilidades com linguagem R e a capacidade de escolher métodos apropriados. Para resolver esse problema, fornecemos uma introdução detalhada aos três métodos mais conhecidos (limma, EdgeR e DESeq2) e tutoriais para a aplicação da análise de expressão dif...

Materiais

Name	Company	Catalog Number	Comments
R		version 3.6.2	free software
Rstudio			free software