Tres métodos de análisis de expresión diferencial para la secuenciación de ARN: limma, EdgeR, DESeq2

Resumen

Se proporcionó un protocolo detallado de métodos de análisis de expresión diferencial para la secuenciación de ARN: limma, EdgeR, DESeq2.

Resumen

La secuenciación de ARN (RNA-seq) es una de las tecnologías más utilizadas en transcriptómica, ya que puede revelar la relación entre la alteración genética y los procesos biológicos complejos y tiene un gran valor en el diagnóstico, pronóstico y terapéutica de tumores. El análisis diferencial de los datos de RNA-seq es crucial para identificar transcripciones aberrantes, y limma, EdgeR y DESeq2 son herramientas eficientes para el análisis diferencial. Sin embargo, el análisis diferencial RNA-seq requiere ciertas habilidades con el lenguaje R y la capacidad de elegir un método apropiado, que falta en el currículo de la educación médica.

Aquí, proporcionamos el protocolo detallado para identificar genes expresados diferencialmente (DEG) entre el colangiocarcinoma (CHOL) y los tejidos normales a través de limma, DESeq2 y EdgeR, respectivamente, y los resultados se muestran en gráficos de volcanes y diagramas de Venn. Los tres protocolos de limma, DESeq2 y EdgeR son similares pero tienen diferentes pasos entre los procesos del análisis. Por ejemplo, se utiliza un modelo lineal para la estadística en limma, mientras que la distribución binomial negativa se utiliza en edgeR y DESeq2. Además, los datos de recuento normalizado de RNA-seq son necesarios para EdgeR y limma, pero no son necesarios para DESeq2.

Aquí, proporcionamos un protocolo detallado para tres métodos de análisis diferencial: limma, EdgeR y DESeq2. Los resultados de los tres métodos se superponen en parte. Los tres métodos tienen sus propias ventajas, y la elección del método solo depende de los datos.

Introducción

La secuenciación de ARN (RNA-seq) es una de las tecnologías más utilizadas en transcriptómica con muchas ventajas (por ejemplo, alta reproducibilidad de datos), y ha aumentado drásticamente nuestra comprensión de las funciones y dinámicas de procesos biológicos complejos^1,2. La identificación de transcripciones aberrantes bajo diferentes contextos biológicos, que también se conocen como genes expresados diferencialmente (DEG), es un paso clave en el análisis de ARN-seq. RNA-seq permite obtener una comprensión profunda de los mecanismos moleculares relacionados con la patogénesis y las funciones biológicas. Por lo tanto, el análisis diferencial ha sido considerado como valioso para el diagnóstico, pronóstico y terapéutica de tumores^3,4,5. Actualmente, se han desarrollado más paquetes R/Bioconductor de código abierto para el análisis de expresión diferencial RNA-seq, particularmente limma, DESeq2 y EdgeR^1,6,7. Sin embargo, el análisis diferencial requiere ciertas habilidades con el lenguaje R y la capacidad de elegir el método apropiado, que falta en el plan de estudios de la educación médica.

En este protocolo, basado en los datos de conteo de ARN-seq de colangiocarcinoma (CHOL) extraídos del Atlas del Genoma del Cáncer (TCGA), se llevaron a cabo tres de los métodos más conocidos (limma^8,EdgeR⁹ y DESeq2^10),respectivamente, por el programa R¹¹ para identificar los DEG entre CHOL y tejidos normales. Los tres protocolos de limma, EdgeR y DESeq2 son similares pero tienen diferentes pasos entre los procesos del análisis. Por ejemplo, los datos de recuento normalizado de RNA-seq son necesarios para EdgeR y limma^8,9, mientras que DESeq2 utiliza sus propias discrepancias de biblioteca para corregir datos en lugar de la normalización¹⁰. Además, edgeR es específicamente adecuado para datos de RNA-seq, mientras que el limma se utiliza para microarrays y RNA-seq. Limma adopta un modelo lineal para evaluar los DEGs^12,mientras que las estadísticas en edgeR se basan en las distribuciones binomiales negativas, incluyendo estimación empírica de Bayes, pruebas exactas, modelos lineales generalizados y pruebas de cuasi-verosimilitud⁹.

