Rup (RNA-seq Usability Assessment Pipeline) - Kwaliteitscontrole voor bulk RNA-seq-experimenten in eukaryoten

Samenvatting

Dit protocol maakt een eerste kwaliteitscontrole mogelijk voor RNA-seq-experimenten voor wet-lab biologen met beperkte bio-informatica-ervaring.

Samenvatting

Moderne benaderingen in de moleculaire plantenwetenschap vereisen vaak bulk RNA-seq-experimenten, bijvoorbeeld om wereldwijde veranderingen in transcriptomen bij behandelingen te volgen of om belangrijke componenten van regulerende routes te identificeren. Bijgevolg vertrouwen verschillende gebieden van de plantenwetenschappen op hoogwaardige en reproduceerbare bulk RNA-seq-gegevens voor wetenschappelijke vooruitgang. Uit onze ervaring blijkt echter dat kennis over en toepassing van kwaliteitsbeheersingsmaatregelen in RNA-seq-datasets vaak ontbreekt. Hier introduceren we Rup (RNA-seq usability assessment pipeline) voor de kwaliteitscontrole van bulk RNA-seq-gegevens, voor daaropvolgende genexpressieanalyses, die op zichzelf staat en direct toepasbaar is voor wet-lab biologen met basiskennis van R. Rup helpt om onderscheid te maken tussen sequentiegegevens van hoge kwaliteit, geschikt voor downstream genexpressie-experimenten, en gegevens die ongeschikt zijn voor algemene verdere analyse. Rup omvat tests voor verschillende vaak voorkomende problemen, zoals onvoldoende leesnummers of mapping, identificatie van contaminaties, kwantificering van rRNA-fracties in de totale RNA-seq-gegevens, replicatietesten op gelijkenis, het gebruik van echte gegevens voor demonstratie en het aanbieden van intuïtieve visualisatie. Rup biedt een reeks hulpmiddelen om experimentele tekortkomingen te identificeren vóór gestandaardiseerde transcriptoomanalyse, waardoor de gegevenskwaliteit voor individuele onderzoekers en het veld wordt verbeterd. Dit vergroot het vertrouwen in de analyse van bulk-RNA-seq-gegevens en biedt een basis voor toekomstige richtlijnen die minimale kwaliteitscontrolecriteria definiëren, waardoor de betrouwbaarheid en transparantie van gepubliceerde RNA-seq-gegevens wordt verbeterd. Rup- en testgegevens zijn beschikbaar op https://github.com/oliverrupp/rup.

Inleiding

