Tre metodi di analisi dell'espressione differenziale per il sequenziamento dell'RNA: limma, EdgeR, DESeq2

Riepilogo

È stato fornito un protocollo dettagliato di metodi di analisi dell'espressione differenziale per il sequenziamento dell'RNA: limma, EdgeR, DESeq2.

Abstract

Il sequenziamento dell'RNA (RNA-seq) è una delle tecnologie più utilizzate nella trascrittomica in quanto può rivelare la relazione tra l'alterazione genetica e processi biologici complessi e ha un grande valore nella diagnostica, nella prognostica e nella terapia dei tumori. L'analisi differenziale dei dati RNA-seq è fondamentale per identificare trascrizioni aberranti e limma, EdgeR e DESeq2 sono strumenti efficienti per l'analisi differenziale. Tuttavia, l'analisi differenziale RNA-seq richiede determinate abilità con il linguaggio R e la capacità di scegliere un metodo appropriato, che manca nel curriculum di educazione medica.

Qui forniamo il protocollo dettagliato per identificare i geni differenzialmente espressi (DEG) tra il colangiocarcinoma (CHOL) e i tessuti normali attraverso limma, DESeq2 e EdgeR, rispettivamente, e i risultati sono mostrati in grafici vulcanici e diagrammi di Venn. I tre protocolli di limma, DESeq2 e EdgeR sono simili ma hanno passaggi diversi tra i processi di analisi. Ad esempio, un modello lineare viene utilizzato per le statistiche in limma, mentre la distribuzione binomiale negativa viene utilizzata in edgeR e DESeq2. Inoltre, i dati normalizzati del conteggio RNA-seq sono necessari per EdgeR e limma, ma non sono necessari per DESeq2.

Qui forniamo un protocollo dettagliato per tre metodi di analisi differenziale: limma, EdgeR e DESeq2. I risultati dei tre metodi sono in parte sovrapposti. Tutti e tre i metodi hanno i loro vantaggi e la scelta del metodo dipende solo dai dati.

Introduzione

Il sequenziamento dell'RNA (RNA-seq) è una delle tecnologie più utilizzate nella trascrittomica con molti vantaggi (ad esempio, un'elevata riproducibilità dei dati) e ha notevolmente aumentato la nostra comprensione delle funzioni e delle dinamiche di processi biologici complessi^1,2. L'identificazione di trascritti ablati in diversi contesti biologici, noti anche come geni differenzialmente espressi (DEG), è un passo chiave nell'analisi RNA-seq. RNA-seq consente di ottenere una profonda comprensione dei meccanismi molecolari e delle funzioni biologiche correlate alla patogenesi. Pertanto, l'analisi differenziale è stata considerata preziosa per la diagnostica, la prognostica e la terapia dei tumori^3,4,5. Attualmente, sono stati sviluppati più pacchetti R/Bioconductor open source per l'analisi dell'espressione differenziale RNA-seq, in particolare limma, DESeq2 e EdgeR^1,6,7. Tuttavia, l'analisi differenziale richiede alcune abilità con il linguaggio R e la capacità di scegliere il metodo appropriato, che manca nel curriculum di educazione medica.