En resumen, proporcionamos los protocolos detallados de análisis de expresión diferencial RNA-seq mediante el uso de limma, DESeq2 y EdgeR, respectivamente. Al consultar este artículo, los usuarios pueden realizar fácilmente el análisis diferencial RNA-seq y elegir los métodos de análisis diferencial apropiados para sus datos.

Protocolo

NOTA: Abra el programa R-studio y cargue el archivo R "DEGs.R", el archivo se puede adquirir desde Archivos suplementarios / Scripts.

1. Descarga y preprocesamiento de datos

1. Descargue los datos de recuento de secuenciación de alto rendimiento (HTSeq) de colangiocarcinoma (CHOL) del Atlas del Genoma del Cáncer (TCGA). Este paso se puede lograr fácilmente con el siguiente código R.
   1. Haga clic en Ejecutar para instalar paquetes de R.
   2. Haga clic en Ejecutar para cargar paquetes de R.
      if(!requireNamespace("BiocManager", quietly=TRUE))
      + install.packages("BiocManager")
      BiocManager::install(c("TCGAbiolinks", "SummarizedExperiment"))
   3. Establezca el directorio de trabajo.
      biblioteca (TCGAbiolinks)
      library(ResumenExperimento)
      setwd("C:/Usuarios/LIUSHIYI/Escritorio")
   4. Elija el tipo de cáncer.
      cáncer <- "TCGA-CHOL"
   5. Ejecute el código R desde el archivo "GDCquery.R" para descargar los datos. El archivo "GDCquery.R" se puede adquirir a partir de archivos complementarios/scripts:
      source("Archivos suplementarios/Scripts/GDCquery.R")
      cabeza(cnt)
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      ##ENSG00000000003 4262
      ##ENSG00000000005 1
      ##ENSG00000000419 1254
      ##ENSG00000000457 699
      ##ENSG00000000460 239
      N.º #ENSG00000000938 334
      NOTA: Después de la ejecución, los datos de recuento de CHOLHTSeq se descargarán y se denominarán "cnt", donde las filas representan los identes de genes de conjunto y las columnas representan los identificados de muestra. Tenga en cuenta los números en las posiciones 14-15 en las identificaciones de muestra; números que van de 01 a 09 indican tumores y que van de 10 a 19 indican tejidos normales.
2. Convierta los ID de genes de conjunto en símbolos genéticos.
   1. Importe el archivo de anotación a R de acuerdo con su ruta de almacenamiento. El archivo de anotación (gencode.v22.annotation.gtf) se puede adquirir desde archivos suplementarios.
      gtf_v22 <- rtracklayer::import('Archivos suplementarios/gencode.v22.annotation.gtf')
   2. Ejecute el código R desde el "gtf_v22. R", que se puede adquirir a partir de archivos complementarios/scripts:
      source("Archivos suplementarios/Scripts/gtf_v22. R")
   3. Aplique la función "ann" para convertir los ID de genes del conjunto en símbolos de genes.
      cnt=ann(cnt,gtf_v22)
3. Filtrado de genes de baja expresación
   1. Haga clic en Ejecutar para instalar el paquete de R "edgeR".
      BiocManager::install("edgeR")
   2. Haga clic en Ejecutar para cargar el paquete de R "edgeR".
      biblioteca(edgeR)
   3. Ejecute el siguiente código R para mantener genes con valores de recuento por millón (CPM) superiores a uno de al menos dos muestras.
      keep <- rowSums(cpm(cnt)>1)>=2
      cnt <- as.matrix(cnt[keep,])
      NOTA: El valor de recuentos por millón (CPM) se utiliza en lugar del recuento de lectura para eliminar la desviación causada por diferentes profundidades de secuenciación.

2. Análisis de expresión diferencial a través de "limma"

