Summary
Ein detailliertes Protokoll der methoden der differentiellen Expressionsanalyse für die RNA-Sequenzierung wurde bereitgestellt: limma, EdgeR, DESeq2.
Abstract
Die RNA-Sequenzierung (RNA-seq) ist eine der am weitesten verbreiteten Technologien in der Transkriptomik, da sie den Zusammenhang zwischen der genetischen Veränderung und komplexen biologischen Prozessen aufdecken kann und einen großen Wert in der Diagnostik, Prognose und Therapeutik von Tumoren hat. Die Differentialanalyse von RNA-seq-Daten ist entscheidend, um aberrante Transkriptionen zu identifizieren, und limma, EdgeR und DESeq2 sind effiziente Werkzeuge für die Differentialanalyse. Die RNA-seq-Differentialanalyse erfordert jedoch bestimmte Fähigkeiten mit der Sprache R und die Fähigkeit, eine geeignete Methode zu wählen, die im Lehrplan der medizinischen Ausbildung fehlt.
Hierin stellen wir das detaillierte Protokoll zur Identifizierung differentiell exprimierter Gene (DEGs) zwischen Cholangiokarzinom (CHOL) und normalem Gewebe durch Limma, DESeq2 bzw. EdgeR zur Verfügung, und die Ergebnisse werden in Vulkandiagrammen und Venn-Diagrammen gezeigt. Die drei Protokolle limma, DESeq2 und EdgeR sind ähnlich, haben aber unterschiedliche Schritte zwischen den Prozessen der Analyse. Beispielsweise wird ein lineares Modell für die Statistik in Limma verwendet, während die negative Binomialverteilung in edgeR und DESeq2 verwendet wird. Darüber hinaus sind die normalisierten RNA-Seq-Zähldaten für EdgeR und Limma notwendig, aber nicht für DESeq2.
Hier stellen wir ein detailliertes Protokoll für drei Differentialanalysemethoden zur Verfügung: limma, EdgeR und DESeq2. Die Ergebnisse der drei Methoden überschneiden sich teilweise. Alle drei Methoden haben ihre eigenen Vorteile, und die Wahl der Methode hängt nur von den Daten ab.
Introduction
Die RNA-Sequenzierung (RNA-seq) ist eine der am weitesten verbreiteten Technologien in der Transkriptomik mit vielen Vorteilen (z. B. hohe Datenreproduzierbarkeit) und hat unser Verständnis der Funktionen und Dynamik komplexer biologischer Prozesse dramatisch verbessert1,2. Die Identifizierung von Aberrat-Transkripten unter verschiedenen biologischen Kontexten, die auch als differentiell exprimierte Gene (DEGs) bezeichnet werden, ist ein wichtiger Schritt in der RNA-seq-Analyse. RNA-seq ermöglicht ein tiefes Verständnis der pathogenesebezogenen molekularen Mechanismen und biologischen Funktionen. Daher wurde die Differentialanalyse als wertvoll für die Diagnostik, Prognose und Therapie von Tumoren3,4,5angesehen. Derzeit wurden weitere Open-Source-R/Bioconductor-Pakete für die RNA-seq-Differentialexpressionsanalyse entwickelt, insbesondere limma, DESeq2 und EdgeR1,6,7. Die Differentialanalyse erfordert jedoch bestimmte Fähigkeiten mit der Sprache R und die Fähigkeit, die geeignete Methode zu wählen, die im Lehrplan der medizinischen Ausbildung fehlt.
In diesem Protokoll wurden basierend auf den Cholangiokarzinom (CHOL) RNA-seq-Zähldaten, die aus dem Cancer Genome Atlas (TCGA) extrahiert wurden, drei der bekanntesten Methoden (Limma8, EdgeR9 und DESeq210) vom R-Programm11 durchgeführt, um die DEGs zwischen CHOL und normalem Gewebe zu identifizieren. Die drei Protokolle limma, EdgeR und DESeq2 sind ähnlich, haben aber unterschiedliche Schritte zwischen den Prozessen der Analyse. Zum Beispiel sind die normalisierten RNA-seq-Zähldaten für EdgeR und Limma8,9notwendig, während DESeq2 seine eigenen Bibliotheksdiskrepanzen verwendet, um Daten anstelle von Normalisierung10zu korrigieren. Darüber hinaus eignet sich edgeR speziell für RNA-seq-Daten, während das limma für Microarrays und RNA-seq verwendet wird. Ein lineares Modell wird von limma zur Bewertung der DEGs12übernommen, während die Statistiken in edgeR auf den negativen Binomialverteilungen basieren, einschließlich empirischer Bayes-Schätzung, exakter Tests, verallgemeinerter linearer Modelle und Quasi-Wahrscheinlichkeitstests9.
Zusammenfassend stellen wir die detaillierten Protokolle der RNA-seq-Differentialexpressionsanalyse unter Verwendung von limma, DESeq2 bzw. EdgeR zur Verfügung. Durch bezugnahmend auf diesen Artikel können Benutzer die RNA-seq-Differentialanalyse einfach durchführen und die geeigneten Differentialanalysemethoden für ihre Daten auswählen.
Subscription Required. Please recommend JoVE to your librarian.
Protocol
HINWEIS: Öffnen Sie das R-studio-Programm und laden Sie die R-Datei "DEGs.R", die Datei kann aus Ergänzenden Dateien/Skripten erworben werden.
1. Herunterladen und Vorverarbeitung von Daten
- Laden Sie die Hochdurchsatz-Sequenzierungsdaten (HTSeq) des Cholangiokarzinoms (CHOL) aus dem Cancer Genome Atlas (TCGA) herunter. Dieser Schritt kann leicht durch den folgenden R-Code erreicht werden.
- Klicken Sie auf Ausführen, um R-Pakete zu installieren.
- Klicken Sie auf Ausführen, um R-Pakete zu laden.
if(!requireNamespace("BiocManager", quietly=TRUE))
+ install.packages("BiocManager")
BiocManager::install(c("TCGAbiolinks", "SummarizedExperiment")) - Legen Sie das Arbeitsverzeichnis fest.
Bibliothek (TCGAbiolinks)
Bibliothek(Zusammengefasstes Experiment)
setwd("C:/Benutzer/LIUSHIYI/Desktop") - Wählen Sie die Krebsart.
Krebs < - "TCGA-CHOL" - Führen Sie den R-Code aus der Datei "GDCquery.R" aus, um die Daten herunterzuladen. Die Datei "GDCquery.R" kann aus Ergänzenden Dateien/Skripten entnommen werden:
source("Ergänzende Dateien/Skripte/GDCquery.R")
Kopf(cnt)
Nr. #TCGA-3X-AAVA-01A-11R-A41I-07
Nr. #ENSG00000000003 4262
Nr. #ENSG00000000005 1
Nr. #ENSG00000000419 1254
Nr. #ENSG00000000457 699
Nr. #ENSG00000000460 239
Nr. #ENSG00000000938 334
HINWEIS: Nach der Ausführung werden die CHOLHTSeq-Zähldaten heruntergeladen und "cnt" genannt, wobei Zeilen Ensemble-Gen-IDs und Spalten Beispiel-IDs darstellen. Bitte beachten Sie die Nummern an den Positionen 14-15 in den Muster-IDs; Zahlen von 01 bis 09 weisen auf Tumore hin und von 10 bis 19 auf normale Gewebe.
- Konvertieren Sie Ensemble-Gen-IDs in Gensymbole.
- Importieren Sie die Anmerkungsdatei entsprechend ihrem Speicherpfad in R. Die Anmerkungsdatei (gencode.v22.annotation.gtf) kann aus Ergänzungsdateien abgebildet werden.
gtf_v22 <- rtracklayer::import('Supplementary files/gencode.v22.annotation.gtf') - Führen Sie den R-Code aus dem "gtf_v22. R"-Datei, die aus Ergänzenden Dateien/Skripten erworben werden kann:
source("Ergänzende Dateien/Skripte/gtf_v22. R") - Wenden Sie die Funktion "ann" an, um die Ensemble-Gen-IDs in Gensymbole umzuwandeln.
cnt=ann(cnt;gtf_v22)
- Importieren Sie die Anmerkungsdatei entsprechend ihrem Speicherpfad in R. Die Anmerkungsdatei (gencode.v22.annotation.gtf) kann aus Ergänzungsdateien abgebildet werden.
- Filtern von niedrig exprimierten Genen
- Klicken Sie auf Ausführen, um das R-Paket "edgeR" zu installieren.
BiocManager::install("edgeR") - Klicken Sie auf Ausführen, um das R-Paket "edgeR" zu laden.
Bibliothek(edgeR) - Führen Sie den folgenden R-Code aus, um Gene mit CPM-Werten (Counts per Million) größer als eins in mindestens zwei Stichproben zu halten.
keep <- rowSums(cpm(cnt)>1)>=2
cnt <- as.matrix(cnt[keep,])
HINWEIS: Der CPM-Wert (Counts per Million) wird anstelle der Leseanzahl verwendet, um die durch unterschiedliche Sequenzierungstiefen verursachte Abweichung zu eliminieren.
- Klicken Sie auf Ausführen, um das R-Paket "edgeR" zu installieren.
2. Differentielle Expressionsanalyse durch "Limma"
- Klicken Sie auf Ausführen, um das R-Paket "limma" zu installieren.
BiocManager::install("limma") - Klicken Sie auf Ausführen, um die R-Pakete "limma", "edgeR" zu laden.
Bibliothek(limma)
Bibliothek(edgeR) - Führen Sie den folgenden R-Code aus, um die Entwurfsmatrix zu erstellen.
Gruppe <- substring(colnames(cnt),14,15) # Extract group information
Gruppe [Gruppe %in% "01"] <- "Cancer" # set '01' as tumor tissue
Gruppe [Gruppe %in% "11"] <- "Normal" # set '11' as normal tissue
Gruppe <- factor (group, levels = c("Normal","Cancer"))- Erstellen Sie die Entwurfsmatrix.
Design <- model.matrix (~group)
Rownames(Design) <- colnames(cnt) - Erstellen Sie das DGEList-Objekt.
dge <- DGEList(counts = cnt, group = group) - Normalisieren Sie die Daten.
dge <- calcNormFactors(dge, method = "TMM") - Führen Sie den folgenden R-Code aus, um die auf der Limma-Trend-Methode basierende Differentielle Ausdrucksanalyse durchzuführen.
dge
##An Objekt der Klasse "DGEList"
##$counts
Nr. #TCGA-3X-AAVA-01A-11R-A41I-07
Nr. #TSPAN6 4262
Nr. #DPM1 1254
Nr. #SCYL3 699
Nr. #C1orf112 239
Nr. #FGR 334 - Berechnen Sie den CPM-Wert.
logdge <- cpm(dge, log=TRUE, prior.count=3) - Klicken Sie auf Ausführen, um ein lineares Modell anzupassen, um die Daten vorherzusagen oder die Beziehung zwischen Variablen abzuleiten.
fit <- lmFit (logdge, Design) - Berechnen Sie den T-Wert, den F-Wert und die Log-Quoten basierend auf Bayesian.
fit <- eBayes(fit, trend=TRUE) - Extrahieren Sie die Ergebnistabelle.
res_limma<- as.data.frame(topTable(fit,n=Inf))
Kopf(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 Ergebnis der differentiellen Expressionsanalyse wird in "res_limma" gespeichert, das die Gen-ID, den log2-Fachänderungswert (logFC), das durchschnittliche log2-Expressionsniveau des Gens im Experiment (AveExpr), die modifizierte t-Statistik (t), den relaventen p-Wert (P.-Wert), den falschen Entdeckungsrates-korrigierten p-Wert (FDR) (adj. P.Val) und die Log-Quoten differentiell exprimierter Gene (B)
HINWEIS: Die Funktion "calcNormFactors()" des "edgeR" wurde verwendet, um die Daten zu normalisieren, um den Einfluss zu eliminieren, der durch Probenvorbereitung oder Bibliotheksaufbau und -sequenzierung verursacht wird. Bei der Konstruktion der Designmatrix ist es notwendig, das experimentelle Design (z. B. Gewebetyp: Normales oder Tumorgewebe) mit den Proben-IDs der Matrix abzugleichen. limma-trend eignet sich für Daten, deren Sequenzierungstiefe gleich ist, während limma-voom geeignet ist: (i) wenn die Größe der Stichprobenbibliothek unterschiedlich ist; ii) Daten, die nicht durch TMM normalisiert wurden; (iii) es gibt viel "Rauschen" in den Daten. Ein positiver logFC bedeutet, dass das Gen im Experiment hochreguliert wird, während die negative Zahl bedeutet, dass das Gen herunterreguliert wird. - Identifizieren Sie die DEGs.
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')) # Der adj.p-Wert < 0,05 und der |log2FC| >= 2 sind Schwellenwerte zur Identifizierung der DEGs
Zusammenfassung(res_limma$sig)
##down nicht oben
##1880 17341 1443 - Geben Sie die Ergebnistabelle in eine Datei aus.
write.csv(res_limma, file = 'result_limma.csv') - Klicken Sie auf Ausführen, um das R-Paket "ggplot2" zu installieren.
install.packages("ggplot2") - Klicken Sie auf Ausführen, um das R-Paket "ggplot2" zu laden.
Bibliothek(ggplot2) - Führen Sie den R-Code vom Vulkan aus. R", um die Vulkanhandlung zu erstellen. Die Datei "volcano. R" kann aus Ergänzungsdateien erworben werden.
Quelle("Ergänzende Dateien/Skripte/Vulkan. R")
vulkan(res_limma,"logFC","adj. P.Val",2,0.05)
HINWEIS: Gene können entsprechend ihren log2FC- und adj-p-Werten auf verschiedene Positionen abgebildet werden, die nach oben regulierten DEGs sind rot gefärbt und die herunterregulierten DEGs sind grün gefärbt. - Klicken Sie auf Exportieren, um das Vulkandiagramm zu speichern.
HINWEIS: Die Vulkanplots können in verschiedenen Formaten (z.B. pdf, TIFF, PNG, JPEG Format) generiert und heruntergeladen werden. Gene können entsprechend ihren log2FC- und adj p-Werten auf verschiedene Positionen abgebildet werden, die hochregulierten DEGs (log2FC > 2, adj p < 0,05) sind rot eingefärbt und die herunterregulierten DEGs (log2FC < -2, adj p < 0,05) sind grün eingefärbt, Nicht-DEGs sind grau gefärbt.
- Erstellen Sie die Entwurfsmatrix.
3. Differentialausdrucksanalyse durch "edgeR"
- Klicken Sie auf Ausführen, um das R-Paket "edgeR" zu laden.
Bibliothek(edgeR) - Führen Sie den folgenden R-Code aus, um eine Entwurfsmatrix zu erstellen.
Gruppe <-Teilzeichenfolge(colnames(cnt),14,15)
Gruppe [Gruppe %in% "01"] <- "Krebs"
Gruppe [Gruppe %in% "11"] <- "Normal"
group=factor(group, levels = c("Normal","Krebs"))
Design <-model.matrix(~group)
Rownames(Design) = Colnames(cnt) - Klicken Sie auf Ausführen, um das DGEList-Objekt zu erstellen.
dge <- DGEList(counts=cnt) - Normalisieren Sie die Daten.
dge <- calcNormFactors(dge, method = "TMM") - Klicken Sie auf Ausführen, um die Streuung der Genexpressionswerte abzuschätzen.
dge <- estimateDisp(dge, design, robust = T) - Klicken Sie auf Ausführen, um das Modell zum Zählen von Daten anzupassen.
fit <- glmQLFit(dge, design) - Führen Sie einen statistischen Test durch.
fit <- glmQLFTest(fit) - Extrahieren Sie die Ergebnistabelle. Das Ergebnis wird in "res_edgeR" gespeichert, das den Log-Fold-Änderungswert, den Log-CPM,F-, p-Wert und den FDR-korrigierten p-Wert enthält.
res_edgeR=as.data.frame(topTags(fit, n=Inf))
Kopf(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 Ergebnis wird in "res_edgeR" gespeichert, das den Log Fold Change Value (logFC), log CPM, F, p-Wert und FDR-korrigierten p-Wert enthält - Identifizieren Sie die DEGs.
res_edgeR$sig = as.factor(
ifelse(res_edgeR$FDR < 0.05 & abs(res_edgeR$logFC) > 2,
ifelse(res_edgeR$logFC > 2 ,'up','down'),'not'))
Zusammenfassung(res_edgeR$sig)
##down nicht oben
##1578 15965 3121 - Geben Sie die Ergebnistabelle in eine Datei aus.
write.csv(res_edgeR, file = 'res_edgeR.csv') - Erstellen Sie das Vulkandiagramm.
Vulkan(res_edgeR,"logFC","FDR",2,0.05) - Klicken Sie auf Exportieren, um das Vulkandiagramm zu speichern.
4. Differentialexpressionsanalyse durch "DESeq2"
- Klicken Sie auf Ausführen, um die R-Pakete "DESeq2" zu installieren.
BiocManager::install("DESeq2") - Klicken Sie auf Ausführen, um R-Pakete "DESeq2" zu laden.
Bibliothek(DESeq2) - Führen Sie den folgenden R-Code aus, um den Gruppierungsfaktor zu bestimmen.
Gruppe <-Teilzeichenfolge(colnames(cnt),14,15)
Gruppe [Gruppe %in% "01"] <- "Krebs"
Gruppe [Gruppe %in% "11"] <- "Normal"
group=factor(group, levels = c("Normal","Krebs")) - Erstellen Sie das DESeqDataSet-Objekt.
dds <-DESeqDataSetFromMatrix (cnt, DataFrame(group), design = ~group)
Dds
##class: DESeqDataSet
Nr. #dim: 20664 45
##metadata(1): Version
##assays(1): zählt
##rownames(20664): TSPAN6 DPM1 ... RP11-274B21.13 LINC01144
##rowData Namen(0):
##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ...
##colData Namen(1): Gruppe - Führen Sie die Analyse durch.
dds <- DESeq(dds) - Generieren Sie die Ergebnistabelle.
res_DESeq2 <- data.frame(results(dds))
Kopf(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
HINWEIS: Das Ergebnis wird in "res_DESeq2" gespeichert, das den Mittelwert der normalisierten Leseanzahl (baseMean), den Log-Fold-Änderungswert (log2FoldChange), den Log-Fold-Change-Standardfehler (lfcSE), die Wald-Statistik (stat), den ursprünglichen p-Wert (pvalue) und den korrigierten p-Wert (padj) enthält. - Identifizieren Sie 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'))
Zusammenfassung(res_DESeq2$sig)
##down nicht oben
##1616 16110 2938 - Geben Sie die Ergebnistabelle in eine Datei aus.
write.csv(res_DESeq2, file = 'res_DESeq2.csv') - Erstellen Sie das Vulkandiagramm.
vulkan(res_DESeq2,"log2FoldChange","padj",2,0.05) - Klicken Sie auf Exportieren, um das Vulkandiagramm zu speichern.
5. Venn-Diagramm
- Klicken Sie auf Ausführen, um das R-Paket "VennDiagram" zu installieren.
install.packages("VennDiagram") - Klicken Sie auf Ausführen, um das R-Paket "VennDiagram" zu laden.
Bibliothek (VennDiagram) - Erstellen Sie ein Venn-Diagramm von regulierten DEGs.
grid.newpage()
grid.draw(venn.diagram(list(Limma=rownames(res_
limma[res_limma$sig=="up",]),
edgeR=rownames(res_edgeR[res_edgeR$sig=="up",]),
DESeq2=Zeilennamen(res_DESeq2[res_DESeq2$sig==
"up",])),
NULL,Höhe = 3,Breite = 3,Einheiten = "in",
col="black",lwd=0.3,fill=c("#FF6666","#FFFF00",
"#993366"),
alpha=c(0,5, 0,5, 0,5),main = "Hochregulierte DEGs")) - Klicken Sie auf Exportieren, um das Venn-Diagramm zu speichern.
- Erstellen Sie ein Venn-Diagramm von herunterregulierten DEGs.
grid.newpage()
grid.draw(venn.diagram(list(Limma=rownames(res_
limma[res_limma$sig=="down",]),
edgeR=Zeilennamen(res_edgeR[res_edgeR$sig==
"unten",]),
DESeq2=rownames(res_DESeq2[res_DESeq2$sig=="down",])),
NULL,Höhe = 3,Breite = 3,Einheiten = "in",
col="black",lwd=0.3,fill=c("#FF6666","#FFFF00",
"#993366"),
alpha=c(0,5, 0,5, 0,5),main = "Down-regulated DEGs")) - Klicken Sie auf Exportieren, um das Venn-Diagramm zu speichern.
Subscription Required. Please recommend JoVE to your librarian.
Representative Results
Es gibt verschiedene Ansätze, um das Ergebnis der Differentialexpressionsanalyse zu visualisieren, unter denen insbesondere das Vulkandiagramm und das Venn-Diagramm verwendet werden. limma identifizierte 3323 DEGs zwischen dem CHOL und normalem Gewebe mit den |logFC|≥2 und adj. P.Val <0,05 als Schwellenwerte, darunter 1880 in CHOL-Geweben herunterreguliert und 1443 hochreguliert wurden (Abbildung 1a). In der Zwischenzeit identifizierte edgeR die 1578 herunterregulierten DEGs und 3121 hochregulierten DEGs (Abbildung 1b); DESeq2 identifizierte die 1616 herunterregulierten DEGs und 2938 hochregulierten DEGs (Abbildung 1c). Vergleicht man die Ergebnisse dieser drei Methoden, so überschnitten sich 1431 hochregulierte DEGs und 1531 herunterregulierte DEGs (Abbildung 2).
Abbildung 1. Identifizierung differentiell exprimierter Gene (DEGs) zwischen CHOL und normalem Gewebe. (a-c) Die Vulkandiagramme aller gene, die von limma, edgeR bzw. DESeq2 erworben wurden, adj p-Wert (-log10) wird gegen die Faltenänderung (log2) aufgezeichnet, rote Punkte repräsentieren die hochregulierten DEGs (angepasster p-Wert<0,05 und log | FC|> 2) und die grünen Punkte repräsentieren die nach unten regulierten DEGs (adjustierter p-Wert< 0,05 und log | FC|< 2). Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.
Abbildung 2. Venn-Diagramme zeigen Überschneidungen zwischen den Ergebnissen, die aus limma, edgeR und DESeq2 abgeleitet wurden. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.
Ergänzende Dateien. Bitte klicken Sie hier, um diese Datei herunterzuladen.
Subscription Required. Please recommend JoVE to your librarian.
Discussion
Reichlich vorhandene Aberrat-Transkripte bei Krebserkrankungen können leicht durch RNA-seq-Differentialanalyse identifiziert werden5. Die Anwendung der RNA-seq-Differentialexpressionsanalyse ist jedoch oft eingeschränkt, da sie bestimmte Kenntnisse mit der R-Sprache und die Fähigkeit, geeignete Methoden zu wählen, erfordert. Um dieses Problem anzugehen, bieten wir eine detaillierte Einführung in die drei bekanntesten Methoden (limma, EdgeR und DESeq2) und Tutorials zur Anwendung der RNA-seq-Differentialexpressionsanalyse. Dies wird das Verständnis der Gemeinsamkeiten und Unterschiede zwischen allen drei Methoden erleichtern, die Auswahl einer geeigneten Methode für einzelne Daten ermöglichen und es uns ermöglichen, die komplexen dynamischen biologischen Prozesse zu verstehen.
Hier präsentieren wir ein detailliertes Protokoll für die RNA-seq-Differentialexpressionsanalyse durch Limma, EdgeR bzw. DESeq2 in fünf Stufen: (i) Herunterladen und Vorverarbeitung von Daten, (ii-iv) differentielle Expressionsanalyse durch Limma, EdgeR bzw. DESeq2, (v) Vergleich der Ergebnisse dieser drei Methoden durch ein Venn-Diagramm.
Die drei Methoden haben ähnliche und unterschiedliche Schritte zwischen den Prozessen der differentialen Expressionsanalyse. Für die Statistik in Limma wird ein lineares Modell verwendet, das für alle Genexpressionstechnologien anwendbar ist, einschließlich Microarrays, RNA-seq und quantitativer PCR8,13, während edgeR und DESeq2 eine Reihe statistischer Methoden implementieren, die auf der negativen Binomialverteilung9,10basieren, und edgeR und DESeq2 für RNA-seq-Daten geeignet sind. Darüber hinaus sind die normalisierten RNA-seq-Zähldaten für EdgeR und Limma notwendig, während DESeq2 seine eigenen Bibliotheksdiskrepanzen verwendet, um Daten anstelle von Normalisierung zu korrigieren, und die Daten in DESeq2 eine ganzzahlige Matrix sein müssen. Die Normalisierungsmethoden umfassen TMM (getrimmtes Mittel der M-Werte), TMMwsp, RLE (relative Log-Expression) und oberes Quartil, unter denen TMM die am häufigsten verwendete Normalisierungsmethode für RNA-Seq-Daten ist. Die Ergebnisse der drei Methoden zeigten, dass DESeq2 und EdgeR mehr DEGs erhalten als Limma. Der Grund für diesen Unterschied ist, dass edgeR und DESeq2 auf dem negativen Binomialmodell basieren, das zu einer großen Anzahl von Fehlalarmen beiträgt. Im Gegenteil, limma-voom verwendet nur die Varianzfunktion und zeigt keine übermäßigen False Positives, wie es bei einer varianzstabilisierenden Transformation der Fall ist, gefolgt von einer linearen Modellanalyse mit Limma14,15,16.
Alle drei Methoden haben ihre eigenen Vorteile, und die Auswahl hängt nur von der Art der Daten ab. Wenn beispielsweise Microarray-Daten vorhanden sind, sollte Limma mit Priorität angegeben werden, aber wenn es sich um die Sequenzierungsdaten der nächsten Generation handelt, werden DESeq2 und EdgeR bevorzugt9,10,17. Zusammenfassend stellen wir hier ein detailliertes Protokoll für die RNA-seq-Differentialexpressionsanalyse mit den R-Paketen limma, edgeR und DESeq2 zur Verfügung. Die Output-Ergebnisse der drei Methoden überschneiden sich teilweise, und diese differentiellen Methoden haben ihre jeweiligen Vorteile. Leider deckt dieses Protokoll nicht die technischen Details für andere Datentypen (z. B. Microarray-Daten) und Methoden (z. B. EBSeq)18ab.
Subscription Required. Please recommend JoVE to your librarian.
Disclosures
Das Manuskript wurde noch nicht veröffentlicht und wird nicht für eine Veröffentlichung an anderer Stelle in Betracht gezogen. Alle Autoren haben zur Erstellung dieses Manuskripts für wichtige intellektuelle Inhalte beigetragen und das endgültige Manuskript gelesen und genehmigt. Wir erklären, dass kein Interessenkonflikt besteht.
Acknowledgments
Diese Arbeit wurde von der National Natural Science Foundation of China (Grant No. 81860276) und Key Special Fund Projects des National Key R&D Program (Grant No. 2018YFC1003200) unterstützt.
Materials
| Name | Company | Catalog Number | Comments |
| R | version 3.6.2 | free software | |
| Rstudio | free software |
References
- Tambonis, T., Boareto, M., Leite, V. B. P. Differential Expression Analysis in RNA-seq Data Using a Geometric Approach. Journal of Computational Biology. 25, 1257-1265 (2018).
- Wang, Z., Gerstein, M., Snyder, M. RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews. Genetics. 10, 57-63 (2009).
- Anders, S., et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nature Protocols. 8, 1765-1786 (2013).
- McDermaid, A., Monier, B., Zhao, J., Liu, B., Ma, Q. Interpretation of differential gene expression results of RNA-seq data: review and integration. Briefings in Bioinformatics. 20, 2044-2054 (2019).
- Costa-Silva, J., Domingues, D., Lopes, F. M. RNA-Seq differential expression analysis: An extended review and a software tool. PloS One. 12, 0190152 (2017).
- Law, C. W., et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research. 5, (2016).
- Varet, H., Brillet-Guéguen, L., Coppée, J. Y., Dillies, M. A. SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data. PloS One. 11, 0157022 (2016).
- Ritchie, M. E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 43, 47 (2015).
- Robinson, M. D., McCarthy, D. J., Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26, Oxford, England. 139-140 (2010).
- Love, M. I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 15, 550 (2014).
- Gentleman, R. C., et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biology. 5, 80 (2004).
- Law, C. W., Chen, Y., Shi, W., Smyth, G. K. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology. 15, 29 (2014).
- Smyth, G. K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology. 3, (2004).
- Lund, S. P., Nettleton, D., McCarthy, D. J., Smyth, G. K. Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates. Statistical Applications in Genetics and Molecular Biology. 11, (2012).
- Reeb, P. D., Steibel, J. P. Evaluating statistical analysis models for RNA sequencing experiments. Frontiers in Genetics. 4, 178 (2013).
- Rocke, D. M., et al. Excess False Positive Rates in Methods for Differential Gene Expression Analysis using RNA-Seq Data. bioRxiv. , (2015).
- Agarwal, A., et al. Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays. BMC genomics. 11, 383 (2010).
- Leng, N., et al. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 29, Oxford, England. 1035-1043 (2013).
