Summary
Een gedetailleerd protocol van differentiële expressie analysemethoden voor RNA-sequencing werd verstrekt: limma, EdgeR, DESeq2.
Abstract
RNA-sequencing (RNA-seq) is een van de meest gebruikte technologieën in transcriptomics omdat het de relatie tussen de genetische verandering en complexe biologische processen kan onthullen en grote waarde heeft in diagnostiek, prognose en therapieën van tumoren. Differentiële analyse van RNA-seq-gegevens is cruciaal om afwijkende transcripties te identificeren, en limma, EdgeR en DESeq2 zijn efficiënte tools voor differentiële analyse. RNA-seq differentiële analyse vereist echter bepaalde vaardigheden met R-taal en het vermogen om een geschikte methode te kiezen, die ontbreekt in het curriculum van medisch onderwijs.
Hierin bieden we het gedetailleerde protocol om differentieel uitgedrukte genen (DEG's) tussen cholangiocarcinoom (CHOL) en normale weefsels te identificeren via respectievelijk limma, DESeq2 en EdgeR, en de resultaten worden weergegeven in vulkaanpercelen en Venn-diagrammen. De drie protocollen van limma, DESeq2 en EdgeR zijn vergelijkbaar, maar hebben verschillende stappen tussen de processen van de analyse. Een lineair model wordt bijvoorbeeld gebruikt voor statistieken in limma, terwijl de negatieve binomiale verdeling wordt gebruikt in edgeR en DESeq2. Bovendien zijn de genormaliseerde RNA-seq-tellingsgegevens noodzakelijk voor EdgeR en limma, maar niet nodig voor DESeq2.
Hier bieden we een gedetailleerd protocol voor drie differentiële analysemethoden: limma, EdgeR en DESeq2. De resultaten van de drie methoden overlappen elkaar gedeeltelijk. Alle drie de methoden hebben hun eigen voordelen en de keuze van de methode hangt alleen af van de gegevens.
Introduction
RNA-sequencing (RNA-seq) is een van de meest gebruikte technologieën in transcriptomics met veel voordelen (bijv. hoge reproduceerbaarheid van gegevens) en heeft ons begrip van de functies en dynamiek van complexe biologische processen drastisch vergroot1,2. Identificatie van aberraattranscripties onder verschillende biologische context, die ook bekend staan als differentieel uitgedrukte genen (DEG's), is een belangrijke stap in de RNA-seq-analyse. RNA-seq maakt het mogelijk om een diepgaand begrip te krijgen van pathogenese gerelateerde moleculaire mechanismen en biologische functies. Daarom is differentiële analyse als waardevol beschouwd voor diagnostiek, prognose en therapieën van tumoren3,4,5. Momenteel zijn er meer open-source R/Bioconductor-pakketten ontwikkeld voor RNA-seq differentiële expressieanalyse, met name limma, DESeq2 en EdgeR1,6,7. Differentiële analyse vereist echter bepaalde vaardigheden met R-taal en het vermogen om de juiste methode te kiezen, die ontbreekt in het curriculum van medisch onderwijs.
In dit protocol, gebaseerd op de cholangiocarcinoom (CHOL) RNA-seq telling gegevens geëxtraheerd uit The Cancer Genome Atlas (TCGA), drie van de meest bekende methoden (limma8, EdgeR9 en DESeq210) werden uitgevoerd, respectievelijk, door het R-programma11 om de DEG's tussen CHOL en normale weefsels te identificeren. De drie protocollen van limma, EdgeR en DESeq2 zijn vergelijkbaar, maar hebben verschillende stappen tussen de processen van de analyse. De genormaliseerde RNA-seq-tellingsgegevens zijn bijvoorbeeld nodig voor EdgeR en limma8,9, terwijl DESeq2 zijn eigen bibliotheekverschillen gebruikt om gegevens te corrigeren in plaats van normaliseren10. Verder is edgeR specifiek geschikt voor RNA-seq data, terwijl de limma wordt gebruikt voor microarrays en RNA-seq. Limma keurt een lineair model goed om de DEG 's12te beoordelen , terwijl de statistieken in edgeR zijn gebaseerd op de negatieve binomiale verdelingen, waaronder empirische Bayes-schatting, exacte tests, gealdaliseerde lineaire modellen en quasi-waarschijnlijkheidstests9.
Samengevat bieden we de gedetailleerde protocollen van RNA-seq differentiële expressieanalyse met behulp van respectievelijk limma, DESeq2 en EdgeR. Door naar dit artikel te verwijzen, kunnen gebruikers eenvoudig de RNA-seq differentiële analyse uitvoeren en de juiste differentiële analysemethoden voor hun gegevens kiezen.
Subscription Required. Please recommend JoVE to your librarian.
Protocol
OPMERKING: Open het R-studio programma en laad R bestand "DEGs.R", het bestand kan worden verkregen uit aanvullende bestanden / scripts.
1. Downloaden en voorverwerking van gegevens
- Download de HTSeq-count-gegevens (high-throughput sequencing) van cholangiocarcinoom (CHOL) uit The Cancer Genome Atlas (TCGA). Deze stap kan eenvoudig worden bereikt met de volgende R-code.
- Klik op Uitvoeren om R-pakketten te installeren.
- Klik op Uitvoeren om R-pakketten te laden.
als(!requireNamespace("BiocManager", quiet=TRUE))
+ install.packages("BiocManager")
BiocManager::install(c("TCGAbiolinks", "SamengevatExperiment")) - Stel de werkmap in.
bibliotheek (TCGAbiolinks)
bibliotheek(SamengevatExperiment)
setwd("C:/Gebruikers/LIUSHIYI/Desktop") - Kies het kankertype.
kanker <- "TCGA-CHOL" - Voer de R-code uit vanuit het bestand "GDCquery.R" om de gegevens te downloaden. Het bestand "GDCquery.R" kan worden verkregen uit aanvullende bestanden / scripts:
bron("Aanvullende bestanden/Scripts/GDCquery.R")
hoofd(cnt)
##TCGA-3X-AAVA-01A-11R-A41I-07
##ENSG00000000003 4262
##ENSG00000000005 1
##ENSG00000000419 1254
##ENSG00000000457 699
##ENSG00000000460 239
##ENSG00000000938 334
OPMERKING: Na uitvoering worden de CHOLHTSeq-tellingsgegevens gedownload en "cnt" genoemd, waarbij rijen ensemblegen-ID's vertegenwoordigen en kolommen voorbeeld-ID's vertegenwoordigen. Let op de nummers op de posities 14-15 in de monster-ID's; getallen variërend van 01 tot 09 wijzen op tumoren en variërend van 10 tot 19 wijzen op normale weefsels.
- Converteer ensemblegen-ID's naar genrembolen.
- Importeer het annotatiebestand in R op basis van het opslagpad. Het annotatiebestand (gencode.v22.annotation.gtf) kan worden verkregen uit aanvullende bestanden.
gtf_v22 <- rtracklayer::import('Aanvullende bestanden/gencode.v22.annotation.gtf') - Voer de R-code uit vanuit de gtf_v22. R"-bestand, dat kan worden verkregen uit aanvullende bestanden/ scripts:
bron("Aanvullende bestanden/Scripts/gtf_v22. R") - Pas de functie "ann" toe om de ensemblegen-ID's om te zetten in genrembolen.
cnt=ann(cnt,gtf_v22)
- Importeer het annotatiebestand in R op basis van het opslagpad. Het annotatiebestand (gencode.v22.annotation.gtf) kan worden verkregen uit aanvullende bestanden.
- Laag-uitgedrukte genen filteren
- Klik op Uitvoeren om het R-pakket "edgeR" te installeren.
BiocManager::installeren("edgeR") - Klik op Uitvoeren om het R-pakket "edgeR" te laden.
bibliotheek(edgeR) - Voer de volgende R-code uit om genen met een CPM-waarde (Counts Per Million) groter dan één op ten minste twee monsters te behouden.
houd <-rowSums(cpm(cnt)>1)>=2
cnt <- as.matrix(cnt[houden,])
OPMERKING: De waarde tellingen per miljoen (CPM) wordt gebruikt in plaats van de leestelling om de afwijking veroorzaakt door verschillende sequencingdiepten te elimineren.
- Klik op Uitvoeren om het R-pakket "edgeR" te installeren.
2. Differentiële expressieanalyse door "limma"
- Klik op Uitvoeren om het R-pakket "limma" te installeren.
BiocManager::installeren("limma") - Klik op Uitvoeren om de R-pakketten "limma", "edgeR" te laden.
bibliotheek(limma)
bibliotheek(edgeR) - Voer de volgende R-code uit om de ontwerpmatrix te maken.
groep <- substring(colnames(cnt),14,15) # Extract group information
groep [groep %in% "01"] <- "Cancer" # set '01' as tumor tissue
groep [groep %in% "11"] <- "Normal" # set '11' as normal tissue
groep <- factor (group, levels = c("Normal","Cancer"))- Maak de ontwerpmatrix.
ontwerp <- model.matrix (~group)
rijnamen(ontwerp) <- colnames(cnt) - Maak het object DGEList.
dge <- DGEList(tellingen = cnt, groep = groep) - Normaliseer de gegevens.
dge <- calcNormFactors(dge, methode = "TMM") - Voer de volgende R-code uit om de op limma-trendmethode gebaseerde differentiële expressieanalyse uit te voeren.
dge
##An object van klasse "DGEList"
##$counts
##TCGA-3X-AAVA-01A-11R-A41I-07
##TSPAN6 4262
##DPM1 1254
##SCYL3 699
##C1orf112 239
##FGR 334 - Bereken de CPM-waarde.
logdge <-cpm(dge, log=TRUE, prior.count=3) - Klik op Uitvoeren om een lineair model te passen om de gegevens te voorspellen of de relatie tussen variabelen af te leiden.
fit <- lmFit (logdge, ontwerp) - Bereken de T-waarde, F-waarde en log-odds op basis van Bayesian.
fit <- eBayes (fit, trend=WAAR) - De resultaattabel extraheren.
res_limma<- as.data.frame (topTabel(fit,n=Inf))
hoofd(res_limma)
## logFC AveExpr t P.Waarde 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 resultaat van differentiële expressieanalyse wordt opgeslagen in "res_limma", waaronder de gen-id, log2 fold change value (logFC), het gemiddelde log2-expressieniveau van het gen in het experiment (AveExpr), de gewijzigde t-statistiek (t), relavent p-waarde (P.Value), de false discovery rate (FDR) gecorrigeerde p-waarde (adj. P.Val) en de log-odds van differentieel uitgedrukte genen (B)
OPMERKING: De functie "calcNormFactors()" van de "edgeR" werd gebruikt om de gegevens te normaliseren om de invloed te elimineren die wordt veroorzaakt door monstervoorbereiding of bibliotheekconstructie en sequencing. Bij de constructie van de ontwerpmatrix is het noodzakelijk om experimenteel ontwerp (bijv. weefseltype: normaal of tumorweefsel) te matchen om ID's van de matrix te bemonsteren. limma-trend is geschikt voor gegevens waarvan de sequencingdiepte hetzelfde is, terwijl limma-voom geschikt is: (i) wanneer de grootte van de voorbeeldbibliotheek anders is; ii) gegevens die niet door TMM zijn genormaliseerd; (iii) er is veel "ruis" in de gegevens. Een positieve logFC betekent dat het gen in het experiment up-gereguleerd is, terwijl negatief getal betekent dat het gen down-gereguleerd is. - Identificeer de DEG's.
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')) # De adj.p Value < 0,05 en de |log2FC| >= 2 zijn drempels om de DEG's te identificeren
samenvatting(res_limma$sig)
##down niet omhoog
##1880 17341 1443 - Voer de resultaattabel uit naar een bestand.
schrijven.csv(res_limma, bestand = 'result_limma.csv') - Klik op Uitvoeren om het R-pakket "ggplot2" te installeren.
install.packages("ggplot2") - Klik op Uitvoeren om het R-pakket "ggplot2" te laden.
bibliotheek(ggplot2) - Voer de R-code uit vanaf de vulkaan. R" om het vulkaanperceel te creëren. Het bestand "vulkaan. R" kan worden verkregen uit aanvullende bestanden.
bron("Aanvullende bestanden/Scripts/volcano. R")
vulkaan(res_limma,"logFC","adj. P.Val",2.0.05)
OPMERKING: Genen kunnen worden toegewezen aan verschillende posities op basis van hun log2FC- en adj-p-waarden, de omhoog gereguleerde DEG's zijn rood gekleurd en de down-gereguleerde DEG's zijn groen gekleurd. - Klik op Exporteren om het vulkaanplot op te slaan.
OPMERKING: De vulkaanpercelen kunnen worden gegenereerd en gedownload in verschillende formaten (bijv. pdf, TIFF, PNG, JPEG-formaat). Genen kunnen worden toegewezen aan verschillende posities op basis van hun log2FC- en adj p-waarden, de up-regulated DEG's (log2FC > 2, adj p < 0,05) zijn rood gekleurd en de down-gereguleerde DEG's (log2FC < -2, adj p < 0,05) zijn groen gekleurd, niet-DEG's zijn grijs gekleurd.
- Maak de ontwerpmatrix.
3. Differentiële expressieanalyse door middel van "edgeR"
- Klik op Uitvoeren om het R-pakket "edgeR" te laden.
bibliotheek(edgeR) - Voer de volgende R-code uit om een ontwerpmatrix te maken.
groep <-substring(colnames(cnt),14,15)
groep [groep %in% "01"] <- "Kanker"
groep [groep %in% "11"] <- "Normaal"
group=factor(groep, niveaus = c("Normaal","Kanker"))
ontwerp <-model.matrix(~group)
rijnamen(ontwerp) = colnames(cnt) - Klik op Uitvoeren om het DGEList-object te maken.
dge <- DGELijst(counts=cnt) - Normaliseer de gegevens.
dge <- calcNormFactors(dge, methode = "TMM") - Klik op Uitvoeren om de spreiding van genexpressiewaarden te schatten.
dge <- schattingDisp(dge, design, robuust = T) - Klik op Uitvoeren om het model aan te passen om gegevens te tellen.
fit <- glmQLFit(dge, ontwerp) - Voer een statistische test uit.
pasvorm <- glmQLFTest(fit) - De resultaattabel extraheren. Het resultaat wordt opgeslagen in "res_edgeR", waaronder de wijzigingswaarde van de logvouwen, log CPM, F, p-waarde en FDR gecorrigeerde p-waarde.
res_edgeR=as.data.frame(topTags(fit, n=Inf))
hoofd(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 resultaat wordt opgeslagen in "res_edgeR", waaronder de log fold change value (logFC), log CPM, F, p-waarde en FDR gecorrigeerde p-waarde - Identificeer de DEG's.
res_edgeR$sig = as.factor(
ifelse (res_edgeR$FDR < 0.05 & abs(res_edgeR$logFC) > 2,
ifelse(res_edgeR$logFC > 2 'up','down'),'not'))
samenvatting(res_edgeR$sig)
##down niet omhoog
##1578 15965 3121 - Voer de resultaattabel uit naar een bestand.
schrijven.csv(res_edgeR, bestand = 'res_edgeR.csv') - Maak het vulkaanperceel.
vulkaan (res_edgeR,"logFC","FDR",2,0.05) - Klik op Exporteren om het vulkaanplot op te slaan.
4. Differentiële expressieanalyse via "DESeq2"
- Klik op Uitvoeren om R-pakketten "DESeq2" te installeren.
BiocManager::installeren("DESeq2") - Klik op Uitvoeren om R-pakketten "DESeq2" te laden.
bibliotheek(DESeq2) - Voer de volgende R-code uit om de groeperingsfactor te bepalen.
groep <-substring(colnames(cnt),14,15)
groep [groep %in% "01"] <- "Kanker"
groep [groep %in% "11"] <- "Normaal"
group=factor(groep, niveaus = c("Normaal","Kanker")) - Maak het OBJECT DESeqDataSet.
dds <-DESeqDataSetFromMatrix (cnt, DataFrame(groep), ontwerp = ~group)
Dds
##class: DESeqDataSet
##dim: 20664 45
##metadata(1): versie
##assays(1): telt
##rownames(20664): TSPAN6 DPM1 ... RP11-274B21.13 LINC01144
##rowData namen(0):
##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ...
##colData namen(1): groep - Voer de analyse uit.
dds <- DESeq(dds) - Genereer de resultatentabel.
res_DESeq2 <- data.frame(resultaten(dds))
hoofd(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
OPMERKING: Het resultaat wordt opgeslagen in "res_DESeq2", dat het gemiddelde van het genormaliseerde aantal lees (baseMean), log fold Change value (log2FoldChange), log fold change standard error (lfcSE), de Wald-statistiek (stat), de oorspronkelijke p-waarde (pvalue) en de gecorrigeerde p-waarde (padj) omvat - Identificeer DEG's.
res_DESeq2$sig = as.factor(
ifelse (res_DESeq2$padj < 0.05 & abs (res_DESeq2$log2FoldChange) > 2,
ifelse(res_DESeq2$log2FoldChange > 2 'up','down'),'not'))
samenvatting(res_DESeq2$sig)
##down niet omhoog
##1616 16110 2938 - Voer de resultaattabel uit naar een bestand.
schrijven.csv(res_DESeq2, bestand = 'res_DESeq2.csv') - Maak het vulkaanperceel.
vulkaan (res_DESeq2,"log2FoldChange","padj",2,0.05) - Klik op Exporteren om het vulkaanplot op te slaan.
5. Venn diagram
- Klik op Uitvoeren om het R-pakket "VennDiagram" te installeren.
install.packages("VennDiagram") - Klik op Uitvoeren om het R-pakket "VennDiagram" te laden.
bibliotheek (VennDiagram) - Maak een Venn-diagram van omhoog gereguleerde DEG's.
grid.newpage()
grid.draw(venn.diagram(list(Limma=rijnamen(res_
limma[res_limma$sig=="omhoog",]),
edgeR=rownames(res_edgeR[res_edgeR$sig=="omhoog",]),
DESeq2=rijnamen(res_DESeq2[res_DESeq2$sig==
"omhoog",])),
NULL,hoogte = 3,breedte = 3,eenheden = "in",
col="zwart",lwd=0,3,fill=c("#FF6666","#FFFF00",
"#993366"),
alfa=c(0,5, 0,5, 0,5),main = "Up-regulated DEGs")) - Klik op Exporteren om het Venn-diagram op te slaan.
- Maak een Venn-diagram van down gereguleerde DEG's.
grid.newpage()
grid.draw(venn.diagram(list(Limma=rijnamen(res_
limma[res_limma$sig=="omlaag",]),
edgeR=rijnamen(res_edgeR[res_edgeR$sig==
"omlaag",]),
DESeq2=rijnamen(res_DESeq2[res_DESeq2$sig=="omlaag",])),
NULL,hoogte = 3,breedte = 3,eenheden = "in",
col="zwart",lwd=0,3,fill=c("#FF6666","#FFFF00",
"#993366"),
alfa=c(0,5, 0,5, 0,5),main = "Down-regulated DEGs")) - Klik op Exporteren om het Venn-diagram op te slaan.
Subscription Required. Please recommend JoVE to your librarian.
Representative Results
Er zijn verschillende benaderingen om het resultaat van differentiële expressieanalyse te visualiseren, waaronder het vulkaanplot en venndiagram in het bijzonder worden gebruikt. limma identificeerde 3323 DEG's tussen de CHOL en normale weefsels met de |logFC|≥2 en adj. P.Val <0,05 als drempels, waaronder 1880 in CHOL-weefsels en 1443 in CHOL-weefsels werden gereguleerd (figuur 1a). Ondertussen identificeerde edgeR de 1578 down-gereguleerde DEG's en 3121 up-regulated DEG's (Figuur 1b); DESeq2 identificeerde de 1616 down-gereguleerde DEG 's en 2938 up-regulated DEG 's (Figuur 1c). Door de resultaten van deze drie methoden te vergelijken, werden 1431 up-regulated DEG's en 1531 down-regulated DEG's overlappend (Figuur 2).
Figuur 1. Identificatie van differentieel uitgedrukte genen (DEG's) tussen CHOL en normale weefsels. (a-c) De vulkaanplots van alle genen verworven door respectievelijk limma, edgeR en DESeq2, adj p-waarde (-log10) worden uitgezet tegen de vouwverandering (log2), rode punten vertegenwoordigen de up-regulated DEGs (aangepaste p-waarde<0,05 en log | FC|> 2) en de groene punten vertegenwoordigen de naar beneden gereguleerde DEG's (aangepaste p-waarde< 0,05 en log | FC|< 2). Klik hier om een grotere versie van deze afbeelding te bekijken.
Figuur 2. Venn diagrammen tonen overlap tussen de resultaten afgeleid van de limma, edgeR en DESeq2. Klik hier om een grotere versie van deze afbeelding te bekijken.
Aanvullende bestanden. Klik hier om dit bestand te downloaden.
Subscription Required. Please recommend JoVE to your librarian.
Discussion
Overvloedige afwijkende transcripties bij kankers kunnen gemakkelijk worden geïdentificeerd door RNA-seq differentiële analyse5. De toepassing van RNA-seq differentiële expressieanalyse is echter vaak beperkt omdat het bepaalde vaardigheden met R-taal vereist en het vermogen om geschikte methoden te kiezen. Om dit probleem aan te pakken, bieden we een gedetailleerde inleiding tot de drie meest bekende methoden (limma, EdgeR en DESeq2) en zelfstudies voor het toepassen van de RNA-seq differentiële expressieanalyse. Dit zal het begrip van de overeenkomsten en verschillen tussen alle drie de methoden vergemakkelijken, de selectie van een geschikte methode voor individuele gegevens mogelijk maken en ons in staat stellen de complexe dynamische biologische processen te begrijpen.
Hier presenteren we een gedetailleerd protocol voor RNA-seq differentiële expressieanalyse via respectievelijk limma, edgeR en DESeq2, in vijf fasen: (i) downloaden en voorverwerking van gegevens, (ii-iv) differentiële expressieanalyse via respectievelijk limma, edgeR en DESeq2, (v) vergelijking van de resultaten van deze drie methoden door middel van een Venn-diagram.
De drie methoden hebben vergelijkbare en verschillende stappen tussen de processen van de differentiële expressieanalyse. Een lineair model wordt gebruikt voor statistieken in limma, dat van toepassing is op alle genexpressietechnologieën, waaronder microarrays, RNA-seq en kwantitatieve PCR8,13, terwijl edgeR en DESeq2 een reeks statistische methodologieën implementeren op basis van de negatieve binomiale verdeling9,10, en edgeR en DESeq2 zijn geschikt voor RNA-seq-gegevens. Bovendien zijn de genormaliseerde RNA-seq-tellingsgegevens nodig voor EdgeR en limma, terwijl DESeq2 zijn eigen bibliotheekverschillen gebruikt om gegevens te corrigeren in plaats van te normaliseren en de gegevens in DESeq2 een gehele matrix moeten zijn. De normalisatiemethoden omvatten TMM (bijgesneden gemiddelde van M-waarden), TMMwsp, RLE (relatieve logboekexpressie) en bovenkwartiel, waaronder TMM de meest gebruikte normalisatiemethode voor RNA-seq-gegevens is. De resultaten van de drie methoden toonden aan dat DESeq2 en EdgeR meer DEG's verkrijgen dan limma. De reden voor dit verschil is dat edgeR en DESeq2 zijn gebaseerd op het negatieve binomiale model, dat bijdraagt aan grote aantallen valse positieven. Integendeel, limma-voom gebruikt alleen de variantiefunctie en vertoont geen overmatige valse positieven, zoals het geval is bij een variantie stabiliserende transformatie gevolgd door lineaire modelanalyse met limma14,15,16.
Alle drie de methoden hebben hun eigen voordelen en de keuze is gewoon afhankelijk van het type gegevens. Als er bijvoorbeeld microarray-gegevens zijn, moet limma prioriteit krijgen, maar wanneer het de volgende generatie sequencinggegevens is, hebben DESeq2 en EdgeR de voorkeur9,10,17. Samengevat bieden we hier een gedetailleerd protocol voor RNA-seq differentiële expressieanalyse met respectievelijk R-pakketten limma, edgeR en DESeq2. De outputresultaten van de drie methoden overlappen elkaar gedeeltelijk en deze differentiële methoden hebben hun respectieve voordelen. Helaas heeft dit protocol geen betrekking op de technische details voor andere gegevenstypen (bv. microarraygegevens) en methoden (bv. EBSeq)18.
Subscription Required. Please recommend JoVE to your librarian.
Disclosures
Het manuscript is nog niet eerder gepubliceerd en wordt niet elders gepubliceerd. Alle auteurs hebben bijgedragen aan de creatie van dit manuscript voor belangrijke intellectuele inhoud en hebben het uiteindelijke manuscript gelezen en goedgekeurd. Wij verklaren dat er geen sprake is van belangenverstrengeling.
Acknowledgments
Dit werk werd ondersteund door de National Natural Science Foundation of China (Grant No. 81860276) en Key Special Fund Projects of National Key R&D Program (Grant No. 2018YFC1003200).
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).