Transcriptomics (RNA-seq) experimenten ondervragen uitgebreid transcriptionele handtekeningen die fenotypes vormgeven. Deze benadering werd onvervangbaar in de moleculaire genetica van planten om enkelvoudige of co-tot expressie gebrachte genen en biologische processen te identificeren die bijvoorbeeld betrokken zijn bij ontwikkeling, interactie tussen plantpathogenen of abiotische stressresistentie^{1, 2 ,3}. Recente ontwikkelingen in RNA-seq-technologie hebben de specificiteit verhoogd en de detectie van verschillende isovormen en varianten met single-base resolutie mogelijk gemaakt, waardoor sequentievariaties van grotere indels tot single-nucleotide polymorfismen (SNP's) kunnen worden geïdentificeerd. De gegevens die door RNA-seq worden verkregen, worden gekenmerkt door een breed dynamisch bereik, waardoor de detectie van zowel zeer overvloedige als laag tot expressie gebrachte transcripten wordt vergemakkelijkt, en vereisen geschikte experimentele omstandigheden voor consistentie. Bovendien zijn RNA-seq-datasets groot en dus rekenintensief om te analyseren en op te slaan, wat een efficiënt gegevensbeheer en uitgebreide rekenbronnen met zich meebrengt. RNA-seq vereist input van hoge kwaliteit bij elke stap van de workflow, aangezien zwakke punten in elke fase de resultaten kunnen verspreiden en in gevaar kunnen brengen. Slechte RNA-integriteit of technische problemen tijdens de voorbereiding van de bibliotheek leiden tot vertekening en een vermindering van de nauwkeurigheid. Parameters voor het verkrijgen van onbewerkte lezingen van hoge kwaliteit kunnen variëren tussen sequencingfaciliteiten van de volgende generatie (NGS).

Volgens onze ervaring raden we aan om alleen RNA te gebruiken met een RNA-integriteitsnummer (RIN) van meer dan 7, wat een indicatie is van een grotendeels intacte mRNA-structuur, als input voor de voorbereiding van de RNA-seq-bibliotheek. Voor een succesvolle sequentiebepaling is ongeveer 2 μg totaal RNA met een concentratie van 50-200 ng/μL vereist voor standaard bibliotheekvoorbereidingsprotocollen. De zuiverheid van het RNA moet worden bevestigd door een OD_{260/280-verhouding} tussen 1,8 en 2,1 en een OD_{260/230-verhouding} van meer dan 1,5 met behulp van een spectrofotometer. Onvoldoende sequencingdiepte, lage mappingsnelheden of verkeerde uitlijning met het referentiegenoom kunnen de genexpressie en splitsingskwantificering verder verstoren. Bovendien kan een ontoereikende experimentele opzet, zoals een gebrek aan discriminatie tussen monsters van verschillende weefsels of behandelingen en een grote variabiliteit tussen duplo's, ruis veroorzaken en de reproduceerbaarheid verminderen. Hoewel veel laboratoria transcriptomen routinematig analyseren, worden strenge kwaliteitscontroles van de eerste essentiële analysestappen vaak niet gerapporteerd of zijn ze volledig afwezig in publicaties. Dit kan leiden tot overinterpretatie van de resultaten van transcriptoomanalyse en bijgevolg tot niet-reproduceerbare resultaten.

Hier bieden we een workflow voor kwaliteitscontrole van de eerste stappen die nodig zijn voor hoogwaardige transcriptoomanalyses van mRNA-profielen om transcriptionele veranderingen te meten. Ons doel is om wet-lab biologen met beperkte kennis in bio-informatica in staat te stellen hun primaire transcriptomics-gegevens te evalueren. Rup (Figuur 1) is toegankelijk voor onderzoekers die bekend zijn met basiskennis van R. Het uitvoeren van de hier gepresenteerde workflow geeft onderzoekers een gedetailleerd inzicht in hun primaire gegevens, inclusief hun mogelijke beperkingen voor latere analyse. Voor zover wij weten, ontbreekt tot nu toe een praktische primaire transcriptoomgegevensevaluatiepijplijn in combinatie met richtlijnen voor het onderscheiden van gegevens van hoge en lage kwaliteit.

Figuur 1: Workflow van de RNA-seq in silico kwaliteitscontrolepijplijn. De invoerbestanden voor kwaliteitscontrole zijn afgeleid van RNA-seq-gegevens die zijn gegenereerd door het sequencen van plantmateriaal, evenals openbaar beschikbare datasets. De meegeleverde R-pijplijn beoordeelt de sequentiekwaliteit aan de hand van drie hoofdbenaderingen: sequencingkwaliteit, mappingkwaliteit en replicate-kwaliteit. Er worden verschillende kwaliteitsstatistieken berekend en statistische resultaten worden gevisualiseerd met behulp van veelgebruikte R-pakketten zoals ggplot2 en pheatmap (bijv. staafdiagrammen en heatmaps). Klik hier om een grotere versie van deze figuur te bekijken.

Eerder gerapporteerde evaluatiepijplijnen vereisten voorbewerkte gegevens (bijv. uitlijning lezen als .bam-bestanden), dekken niet alle metrische gegevens of worden niet langer onderhouden ^4,5,6. Het voordeel van de hier gepresenteerde workflow ligt in de volledigheid ervan door de meest voorkomende kwaliteitscontroleproblemen in één pijplijn te consolideren. Verder geven we voorbeelden van gegevens van hoge kwaliteit die geschikt zijn voor alle downstream-analysetoepassingen, maar ook voorbeelden van gegevens van lage kwaliteit, en bespreken we hun specifieke beperkingen voor verdere analyse. Verschillende eerder gepubliceerde rapporten beschrijven het doel van RNA-seq-analysetools en bieden vergelijkende evaluaties van hun prestaties ^7,8. De normen voor RNA-seq-kwaliteitscontrole en methodologierapportage zijn echter niet gestandaardiseerd en verslechteren de reproduceerbaarheid en biologisch zinvolle interpretaties van transcriptomics-experimenten.

Rup integreert standaard hoogwaardige tools, kwaliteitscontrole-analyse en visualisatie van de resultaten. Rup vereist ruwe sequencing-lezingen en een geannoteerd genoom als invoer en draait op Mac OS, Linux en "Windows Subsystem for Linux" (WSL) op Windows-systemen. Het integreert tools voor de kwaliteit van de sequencing en kwantificeert read mapping naar een genoom, omvat rRNA-inhoudsmeting in de monsters en biedt monstercorrelatie om replicaatcorrelatie te evalueren. De testgegevens werden verkregen met behulp van lasermicrodissectie van het centrale meristeemweefsel van twee verschillende stadia van de plantensoort Eschscholzia californica, een protocol met ultralage invoer voor bibliotheekbereiding, en gesequenced op Novaseq 6000.

Protocol

Rup kan worden uitgevoerd als een enkel script met een minimum aan vereiste invoerbestanden. De pijplijn, zoals hier weergegeven, vereist Illumina-gepaarde RNA-sequencinggegevens. Een aangepaste pijplijn voor single-end sequencing met dezelfde analysestappen wordt ook gedeponeerd in de GitHub-opslagplaats. Alleen de genoomsequentie als een fasta-bestand, het genmodel en de rRNA-annotaties als gtf-bestanden, en de onbewerkte sequencing leest als fastq.gz-bestanden zijn nodig. Zorg ervoor dat de genoom- en annotatiebestanden worden verstrekt als genome.fa, annotation.gtf en rRNA.gtf in een map die is opgegeven in de variabele reference_folder. Wanneer alle opeenvolgende bestanden als .fq.gz bestanden zijn opgeslagen in de variabele read_file_folder, kan Rup als volgt worden uitgevoerd:

1. Voorbereidingen:

1. Installeer alle benodigde R-pakketten met behulp van Bioconductor.

   BiocManager::install(c("getopt", "ggplot2", "reshape2", "pheatmap", "fastqcr", "Rfastp", "Rsubread", "Rsamtools"))

2. Stel alle benodigde bestanden in.
   1. Maak een bronmap die alle benodigde bestanden bevat. Voeg de referentiegenoomsequentie in fasta bestand met de naam "reference/genome.fa" toe aan de bronmap. Geef de annotatie van het genmodel op als een gtf-bestand: "reference/annotation.gtf", Geef optioneel de annotatie van het rRNA-gen op als gtf-bestand: "reference/rRNA.gtf".
   2. Voeg alle sequencing reads als fastq-bestanden toe aan de map "reads". Zorg ervoor dat alle fastq-bestandsnamen hetzelfde patroon volgen; _1.fastq.gz en _2.fastq.gz voor de vooruit- en achteruitleesbewerkingen.

      # define source folder and all derived subfolders for input and output files
source_folder <- "data"
reference_folder <- file.path(source_folder, "reference")
read_file_folder <- file.path(source_folder, "reads")
results_folder <- file.path(source_folder, "results")

# define input reference files
genome_fasta_file <- file.path(reference_folder, "genome.fa")
annotation_file <- file.path(reference_folder, "annotation.gtf")
rrna_file <- file.path(reference_folder, "rRNA.gtf")

# define result subfolders
fastqc_folder <- file.path(results_folder, "fastqc")
trimmed_fastqc_folder <- file.path(results_folder, "trimmed_fastqc")
trimmed_read_folder <- file.path(results_folder, "trimmed")
bam_folder <- file.path(results_folder, "bam")
sorted_bam_folder <- file.path(results_folder, "sorted_bam")
counts_folder <- file.path(results_folder, "counts")

# create results folders
for(folder in c(results_folder, fastqc_folder, trimmed_fastqc_folder, trimmed_read_folder, bam_folder, sorted_bam_folder, counts_folder)) {
 if(!dir.exists(folder)) {
  dir.create(folder)
 }
}

# get sample prefixes from input fastq files
fastq_files <- list.files(read_file_folder, pattern = "*_1.f(ast)?q.gz")
sample_prefixes <- gsub("_1.f(ast)?q.gz", "", fastq_files)

3. Stel de parameter in op basis van de rekenbronnen.

   n_threads <- 8 # the number of available CPU cores
bamSortMemory <- "1024" # maximum memory for bam file sorting

4. Stel de parameter in volgens de volgordemethode.
    

   MinReadLength <- 25 # minimum read length, should not be less than 25
minFragLength <- 0 # minimum fragment length distribution
maxFragLength <- 300 # maximum fragment length distribution
orientation <- "fr" # read orientation for read mapping ("fr", "rf", "ff")
stranded <- 0 # stranded sequencing
         # (0 (unstranded), 1 (stranded) and 2 (reversely stranded))

   
   
   OPMERKING: De minimale leeslengte mag niet kleiner zijn dan 25 bp, omdat kleinere leestijden ertoe kunnen leiden dat de mapping wordt verbroken. Grotere waarden resulteren in unieker in kaart gebrachte lezingen, maar ook in meer weggegooide lezingen tijdens het bijsnijden, afhankelijk van de volgordekwaliteit. Meer details over deze parameters worden beschreven in de supplementen.

2. Kwaliteitsbeoordeling in volgorde

OPMERKING: Met deze stap wordt een staafdiagram gemaakt met het aantal lezingen voor en na het bijsnijden van elk voorbeeld.

1. Voer fastqc⁹ uit, een hulpmiddel voor algemene kwaliteitscontrole van sequentiegegevens met hoge doorvoer, dat informatie biedt over het totale aantal gesequente lezingen, de leeslengte, de gemiddelde GC-inhoud, vervuiling van de sequentieadapter en oververtegenwoordigde sequenties (bijv. rRNA-lezingen). Analyseer de sequentiekwaliteit volgens basis- en sequentiekwaliteitsrapporten. Voer fastqc uit op de raw sequencing bestanden met behulp van het R package fastqcr¹⁰ en plaats alle invoerbestanden in dezelfde directory (read_file_folder). Vat de uitvoerbestanden van fastqc_folder samen in een R-object met behulp van de functie qc_aggregate.
   OPMERKING: De resultaten worden opgeslagen in de map fastqc_folder. De .html bestanden die door fastqc in de uitvoermap zijn gemaakt, zijn toegankelijk voor aanvullende statistieken.
    

   # load the fastqc library
library(fastqcr)
# run fastqc on all raw fastq files
fastqc(fq.dir=read_file_folder, qc.dir=fastqc_folder, threads=n_threads)

# aggregate the fastqc statistics
qc <- qc_aggregate(fastqc_folder, progress=F)
qc$tot.seq <- as.numeric(qc$tot.seq)
# remove suffixes from sample names
qc$sample = gsub(".f(ast)?q.gz", "" , qc$sample)

2. Met alle voorbeeldvoorvoegsels die zijn opgeslagen in de vector sample_prefixes, herhaalt u de voorvoegsels en voert u de kwaliteitsbijsnijdtool fastp uit op alle voorbeelden. Trimmen om adaptersequenties, primers, basissen van lage kwaliteit, poly-A/T-sequenties te verwijderen en wordt uitgevoerd door Fastp¹¹, dat beschikbaar is in het R-pakket Rfastp¹². De bijgesneden fastq-bestanden worden opgeslagen in de map trimmed_read_folder.
    

   # load the rfastp library
library(Rfastp)

# iterate over sample prefixes
for(prefix in sample_prefixes) {
 # create output prefix
 # !!! Rfastp automatically adds _R1.fastq.gz and _R2.fastq.gz to the prefix
outputPrefix <- file.path(trimmed_read_folder, prefix)

# get input reads based on sample prefix
read1 = file.path(read_file_folder, paste(prefix, "_1.fastq.gz", sep=""))
read2 = file.path(read_file_folder, paste(prefix, "_2.fastq.gz", sep=""))
if(!file.exists(read1)) { read1 = file.path(read_file_folder, paste(prefix, "_1.fq.gz", sep="")) }
if(!file.exists(read2)) { read2 = file.path(read_file_folder, paste(prefix, "_2.fq.gz", sep="")) }

# run fastp trimmig with minimum read length
fastp_stats <- rfastp(read1 = read1,
     read2 = read2,
     minReadLength = minReadLength,
     outputFastq = outputPrefix,
     thread = n_threads)
}

3. Voer fastqc opnieuw uit op de bijgesneden reads.
    

   # run fastqc on the trimmed reads
fastqc(fq.dir=trimmed_read_folder, qc.dir=trimmed_fastqc_folder, threads=n_threads)

# aggreagate the fastqc statistics
qc_trimmed <- qc_aggregate(trimmed_fastqc_folder, progress=F)
qc_trimmed$tot.seq <- as.numeric(qc_trimmed$tot.seq)
# correct sample names (change the fastp "R1","R2" suffix to "1","2")
qc_trimmed$sample = gsub("_R([12])$", "_\\1" , qc_trimmed$sample)

4. Evalueer het bijsnijden door het aantal lezingen voor en na het bijsnijden voor alle monsters te verzamelen. Maak een staafdiagram met de gelezen getallen met behulp van ggplot2¹³.
    

   # load the ggplot2 library for plotting
library(ggplot2)

# add the trimming status to the fastqc results
qc$Trimming = "raw"
qc_trimmed$Trimming = "trimmed"

# combine the results into one vector
qc_all = rbind(qc, qc_trimmed)

# plot the read number before and after trimming as barplot
read_number_plot <- ggplot(qc_all, aes(x=sample, y=tot.seq, fill=Trimming)) +
geom_bar(stat="identity", position = "dodge") +
xlab("Samples") + ylab("Number of Reads") +
theme(axis.text.x = element_text(angle = 90, hjust = 0)) +
theme(text = element_text(size = 18)) +
ggtitle("Number of reads before and afer trimming")
print(read_number_plot)

Figuur 2: Resultaten van sequentiebepaling en bijsnijden. Leestellingen voor (rood) en na het bijsnijden (groen). Sample s2_r1 heeft al een laag aantal gelezen lezen voorafgaand aan het trimmen, terwijl bij het bijsnijden een groot deel van de sample s2_r2 lezingen werd verwijderd. Alle andere monsters vertonen een acceptabel verlies van aflezingen door bijsnijden. Klik hier om een grotere versie van deze figuur te bekijken.

3. Kwaliteit in kaart brengen

OPMERKING: Deze stap berekent staafdiagrammen om onderscheid te maken tussen enkelvoudig in kaart gebrachte lezingen (bijv. transcripten van eiwitcoderende genen) en niet-in kaart gebrachte lezingen (bijv. rRNA-lezingen) en niet-in kaart gebrachte lezingen (bijv. contaminaties).

1. Gebruik het Rsubread^14-pakket om de lezingen toe te wijzen aan het referentiegenoom. Bouw eerst een index van het referentiegenoom fasta-bestand (genome_fasta_file).
   OPMERKING: Deze stap wordt één keer uitgevoerd voor elk referentiegenoom.
    

   # load the Rsubread library for read mapping and read counting
library(Rsubread)

# create a subread index from the reference genome sequence
buildindex(file.path(reference_folder,"subread.index"), genome_fasta_file)

2. Herhaal alle samples en lijn de lezingen uit met het referentiegenoom met behulp van de align() -functie van het Rsubread-pakket. Met deze functie worden .bam-bestanden gemaakt in de uitvoermap.
    

   # load the Rsamtools library for bam file sorting and indexing
library(Rsamtools)

# iterate over sample names
for(prefix in sample_prefixes) {
# create the bam output file name
output_bam = file.path(bam_folder, paste(prefix, ".bam", sep=""))

# align the reads to the reference genome index
mapping_stats <- align(index=file.path(reference_folder,"subread.index"),
     # create the forward and reverse read file names based on the prefix
     # "_R1.fastq.gz" and "_R2.fastq.gz" are forced by fastp
     readfile1=file.path(trimmed_read_folder, paste(prefix, "_R1.fastq.gz", sep="")),
     readfile2=file.path(trimmed_read_folder, paste(prefix, "_R2.fastq.gz", sep="")),
     output_file=output_bam, # set the output file name
     type=0, # the reads are RNA-seq
     minFragLength=minFragLength, # the minimal allowed fragment length
     maxFragLength=maxFragLength, # the minimal allowed fragment length
     PE_orientation=orientation, # read orientation
     nthreads=n_threads,
     # use a GTF annotation file to support the mapping
     useAnnotation=TRUE,
     annot.ext=annotation_file,
     isGTF=TRUE,
     nBestLocations=2) # to distinguish between unique and multi-mapping reads

# create the sorted output file name (the .bam suffix will be added automatically)
output_sorted_bam = file.path(sorted_bam_folder, prefix)

# sort and index the bam file
sortBam(output_bam, output_sorted_bam, maxMemory=bamSortMemory, nThreads=n_threads)
indexBam(paste0(output_sorted_bam, ".bam"))
}

3. Tel de lezingen voor elk gen met behulp van de functie featureCounts() uit het Rsubread-pakket. Zorg ervoor dat het annotatiebestand (annotation_file) de GTF-indeling heeft. Count leest alleen met een enkele overeenkomst met het genoom.
    

   # collect all bam files
bam_files <- list.files(bam_folder, pattern = "*.bam$")

# count the number of reads for each gene
     gene_feature_counts = featureCounts(file.path(bam_folder, bam_files),
     annot.ext = annotation_file, # the gene models
     isGTFAnnotationFile = TRUE,
     countMultiMappingReads = FALSE, # only count unique reads
     strandSpecific = stranded, # for strand specific data
     isPairedEnd = TRUE, # paired-end sequencing
# the following paramter allow to control for false-positive read assignments
     requireBothEndsMapped = TRUE,
     checkFragLength = TRUE,
     minFragLength = minFragLength,
     maxFragLength = maxFragLength,
     nthreads = n_threads)

4. Tel de lezingen die in kaart worden gebracht op rRNA-genen voor het beoordelen van het rRNA-gehalte in de monsters. Sta toe dat multimapped reads hier worden meegeteld. Specificeer de rRNA-genloci in een GTF-bestand (rrna_file).
    

   # count the reads mapping to rRNA genes
rrna_feature_counts = featureCounts(file.path(bam_folder, bam_files),
     annot.ext = rrna_file, # the rRNA gene model file
     isGTFAnnotationFile = TRUE,
     countMultiMappingReads = TRUE, # rRNA reads usually map to multiple loci
     fraction = TRUE,
     strandSpecific = stranded,
     isPairedEnd = TRUE,
     nthreads = n_threads)

5. Verzamel de statistieken van leesopdrachten die zijn gemaakt door featureCounts. De statistieken omvatten het aantal uitlijningen voor elke steekproef in de verschillende categorieën (bijv. Assigned, UnMapped, MultiMapped).
    

   # load the reshape2 library for data restructuring
library(reshape2)

# get the assignement stats
gene_count_stats <- gene_feature_counts$stat
# set the sample names as column names
colnames(gene_feature_counts$counts) = gsub(".bam", "", colnames(gene_feature_counts$counts))

# set the assignment type as row names
rownames(gene_count_stats) <- gene_count_stats$Status
# select the columns with the assignment statistics
gene_count_stats <- gene_count_stats[,seq(2, ncol(gene_count_stats))]
# set the sample names as column names
colnames(gene_count_stats) <- gsub(".bam", "", colnames(gene_count_stats))

# transform the dataframe for plotting
transformed_stats <- melt(t(gene_count_stats))

# set the column names
colnames(transformed_stats) <- c("Sample", "Group", "Alignments")

# adjust the groupa and sample ordering for plotting
transformed_stats$Group <- factor(transformed_stats$Group,
          levels = rev(levels(transformed_stats$Group)[order(levels(transformed_stats$Group))]))
transformed_stats$Sample <- factor(transformed_stats$Sample,
          levels = rev(levels(transformed_stats$Sample)[order(as.character(transformed_stats$Sample))]))

# remove assignment classes without any alignments
transformed_stats <- transformed_stats[transformed_stats$Alignments > 0,]
# the statistics refer to proteins coding genes
transformed_stats$Reference = "Genes"

6. Verzamel de statistieken van de rRNA-toewijzing.
    

   # collect the rrna assignment statistics
rrna_count_stats <- rrna_feature_counts$stat
# set the sample names as column names
colnames(rrna_count_stats) <- gsub(".bam", "", colnames(rrna_count_stats))

# only the number of assigned reads are important in this case
rrna_counts = as.data.frame(t(rrna_count_stats[rrna_count_stats$Status == "Assigned",2:ncol(rrna_count_stats)]))
colnames(rrna_counts) = "Alignments"

# setup names and categories for plotting
rrna_counts$Sample = rownames(rrna_counts)
rrna_counts$Group = "rRNA"
rrna_counts$Reference = "rRNA"

7. Teken de gelezen toewijzingsstatistieken als een staafdiagram.
    

   # combine the proteind coding and rRNA gene statistics
mapping_stats = rbind(transformed_stats, rrna_counts)

# plot the assignment statistics
mapping_stats_plot = ggplot(data = mapping_stats, aes(x = Sample, y = Alignments)) +
 geom_col(aes(fill = Group), width = 0.7) +
 theme_bw() + facet_wrap(~Reference) +
 ylab("Number of Alignments") + xlab("Samples") +
 ggtitle("Read Mapping Numbers") +
 theme(text = element_text(size = 18)) +
 theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(mapping_stats_plot)

8. Classificeer genen in groepen op basis van het aantal lezingen dat aan hen is toegewezen en plot de resultaten als een staafdiagram.
    

   # get the reads per gene counts
gene_count_matrix = gene_feature_counts$counts
# set the sample names as column names
colnames(gene_count_matrix) = gsub(".bam", "", colnames(gene_count_matrix))

# group and count genes in classes of specific read counts
read_count_classes = data.frame("no_reads"=colSums(gene_count_matrix == 0),
     "at_least_1_read"=colSums(gene_count_matrix >= 1 & gene_count_matrix < 10),
     "at_least_10_reads"=colSums(gene_count_matrix >= 10 & gene_count_matrix < 100),
     "at_least_100_reads"=colSums(gene_count_matrix >= 100 & gene_count_matrix < 1000),
     "at_least_1000_reads"=colSums(gene_count_matrix >= 1000))

# setup data.frame for plotting
read_count_classes$Sample = rownames(read_count_classes)
melt_rcc = melt(read_count_classes)
melt_rcc$variable = as.character(melt_rcc$variable)

# sort the gene groups
melt_rcc$variable = factor(melt_rcc$variable, levels=c("no_reads",
         "at_least_1_read",
         "at_least_10_reads",
         "at_least_100_reads",
         "at_least_1000_reads"))

# plot the data as bar plot
gene_coverage_plot <- ggplot(melt_rcc, aes(x=Sample, y=value, fill=variable)) +
geom_bar(stat="identity") +
ylab("Number of Genes") + xlab("Samples") +
ggtitle("Number of Reads per Gene") +
guides(fill=guide_legend(title="Number of assigned reads")) +
theme(text = element_text(size = 18)) +
theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(gene_coverage_plot)

Figuur 3: Lees het kaartoverzicht en de statistieken. De meeste monsters vertonen een hoog aantal toegewezen lezingen (linkerkant, bruin) en een laag aantal rRNA-lezingen (rechterkant, blauw). Steekproef s2_r3 heeft een ongewoon hoog aantal multi-mapped lezingen (linkerkant, groen) die overeenkomen met een hoog rRNA-leesnummer (rechterkant). Steekproef s2_r4 toont een groot aantal niet-in kaart gebrachte lezingen in combinatie met een verwacht rRNA-afleesnummer, wat wijst op besmetting met aflezingen van een ander organisme. Klik hier om een grotere versie van deze figuur te bekijken.

4. Repliceer kwaliteit.

1. Zorg ervoor dat duplo's van hetzelfde monster of weefsel een hogere correlatie met elkaar vertonen dan met duplo's van verschillende monsters of weefsels. Teken een geclusterde heatmap van de gerepliceerde correlaties om de gegevens te visualiseren met behulp van de onbewerkte leestellingen of de genormaliseerde tellingen per miljoen (TPM)¹⁵ .
    

   # load the pheatmap library
library(pheatmap)

# compute TPM values
Length_kb = gene_feature_counts$annotation$Length / 1000 # get gene lengths in kb
RPK = gene_feature_counts$counts / Length_kb # normalize read counts by gene length
scaling_factors = colSums(RPK) / 1e6 # get scaling factor based on the sum normalized read counts
TPM = RPK / scaling_factors # scale the normalized read counts to 1e6

# plot the sample/replicate correlation heatmap
# the pearson correlation is computed on the log2 transformed TPM values
pheatmap(cor(log2(TPM+1)), fontsize = 18, main="Sample Correlation Heatmap")

Figuur 4: Classificatie van het aantal gelezen genen. Alle genen in een monster worden geclassificeerd in een van de vijf categorieën op basis van het aantal toegewezen lezingen. In rood wordt het aantal genen weergegeven waaraan geen enkele reading is toegewezen. Monsters met een laag aantal totale toegewezen lezingen (s2_r1 tot s1_r4) hebben een hoger aantal genen met 10 tot 100 toegewezen lezingen en een lager aantal genen met meer dan 1000 lezingen. Genen met een lage expressie kunnen in deze monsters worden gemist. Klik hier om een grotere versie van deze figuur te bekijken.

Resultaten

De pipeline is volledig geïmplementeerd als een R-script en is getest op Linux- en Mac OS-besturingssystemen. Windows-gebruikers kunnen Windows Subsystem for Linux (WSL) gebruiken. De code en testgegevens zijn beschikbaar als een GitHub-opslagplaats: https://github.com/oliverrupp/rup. De sequentiegegevens zijn beschikbaar in het kader van het ENA EBI-project PRJEB96400.

Tien monsters werden kunstmatig gemaakt van twee echte monsters om diverse problemen te illustreren die zich kunnen voordoen tijdens kwaliteitscontrole van bulk-RNA-seq met behulp van Rup voor analyse. Voorbeeld s2_r1 is ontworpen om een laag totaal leesgetal te vertonen, monster s2_r2 bevat een groot deel van de lezingen van lage kwaliteit die door het bijsnijdproces moeten worden weggegooid. Monster s2_r3 bevat een groot deel van de rRNA-lezingen en monster s2_r4 bevat contaminante lezingen die niet in kaart konden worden gebracht naar het referentiegenoom. De namen van de steekproeven s1_r5 en s2_r5 werden verwisseld om een lage replicaatcorrelatie te illustreren.

Protocolsectie 2 maakt de identificatie mogelijk van monsters met lage leesnummers voor of na het bijsnijden (Figuur 2). De staafdiagram toont het onderste leesnummer in het voorbeeld s1_r1 zowel voor als na het bijsnijden. Hier was het aanvankelijke leesnummer laag. Het lage leesaantal van het monster dat na het bijsnijden s1_r2, suggereert een grote hoeveelheid adaptersequenties, sequentiefouten, primersequenties, poly-A/T-rekken die tijdens het trimproces zijn verwijderd. Degradatie van invoer-RNA kan ook het aantal lezingen van hoge kwaliteit verminderen.

Protocol sectie 3 identificeert problemen in leesopdrachten. Figuur 3 toont het verhoogde aantal multi-mapped lezingen in steekproef s2_r3 (groen), evenals het hoge aantal rRNA-lezingen. Contaminatie van het s2_r4 van het monster blijkt uit het grote deel van de lezingen dat niet is toegewezen aan het referentiegenoom (roze). Deze lezingen correleren niet met rRNA, maar met sequenties van een niet-doelorganisme. Figuur 4 toont algemene problemen die verband houden met een laag aantal toegewezen lezingen in transcriptomen. In monsters met een laag aantal reads die uniek zijn toegewezen aan het referentiegenoom (s2_r1 tot r4), heeft slechts ongeveer een derde van de genen meer dan 100 reads toegewezen, terwijl in de andere monsters ongeveer de helft van de genen in deze klassen valt. Het is mogelijk dat lezingen van genen met een zeer lage expressie niet worden gevonden in de monsters s2_r1 tot r4, en daarom zal vergelijkende expressieanalyse met deze monsters zeer onbetrouwbaar zijn en moet worden vermeden.

Protocolsectie 4 kan worden gebruikt om uitschieters van replica's te identificeren. Van replicaten van hetzelfde monster/dezelfde aandoening/weefsel wordt verwacht dat ze een hogere correlatie tussen elkaar vertonen dan replicaten van andere monsters/aandoeningen/weefsels. Figuur 5 toont een correlatie-heatmap van de twee monsters (S1 en S2) met elk vijf replicaten. De dendrogrammen aan de boven- en zijkant van de grafiek tonen twee clusters met vijf replicaten in elke cluster. Het linkercluster bevat vier replicaten van monster 1 en één replicaat van monster 2 (s2_r5), het rechtercluster bevat vier replicaten van monster 2 en één replicaat van monster 1 (s1_r5). In dit geval, wanneer de voorbeeldnamen s1_r5 en s2_r5 weer worden verwisseld, bevat elk cluster alle replicaten van één voorbeeld, wat kan duiden op een replicatie-labelfout. Andere redenen waarom replicaten niet samenklonteren, kunnen een gebrek aan differentiatie tussen de monsters/aandoeningen/weefsels zijn, of clustering van replicaten alleen om redenen van sequentiekwaliteit. Dit laatste kan optreden wanneer alle replicaten met een uitzonderlijk laag aantal of die met een uitzonderlijk hoog aantal leesbewerkingen na het bijsnijden en in kaart brengen een cluster vormen.

Figuur 5: Voorbeeld correlatie heatmap. De heatmap voor voorbeeldcorrelatie is gebaseerd op log₂ getransformeerde TPM-waarden en toont twee afzonderlijke clusters van elk vijf monsters. Van biologische/technische replicaten wordt verwacht dat ze een hogere correlatie tussen elkaar vertonen dan die van andere weefsels/behandelingen. Het linkercluster bevat vier replicaten van monster 1 en één replicaat van monster 2 (s2_r5), het rechtercluster bevat vier replicaten van monster 2 en één replicaat van monster 1 (s1_r5). De afzonderlijke replicaten moeten worden gecontroleerd op mogelijke sample-swapping of batch-effecten om hun clustering in de heatmap te verklaren. Klik hier om een grotere versie van deze figuur te bekijken.

Tabel 1: Vergelijking van RNA-seq-kwaliteitsbeoordelingspijplijnen. Voor een gedetailleerde analyse, zie Aanvullend Dossier 1. Klik hier om deze tabel te downloaden.

Aanvullend bestand 1: Selectie van parameters en referenties. De analyses tonen de invloed van parameter- en referentieselectie op de algemene QC-resultaten. Klik hier om dit bestand te downloaden.

Discussie

De kwaliteit van differentiële genexpressieanalyses hangt sterk af van twee factoren: het aantal lezingen dat in elk monster en elke replicatie wordt gesequenced¹⁶ en het aantal replicaten per monster^17,18. Hier presenteren we de gebruiksvriendelijke Rup-pijplijn om het aantal lezingen te bepalen dat geschikt is voor kwantificering van genexpressie in elke replicaat. Verschillende statistieken stellen onderzoekers in staat om te begrijpen waarom replicaten lage toegewezen leesnummers vertonen en om problemen in de correlatie van monsters te identificeren. Hoewel Rup is ontwikkeld voor kwaliteitsbeoordeling van RNA-seq-monsters van planten, is het even geschikt voor andere eukaryote organismen, waarvoor geen verdere aanpassingen nodig zijn voor niet-plantmonsters (aanvullend bestand 1).

De eerste stap van Rup bepaalt het totale aantal gelezen lezen voor en na het bijsnijden van het lezen. Het aantal gesequente lezingen bepaalt het aantal detecteerbare genen, en als het afleesaantal te laag is, zullen veel differentieel tot expressie gebrachte genen onopgemerkt blijven. Een groot deel van de lezingen die tijdens het kwaliteitstrim- en filterproces worden weggegooid, kan duiden op RNA-degradatie, wat zou kunnen leiden tot een onderschatting van de expressie voor sterk gedegradeerde genen. Of een monster kan worden gebruikt voor downstream-toepassingen hangt echter af van het doelorganisme en het doel van het onderzoek. Om bijvoorbeeld de expressie van de meeste plantengenen vast te leggen, zijn naar onze ervaring 30 miljoen tot 50 miljoen reads nodig, terwijl voor schimmels slechts 10 miljoen voldoende zou kunnen zijn. Rup zal dus geen kwaliteitsdrempels definiëren voor de uitsluiting van problematische steekproeven, maar eerder maatstaven bieden om diverse problemen te identificeren die zich kunnen voordoen.

De tweede stap (karteringskwaliteit) evalueert de nauwkeurigheid en betrouwbaarheid van de uitlijning van het lezen met het referentiegenoom en specificeert de geschiktheid ervan voor stroomafwaartse analyses. Niet alle gesequente lezingen kunnen worden gebruikt voor de berekening van genovervloed; Leesbewerkingen die niet overeenkomen met het genoom of transcriptoom, of lezingen die worden toegewezen aan meerdere locaties op het genoom, worden in de meeste gevallen genegeerd¹⁹ en dragen niet bij aan het totale aantal gelezen lezen. De resultaten van de tweede stap van de pijplijn kunnen worden gebruikt om te begrijpen waarom lezingen niet worden gebruikt in de berekening van de abundantie. Een groot aantal niet-in kaart gebrachte lezingen kan duiden op besmetting tijdens RNA-extractie (bijv. met plantpathogenen of herbivoren) of een onvolledig referentiegenoom. Onvolledige genmodellen van het referentiegenoom kunnen resulteren in een groot aantal "no feature" reads. Een groot deel van de multi-mapped reads kan worden veroorzaakt door een groot aantal rRNA-genen in de bibliotheek, wat wijst op onvoldoende rRNA-verwijdering tijdens het voorbereidingsproces van de sequencingbibliotheek (in het geval dat de multi-mapped reads echt zijn afgeleid van rRNA's, kunnen ze verder worden geanalyseerd door een rRNA-annotatiebestand aan de pijplijn te verstrekken, die vervolgens automatisch het aantal mogelijke rRNA-lezingen berekent).

De derde, even belangrijke kwaliteitsmaatstaf is de correlatie tussen replicaten van hetzelfde monster. Over het algemeen zou de correlatie tussen replicaten van hetzelfde monster hoger moeten zijn dan de correlatie tussen replicaten van verschillende monsters. Een lage correlatie tussen replicaten kan duiden op grote biologische variabiliteit tussen de replicaten, hoge gelijkenis tussen de monsters/omstandigheden/weefsels, andere batcheffecten, of zelfs monsterruil of verkeerde etikettering. Rup berekent paarsgewijze correlaties tussen alle steekproeven en produceert een geclusterde heatmap van de steekproefcorrelaties. Bovendien kan een hoofdcomponentenanalyse (PCA) op de monsters worden berekend om mogelijke batcheffecten te identificeren die in verdere stroomafwaartse analyses moeten worden erkend. Hoewel batcheffecten kunnen worden gecorrigeerd in downstream-analyses, is het misschien beter om problematische monsters te verwijderen als het gemeten verschil te groot is.

Invoermateriaal heeft een directe invloed op de kwaliteit van de sequencing: weefselbemonstering, RNA-extractie en bibliotheekvoorbereiding zijn cruciale stappen om problemen met de RNA-seq-kwaliteit tot een minimum te beperken. Indien mogelijk moeten consistente omstandigheden worden gehandhaafd, bijvoorbeeld door gebruik te maken van groeikamers. Ongediertebestrijdingsmaatregelen moeten tijdig worden uitgevoerd, aangezien alleen gezonde personen voor bemonstering mogen worden geselecteerd. Het is verder raadzaam om monsters op dezelfde dag en/of op hetzelfde tijdstip te nemen om de variatie in het circadiane transcriptoom te verminderen. Er zijn talloze RNA-extractiekits beschikbaar en het selecteren van een geschikte kit voor de doelsoort kan de mRNA-kwaliteit verbeteren. De voorbereidingsprotocollen van de bibliotheek moeten stappen bevatten voor polyA+ RNA-verrijking om de rRNA-fractie te minimaliseren.

Rup kan worden gebruikt als een eerste kwaliteitscontrolestap in de analyse van differentiële genexpressie. Deze kwaliteitscontrole is om verschillende kritieke redenen vereist, omdat het de kwaliteit van de invoergegevens waarborgt (replicaten/samples van lage kwaliteit identificeert, kwaliteitsdrempels voor leesdiepte en toewijzingssnelheden kan instellen) en technische problemen identificeert, zoals volgordefouten, batcheffecten of verkeerde etikettering. Als de RNA-invoer van voldoende kwaliteit is, kan Rup helpen door lezingen van lage kwaliteit bij te snijden of verkeerd gelabelde monsters te identificeren. Rup kan echter geen slechte input-RNA-kwaliteit compenseren. In het geval van RNA van lage kwaliteit of foutieve sequencinggegevens, kan het nodig zijn om monsters opnieuw te verzamelen, het RNA-extractieprotocol aan te passen en/of de sequencing te herhalen. Hetzelfde geldt voor monsters met een groot deel van niet-in kaart gebrachte lezingen die mogelijk zijn veroorzaakt door vervuiling. Uit onze ervaring blijkt echter dat niet alle RNA-extracties zomaar kunnen worden herhaald vanwege een gebrek aan materiaalbeschikbaarheid. In dit geval zijn ze ongeschikt voor sommige downstream-toepassingen: monsters met een laag aantal single-mapped reads mogen niet worden geanalyseerd in differentiële genexpressieanalyses, maar ze bevatten nog steeds informatie voor de analyse van de aanwezigheid van transcripten. In dit geval kan de afwezigheid van transcripten in weefsels/behandelingen/aandoeningen niet in aanmerking worden genomen voor de analyse en kwantificering van de overvloed aan transcripten.

Rup bevat enkele beperkingen. Het test bijvoorbeeld niet direct op RNA-degradatie, omdat directe metingen van de RNA-integriteit vereist zijn vóór de voorbereiding en sequencing van de bibliotheek. Een script om RNA-degradatie te identificeren met behulp van de RSeQC-modules geneBodyCoverage.py en tin.py is echter opgenomen in de repository van deze pijplijn. Bovendien test Rup niet op GC-inhoud en transcriptielengtebias. De referentiegenoomgrootte is momenteel beperkt tot 4 Gb, en volgens onze ervaring overschat de read mapping-module multi-mapped reads in polyploïden, zoals geïllustreerd in het aanvullende bestand in een vergelijking van read mapping met een haploïde versus de diploïde genoomversies van E. californica, een haploïde genoom is dus de voorkeursinput voor deze pijplijn. De kwaliteit van de Rup-output hangt grotendeels af van de kwaliteit van het referentiegenoom en de annotatie, zodat een onvolledige of sterk gefragmenteerde genoomsequentie kan leiden tot een overschatting van niet-in kaart gebrachte lezingen. Bovendien leidt onvolledige annotatie van genmodellen tot een overschatting van niet-toegewezen lezingen. Volledigheid en duplicatie van de genoomsequentie en de genannotatie kunnen worden afgeleid met tools zoals BUSCO²⁰.

Rup is een stand-alone tool voor alle initiële kwaliteitscontrolestappen die essentieel zijn in RNA-seq-experimenten. In tegenstelling tot RSeQC en RNA-SeQC is voorverwerking van de onbewerkte sequencing-lezingen niet vereist. Kwaliteitsanalyse van sequencing kan niet worden uitgevoerd met RNA-SeQC en heatmaps voor monstercorrelatie worden niet berekend door RSeQC en RNA-SeQC (Tabel 1). Verder is de normalisatie-uitvoer in FPKM in RSeQC en is uitvoervisualisatie niet geïmplementeerd als standaard voor alle modules in RSeQC en RNA-SeQC. Verschillende kwaliteitscontroleproblemen worden dus niet getest in de twee alternatieve tools, waaronder een laag leesnummer, een hoog deel van de bijgesneden leeswaarden en de identificatie van replicate uitschieters. Een vergelijking van Rup, RSeQC en RNA-SeQC is beschikbaar in aanvullend bestand 1. Bovendien kan Rup worden gebruikt als aanvulling op RNA-SeQC of RSeQC, bijvoorbeeld de gesorteerde BAM-bestanden die door deze pijplijn worden geproduceerd, kunnen worden gebruikt als input voor deze tools.

Momenteel is Rup geoptimaliseerd voor een klein aantal samples, maar de meest tijdrovende stappen, zoals quality trimming en read mapping, kunnen vooraf worden berekend, bijvoorbeeld op een computercluster of een cloudinfrastructuur. Voor genomen groter dan 4 Gb is dit verplicht, aangezien de Rsubread-aligner beperkt is tot genomen kleiner dan 4 Gb. rRNA-annotatie wordt gedaan met barrnap, die de sterk geconserveerde rRNA-genen identificeert. Volgens onze ervaring vereist annotatie van rRNA-gen geen volledigheid, want zelfs als sommige abnormale rRNA-genen ontbreken, is een bruto overzicht van de aanwezigheid van rRNA in de dataset voldoende voor de kwaliteitsbeoordeling van de dataset. Toekomstige versies van Rup kunnen overschakelen naar een andere mapping-methode, zoals de pseudo-uitlijningstool salmon²¹, om de runtime te verkorten. Bovendien zouden meer normalisatie- en correctiemethoden, zoals GC-bias of lengtebiascorrectie, aan Rup kunnen worden toegevoegd.

Samenvattend biedt Rup essentiële informatie om de betrouwbaarheid en reproduceerbaarheid van de RNA-seq-gegevensanalyse te waarborgen. Deze pijplijn rapporteert uitgebreid de belangrijkste RNA-seq-kwaliteitsstatistieken en produceert uitvoerbestanden voor direct gebruik in downstream-analyses. Het is ontworpen als een op zichzelf staand hulpmiddel voor onderzoekers met minimale kennis van bio-informatica om de kwaliteit van hun primaire sequencinggegevens te beoordelen met intuïtieve visualisatie.

Materialen

Lijst van materialen gebruikt in dit artikel

Naam	Bedrijf	Catalogusnummer	Opmerkingen
fastqcr	R	0.1.3
FUJITSU ESPRIMO D958 Desktop	FUJITSU		de pipeline werd ontwikkeld en getest op Ubuntu Linux 24.02 (32 Gb RAM)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	Apple		de pipeline werd getest op maxOS 15.4.1 (16 Gb RAM)
pheatmap	R	1.0.13
R		4.4.3
reshape2	R	1.4.4
rfastp	bioconductor	1.16.0
rsamtools	bioconductor	2.22.0
rsubread	bioconductor	2.20.0