1. Haga clic en Ejecutar para instalar el paquete R "limma".
   BiocManager::install("limma")
2. Haga clic en Ejecutar para cargar los paquetes de R "limma", "edgeR".
   biblioteca(limma)
   biblioteca(edgeR)
3. Ejecute el siguiente código de R para crear la matriz de diseño.
   grupo <- substring(colnames(cnt),14,15) # Extract group information
   grupo [grupo %in% "01"] <- "Cancer" # set '01' as tumor tissue
   grupo [grupo %en% "11"] <- "Normal" # set '11' as normal tissue
   grupo <- factor (group, levels = c("Normal","Cancer"))
   1. Cree la matriz de diseño.
      design <- model.matrix (~group)
      rownames(design) <- colnames(cnt)
   2. Cree el objeto DGEList.
      dge <- DGEList(counts = cnt, group = group)
   3. Normalizar los datos.
      dge <- calcNormFactors(dge, method = "TMM")
   4. Ejecute el siguiente código de R para realizar el análisis de expresión diferencial basado en el método limma-trend.
      dge
      ##An objeto de la clase "DGEList"
      ##$counts
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      ##TSPAN6 4262
      ##DPM1 1254
      ##SCYL3 699
      ##C1orf112 239
      N.º #FGR 334
   5. Calcule el valor de CPM.
      logdge <- cpm(dge, log=TRUE, prior.count=3)
   6. Haga clic en Ejecutar para ajustar un modelo lineal para predecir los datos o inferir la relación entre variables.
      ajuste <- lmFit (logdge, diseño)
   7. Calcule el valor T, el valor F y las probabilidades logarítmos basadas en bayesiano.
      fit <- eBayes(fit, trend=TRUE)
   8. Extraiga la tabla de resultados.
      res_limma<- as.data.frame(topTable(fit,n=Inf))
      
      cabeza(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 del análisis de expresión diferencial se guarda en "res_limma", que incluye el id del gen, el valor de cambio log2 fold (logFC), el nivel medio de expresión log2 del gen en el experimento (AveExpr), la estadística t modificada (t), el valor relavent p (P.Value), la tasa de descubrimiento falso (FDR) corregido p valor (adj. P.Val) y las probabilidades logarítmosa de genes expresados diferencialmente (B)
      NOTA: La función "calcNormFactors()" del "edgeR" se utilizó para normalizar los datos para eliminar la influencia causada por la preparación de muestras o la construcción y secuenciación de bibliotecas. En la construcción de la matriz de diseño, es necesario hacer coincidir el diseño experimental (por ejemplo, tipo de tejido: tejidos normales o tumorales) con las muestras de las direcciones de la matriz. limma-trend es adecuado para datos cuya profundidad de secuenciación es la misma, mientras que limma-voom es adecuado: (i) cuando el tamaño de la biblioteca de muestras es diferente; ii) datos no normalizados por TMM; (iii) hay mucho "ruido" en los datos. Un logFC positivo significa que el gen está regulado al pie en el experimento, mientras que el número negativo significa que el gen está regulado a la baja.
   9. Identificar los DEG.
      res_limma$sig <- as.factor(
      ifelse(res_limma$adj. P.Val < 0.05 & abs(res_limma$logFC) > 2,
      ifelse(res_limma$logFC > 2 ,'up','down'),'not')) # El valor de adj.p < 0.05 y el |log2FC| > = 2 son umbrales para identificar los DEG
      resumen(res_limma$sig)
      ##down no arriba
      ##1880 ​17341 1443
   10. Genere la tabla de resultados en un archivo.
       write.csv(res_limma, archivo = 'result_limma.csv')
   11. Haga clic en Ejecutar para instalar el paquete de R "ggplot2".
       install.packages("ggplot2")
   12. Haga clic en Ejecutar para cargar el paquete R "ggplot2".
       biblioteca(ggplot2)
   13. Ejecute el código R desde el "volcán". R" para crear la trama del volcán. El archivo "volcán. R" se puede adquirir a partir de archivos complementarios.
       source("Archivos suplementarios/Scripts/volcano. R")
       volcano(res_limma,"logFC","adj. P.Val",2,0.05)
       NOTA: Los genes se pueden mapear a diferentes posiciones de acuerdo con sus valores log2FC y adj-p, los DEG regulados hacia arriba están coloreados en rojo y los DEG regulados hacia abajo están coloreados en verde.
   14. Haga clic en Exportar para guardar el gráfico del volcán.
       NOTA: Los diagramas de volcanes se pueden generar y descargar en diferentes formatos (por ejemplo, pdf, TIFF, PNG, formato JPEG). Los genes se pueden mapear a diferentes posiciones de acuerdo con sus valores log2FC y adj p, los DEG regulados al pie (log2FC > 2, adj p < 0.05) están coloreados en rojo, y los DEG regulados a la baja (log2FC < -2, adj p < 0.05) están coloreados en verde, los no DEG están coloreados en gris.