In questo protocollo, basato sui dati di conteggio RNA-seq del colangiocarcinoma (CHOL) estratti da The Cancer Genome Atlas (TCGA), tre dei metodi più noti (limma^8,EdgeR⁹ e DESeq2¹⁰⁾sono stati effettuati, rispettivamente, dal programma R¹¹ per identificare i DEG tra CHOL e tessuti normali. I tre protocolli di limma, EdgeR e DESeq2 sono simili ma hanno passaggi diversi tra i processi di analisi. Ad esempio, i dati normalizzati del conteggio RNA-seq sono necessari per EdgeR e limma^8,9, mentre DESeq2 utilizza le proprie discrepanze di libreria per correggere i dati invece della normalizzazione¹⁰. Inoltre, edgeR è specificamente adatto per i dati RNA-seq, mentre il limma viene utilizzato per microarray e RNA-seq. Un modello lineare è adottato da limma per valutare i DEG¹², mentre le statistiche in edgeR si basano sulle distribuzioni binomiali negative, tra cui stima empirica di Bayes, test esatti, modelli lineari generalizzati e test di quasi-verosimiglianza⁹.

In sintesi, forniamo i protocolli dettagliati dell'analisi dell'espressione differenziale RNA-seq utilizzando rispettivamente limma, DESeq2 e EdgeR. Facendo riferimento a questo articolo, gli utenti possono facilmente eseguire l'analisi differenziale RNA-seq e scegliere i metodi di analisi differenziale appropriati per i loro dati.

Protocollo

NOTA: Aprire il programma R-studio e caricare il file R "DEGs.R", il file può essere acquisito da file / script supplementari.

1. Download e pre-elaborazione dei dati

1. Scarica i dati di conteggio del sequenziamento ad alto rendimento (HTSeq) del colangiocarcinoma (CHOL) da The Cancer Genome Atlas (TCGA). Questo passaggio può essere facilmente ottenuto con il seguente codice R.
   1. Fare clic su Esegui per installare i pacchetti R.
   2. Fare clic su Esegui per caricare i pacchetti R.
      if(!requireNamespace("BiocManager", quietly=TRUE))
      + install.packages("BiocManager")
      BiocManager::install(c("TCGAbiolinks", "SummarizedExperiment"))
   3. Impostare la directory di lavoro.
      libreria (TCGAbiolinks)
      libreria(SummarizedExperiment)
      setwd("C:/Utenti/LIUSHIYI/Desktop")
   4. Scegli il tipo di cancro.
      cancro <- "TCGA-CHOL"
   5. Eseguire il codice R dal file "GDCquery.R" per scaricare i dati. Il file "GDCquery.R" può essere acquisito da file/script supplementari:
      source("File/script/GDCquery.R supplementari")
      testa(cnt)
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      N. #ENSG00000000003 4262
      ##ENSG00000000005 1
      N. #ENSG00000000419 1254
      N. #ENSG00000000457 699
      N. #ENSG00000000460 239
      N. #ENSG00000000938 334
      NOTA: dopo l'esecuzione, i dati di conteggio CHOLHTSeq verranno scaricati e denominati "cnt", dove le righe rappresentano gli ID dei geni dell'insieme e le colonne rappresentano gli ID campione. Si prega di notare i numeri nelle posizioni 14-15 negli ID di esempio; numeri che vanno da 01 a 09 indicano tumori e vanno da 10 a 19 indicano tessuti normali.
2. Convertire gli ID dei geni dell'insieme in simboli genetici.
   1. Importare il file di annotazione in R in base al relativo percorso di archiviazione. Il file di annotazione (gencode.v22.annotation.gtf) può essere acquisito da file supplementari.
      gtf_v22 <- rtracklayer::import('File supplementari/gencode.v22.annotation.gtf')
   2. Eseguire il codice R dal "gtf_v22. R", che può essere acquisito da file/script supplementari:
      source("File/Script/gtf_v22 supplementari. R")
   3. Applicare la funzione "ann" per convertire gli ID del gene dell'insieme in simboli genici.
      cnt=ann(cnt;gtf_v22)
3. Filtraggio di geni a bassa esprime
   1. Fare clic su Esegui per installare il pacchetto R "edgeR".
      BiocManager::install("edgeR")
   2. Fare clic su Esegui per caricare il pacchetto R "edgeR".
      libreria(edgeR)
   3. Eseguire il codice R seguente per mantenere i geni con valori di conteggio per milione (CPM) maggiori di uno in almeno due campioni.
      mantenere <- rowSums(cpm(cnt)>1)>=2
      cnt <- as.matrix(cnt[keep,])
      NOTA: il valore dei conteggi per milione (CPM) viene utilizzato al posto del conteggio delle lette per eliminare la deviazione causata da diverse profondità di sequenziamento.

2. Analisi dell'espressione differenziale attraverso "limma"

1. Fare clic su Esegui per installare il pacchetto R "limma".
   BiocManager::install("limma")
2. Fare clic su Esegui per caricare i pacchetti R "limma", "edgeR".
   biblioteca(limma)
   libreria(edgeR)
3. Eseguire il codice R seguente per creare la matrice di progettazione.
   gruppo <- substring(colnames(cnt),14,15) # Extract group information
   gruppo [gruppo %in% "01"] <- "Cancer" # set '01' as tumor tissue
   gruppo [gruppo %in% "11"] <- "Normal" # set '11' as normal tissue
   gruppo <- factor (group, levels = c("Normal","Cancer"))
   1. Creare la matrice di progettazione.
      design <- model.matrix (~gruppo)
      rownames(design) <- colnames(cnt)
   2. Creare l'oggetto DGEList.
      dge <- DGEList(counts = cnt, group = group)
   3. Normalizzare i dati.
      dge <- calcNormFactors(dge, method = "TMM")
   4. Eseguire il codice R riportato di seguito per eseguire l'analisi dell'espressione differenziale basata sul metodo limma-trend.
      dge
      ##An oggetto di classe "DGEList"
      ##$counts
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      N. #TSPAN6 4262
      N. #DPM1 1254
      N. #SCYL3 699
      N. #C1orf112 239
      N. #FGR 334
   5. Calcolare il valore CPM.
      logdge <- cpm(dge, log=TRUE, prior.count=3)
   6. Fare clic su Esegui per adattare un modello lineare per prevedere i dati o dedurre la relazione tra le variabili.
      fit <- lmFit (logdge, design)
   7. Calcola il valore T, il valore F e le quote di log in base al bayesiano.
      fit <- eBayes(fit, trend=TRUE)
   8. Estrarre la tabella dei risultati.
      res_limma<- as.data.frame(topTable(fit,n=Inf))
      
      testa(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 risultato dell'analisi dell'espressione differenziale viene salvato in "res_limma", che include l'id del gene, il valore di modifica della piega log2 log2 (logFC), il livello medio di espressione log2 del gene nell'esperimento (AveExpr), la statistica t modificata (t), il valore p relavent (P.Value), il tasso di falsa scoperta (FDR) il valore p corretto (adj. P.Val) e le log-odds dei geni differenzialmente espressi (B)
      NOTA: la funzione "calcNormFactors()" di "edgeR" è stata utilizzata per normalizzare i dati per eliminare l'influenza causata dalla preparazione del campione o dalla costruzione e sequenziamento della libreria. Nella costruzione della matrice di progettazione, è necessario abbinare il design sperimentale (ad esempio, tipo di tessuto: tessuti normali o tumorali) agli ID campione della matrice. limma-trend è adatto a dati la cui profondità di sequenziamento è la stessa, mentre limma-voom è adatto: (i) quando la dimensione della libreria di campioni è diversa; (ii) dati non normalizzati da TMM; (iii) c'è molto "rumore" nei dati. Un logFC positivo significa che il gene è up-regolato nell'esperimento, mentre il numero negativo significa che il gene è down-regolato.
   9. Identificare i 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')) # Il valore adj.p < 0.05 e il valore |log2FC| >= 2 sono soglie per identificare i DEG
      sommario(res_limma$sig)
      ##down non su
      ##1880 ​17341 1443
   10. Eseguire l'output della tabella dei risultati in un file.
       write.csv(res_limma, file = 'result_limma.csv')
   11. Fare clic su Esegui per installare il pacchetto R "ggplot2".
       install.packages("ggplot2")
   12. Fare clic su Esegui per caricare il pacchetto R "ggplot2".
       biblioteca(ggplot2)
   13. Esegui il codice R dal "vulcano. R" per creare la trama del vulcano. Il file "vulcano. R" può essere acquisito da file supplementari.
       source("File supplementari/Script/vulcano. R")
       volcano(res_limma,"logFC","adj. P.Val",2.0.05)
       NOTA: i geni possono essere mappati in posizioni diverse in base ai loro valori log2FC e adj-p, i DEG regolati verso l'alto sono colorati in rosso e i DEG regolati verso il basso sono colorati in verde.
   14. Fare clic su Esporta per salvare la trama del vulcano.
       NOTA: Le trame del vulcano possono essere generate e scaricate in diversi formati (ad esempio, pdf, TIFF, PNG, JPEG). I geni possono essere mappati in posizioni diverse in base ai loro valori log2FC e adj p, i DEG up-regolati (log2FC > 2, adj p < 0,05) sono colorati in rosso e i DEG down-regolati (log2FC < -2, adj p < 0,05) sono colorati in verde, i non-DEG sono colorati in grigio.

3. Analisi dell'espressione differenziale tramite "edgeR"

1. Fare clic su Esegui per caricare il pacchetto R "edgeR".
   libreria(edgeR)
2. Eseguire il codice R riportato di seguito per creare una matrice di progettazione.
   gruppo <-sottostringa(colnames(cnt),14,15)
   gruppo [gruppo %in% "01"] <- "Cancro"
   gruppo [gruppo %in% "11"] <- "Normale"
   gruppo=fattore(gruppo, livelli = c("Normale","Cancro"))
   design <-model.matrix(~group)
   rownames(design) = colnames(cnt)
3. Fare clic su Esegui per creare l'oggetto DGEList.
   dge <- DGEList(counts=cnt)
4. Normalizzare i dati.
   dge <- calcNormFactors(dge, method = "TMM")
5. Fare clic su Esegui per stimare la dispersione dei valori di espressione genica.
   dge <- estimateDisp(dge, design, robust = T)
6. Fare clic su Esegui per adattare il modello al conteggio dei dati.
   fit <- glmQLFit (dge, design)
7. Condurre un test statistico.
   fit <- glmQLFTest(fit)
8. Estrarre la tabella dei risultati. Il risultato viene salvato in "res_edgeR", che include il valore di modifica della piega del registro, il valore CPM del registro, F, il valore p e il valore p corretto FDR.
   res_edgeR=as.data.frame(topTags(fit, n=Inf))
   testa(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 risultato viene salvato in "res_edgeR", che include il valore di modifica della piega del registro (logFC), il valore CPM del registro, il valore F, p e il valore p corretto FDR
9. Identificare i 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'))
   sommario(res_edgeR$sig)
   ##down non su
   ##1578 15965 3121
10. Eseguire l'output della tabella dei risultati in un file.
    write.csv(res_edgeR, file = 'res_edgeR.csv')
11. Crea la trama del vulcano.
    vulcano(res_edgeR,"logFC","FDR",2,0.05)
12. Fare clic su Esporta per salvare la trama del vulcano.

4. Analisi dell'espressione differenziale tramite "DESeq2"

1. Fare clic su Esegui per installare i pacchetti R "DESeq2".
   BiocManager::install("DESeq2")
2. Fare clic su Esegui per caricare i pacchetti R "DESeq2".
   biblioteca(DESeq2)
3. Eseguire il codice R riportato di seguito per determinare il fattore di raggruppamento.
   gruppo <-sottostringa(colnames(cnt),14,15)
   gruppo [gruppo %in% "01"] <- "Cancro"
   gruppo [gruppo %in% "11"] <- "Normale"
   gruppo=fattore(gruppo, livelli = c("Normale","Cancro"))
4. Creare l'oggetto DESeqDataSet.
   dds <-DESeqDataSetFromMatrix (cnt, DataFrame(group), design = ~group)
   Dds
   ##class: DESeqDataSet
   ##dim: 20664 45
   ##metadata(1): versione
   ##assays(1): conta
   # #rownames(20664): TSPAN6 DPM1 ... RP11-274B21,13 LINC01144
   ##rowData nomi(0):
   # #colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ...
   ##colData nomi(1): gruppo
5. Eseguire l'analisi.
   dds <- DESeq(dds)
6. Generare la tabella dei risultati.
   res_DESeq2 <- data.frame(results(dds))
   
   testa(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: il risultato viene salvato in "res_DESeq2", che include la media del conteggio di lettura normalizzato (baseMean), il valore di modifica della piega del registro (log2FoldChange), l'errore standard di modifica della piega del registro (lfcSE), la statistica di Wald (stat), il valore p originale (pvalue) e il valore p corretto (padj)
7. Identificare i DEG.
   res_DESeq2$sig = as.factor(
   ifelse(res_DESeq2$padj < 0.05 & abs(res_DESeq2$log2FoldChange) > 2,
   ifelse(res_DESeq2$log2FoldChange > 2 ,'up','down'),'not'))
   sommario(res_DESeq2$sig)
   ##down non su
   ##1616 16110 2938
8. Eseguire l'output della tabella dei risultati in un file.
   write.csv(res_DESeq2, file = 'res_DESeq2.csv')
9. Crea la trama del vulcano.
   volcano(res_DESeq2,"log2FoldChange","padj",2,0.05)
10. Fare clic su Esporta per salvare la trama del vulcano.

5. Diagramma di Venn

1. Fare clic su Esegui per installare il pacchetto R "VennDiagram".
   install.packages("VennDiagram")
2. Fare clic su Esegui per caricare il pacchetto R "VennDiagram".
   biblioteca (VennDiagram)
3. Crea un diagramma di Venn dei DEG regolamentati.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=rownames(res_
   limma[res_limma$sig=="up",]),
   edgeR=rownames(res_edgeR[res_edgeR$sig=="up",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig==
   "su",])),
   NULL,altezza = 3,larghezza = 3,unità = "in",
   col="nero",lwd=0.3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alpha=c(0.5, 0.5, 0.5),main = "DEG up-regulated"))
4. Fare clic su Esporta per salvare il diagramma di Venn.
5. Crea un diagramma di Venn dei DEG down regolamentati.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=rownames(res_
   limma[res_limma$sig=="giù",]),
   edgeR=rownames(res_edgeR[res_edgeR$sig==
   "giù",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig=="down",])),
   NULL,altezza = 3,larghezza = 3,unità = "in",
   col="nero",lwd=0.3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alpha=c(0.5, 0.5, 0.5),main = "DEG down-regulated"))
6. Fare clic su Esporta per salvare il diagramma di Venn.

Risultati

Esistono vari approcci per visualizzare il risultato dell'analisi dell'espressione differenziale, tra cui il diagramma del vulcano e il diagramma di Venn sono particolarmente utilizzati. limma ha identificato 3323 DEG tra il CHOL e i tessuti normali con il |logFC|≥2 e adj. P.Val <0,05 come soglie, tra cui 1880 erano down-regolati nei tessuti CHOL e 1443 erano up-regolati (Figura 1a). Nel frattempo, edgeR ha identificato i 1578 DEG down-regulated e i 3121 DEG up-regulated (

Discussione

Abbondanti trascrizioni ablate nei tumori possono essere facilmente identificate dall'analisi differenziale RNA-seq⁵. Tuttavia, l'applicazione dell'analisi dell'espressione differenziale RNA-seq è spesso limitata in quanto richiede determinate abilità con linguaggio R e la capacità di scegliere metodi appropriati. Per affrontare questo problema, forniamo un'introduzione dettagliata ai tre metodi più noti (limma, EdgeR e DESeq2) e tutorial per l'applicazione dell'analisi dell'espressione differ...

Materiali

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