3. Análisis de expresión diferencial a través de "edgeR"

1. Haga clic en Ejecutar para cargar el paquete de R "edgeR".
   biblioteca(edgeR)
2. Ejecute el siguiente código de R para crear una matriz de diseño.
   grupo <-subcadena(colnames(cnt),14,15)
   grupo [grupo %en% "01"] <- "Cáncer"
   grupo [grupo %en% "11"] <- "Normal"
   grupo=factor(grupo, niveles = c("Normal","Cáncer"))
   design <-model.matrix(~group)
   rownames(diseño) = colnames(cnt)
3. Haga clic en Ejecutar para crear el objeto DGEList.
   dge <- DGEList(counts=cnt)
4. Normalizar los datos.
   dge <- calcNormFactors(dge, method = "TMM")
5. Haga clic en Ejecutar para estimar la dispersión de los valores de expresión génica.
   dge <- estimateDisp(dge, diseño, robusto = T)
6. Haga clic en Ejecutar para ajustar el modelo y contar los datos.
   ajuste <- glmQLFit(dge, diseño)
7. Realizar una prueba estadística.
   ajuste <- glmQLFTest(fit)
8. Extraiga la tabla de resultados. El resultado se guarda en "res_edgeR", que incluye el valor de cambio de pliegue de registro, CPM de registro, F, valor de p y valor de p corregido de FDR.
   res_edgeR=as.data.frame(topTags(fit, n=Inf))
   cabeza(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 se guarda en "res_edgeR", que incluye el valor de cambio de pliegue de registro (logFC), CPM de registro, F, valor de p y valor de p corregido de FDR
9. Identificar los DEG.
   res_edgeR$sig = as.factor(
   ifelse(res_edgeR$FDR < 0.05 & abs(res_edgeR$logFC) > 2,
   ifelse(res_edgeR$logFC > 2 ,'up','down'),'not'))
   resumen(res_edgeR$sig)
   ##down no arriba
   ##1578 15965 3121
10. Genere la tabla de resultados en un archivo.
    write.csv(res_edgeR, archivo = 'res_edgeR.csv')
11. Crea la trama del volcán.
    volcano(res_edgeR,"logFC","FDR",2,0.05)
12. Haga clic en Exportar para guardar el gráfico del volcán.

4. Análisis de expresión diferencial a través de "DESeq2"

1. Haga clic en Ejecutar para instalar los paquetes de R "DESeq2".
   BiocManager::install("DESeq2")
2. Haga clic en Ejecutar para cargar los paquetes de R "DESeq2".
   biblioteca(DESeq2)
3. Ejecute el siguiente código de R para determinar el factor de agrupación.
   grupo <-subcadena(colnames(cnt),14,15)
   grupo [grupo %en% "01"] <- "Cáncer"
   grupo [grupo %en% "11"] <- "Normal"
   grupo=factor(grupo, niveles = c("Normal","Cáncer"))
4. Cree el objeto DESeqDataSet.
   dds <-DESeqDataSetFromMatrix (cnt, DataFrame(group), design = ~group)
   Dds
   ##class: DESeqDataSet
   ##dim: 20664 45
   ##metadata(1): versión
   ##assays(1): cuenta
   ##rownames(20664): TSPAN6 DPM1 ... RP11-274B21.13 LINC01144
   ##rowData nombres(0):
   ##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ...
   ##colData nombres(1): grupo
5. Realizar el análisis.
   dds <- DESeq(dds)
6. Genere la tabla de resultados.
   res_DESeq2 <- data.frame(results(dds))
   
   cabeza(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: El resultado se guarda en "res_DESeq2", que incluye la media del recuento de lectura normalizado (baseMean), el valor de cambio de pliegue de registro (log2FoldChange), el error estándar de cambio de pliegue de registro (lfcSE), la estadística de Wald (estadística), el valor p original (pvalue) y el valor p corregido (padj)
7. Identificar DEGs.
   res_DESeq2$sig = as.factor(
   ifelse(res_DESeq2$padj < 0.05 & abs(res_DESeq2$log2FoldChange) > 2,
   ifelse(res_DESeq2$log2FoldChange > 2 ,'up','down'),'not'))
   resumen(res_DESeq2$sig)
   ##down no arriba
   ##1616 16110 2938
8. Genere la tabla de resultados en un archivo.
   write.csv(res_DESeq2, archivo = 'res_DESeq2.csv')
9. Crea la trama del volcán.
   volcano(res_DESeq2,"log2FoldChange","padj",2,0.05)
10. Haga clic en Exportar para guardar el gráfico del volcán.

5. Diagrama de Venn

1. Haga clic en Ejecutar para instalar el paquete R "VennDiagram".
   install.packages("VennDiagram")
2. Haga clic en Ejecutar para cargar el paquete R "VennDiagram".
   biblioteca (VennDiagram)
3. Haga un diagrama de Venn de los DEG regulados.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=rownames(res_
   limma[res_limma$sig=="arriba",]),
   edgeR=rownames(res_edgeR[res_edgeR$sig=="up",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig==
   "arriba",])),
   NULL,altura = 3,ancho = 3,unidades = "en",
   col="negro",lwd=0.3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alfa=c(0.5, 0.5, 0.5),main = "DEGs regulados al up"))
4. Haga clic en Exportar para guardar el diagrama de Venn.
5. Haga un diagrama de Venn de DEG regulados hacia abajo.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=rownames(res_
   limma[res_limma$sig=="abajo",]),
   edgeR=rownames(res_edgeR[res_edgeR$sig==
   "abajo",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig=="abajo",])),
   NULL,altura = 3,ancho = 3,unidades = "en",
   col="negro",lwd=0.3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alfa=c(0.5, 0.5, 0.5),main = "DEGs regulados a la baja"))
6. Haga clic en Exportar para guardar el diagrama de Venn.

Resultados

Existen diversos enfoques para visualizar el resultado del análisis de expresión diferencial, entre los que se utilizan particularmente la gráfica del volcán y el diagrama de Venn. limma identificó 3323 DEGs entre el CHOL y los tejidos normales con el |logFC|≥2 y adj. P.Val <0,05 como umbrales, entre los cuales 1880 fueron regulados a la baja en los tejidos CHOL y 1443 fueron regulados al pie(Figura 1a). Mientras tanto, edgeR identificó los 1578 DEG regulados a la baja y 3121 DEG reg...

Discusión

Las abundantes transcripciones aberrantes en los cánceres se pueden identificar fácilmente mediante el análisis diferencial de ARN-seq⁵. Sin embargo, la aplicación del análisis de expresión diferencial RNA-seq a menudo está restringida, ya que requiere ciertas habilidades con el lenguaje R y la capacidad de elegir métodos apropiados. Para abordar este problema, proporcionamos una introducción detallada a los tres métodos más conocidos (limma, EdgeR y DESeq2) y tutoriales para aplicar el...

Materiales

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