Dieses Protokoll ermöglicht eine erste Qualitätskontrolle für RNA-seq-Experimente für Nasslaborbiologen mit begrenzter Erfahrung in der Bioinformatik.
Methodenartikel
Dieses Protokoll ermöglicht eine erste Qualitätskontrolle für RNA-seq-Experimente für Nasslaborbiologen mit begrenzter Erfahrung in der Bioinformatik.
Moderne Ansätze in der molekularen Pflanzenwissenschaft erfordern häufig RNA-seq-Experimente, z. B. um globale Veränderungen in Transkriptomen nach Behandlungen zu verfolgen oder Schlüsselkomponenten von Regulationswegen zu identifizieren. Folglich sind verschiedene Bereiche der Pflanzenwissenschaften auf qualitativ hochwertige und reproduzierbare RNA-seq-Massendaten angewiesen, um den wissenschaftlichen Fortschritt voranzutreiben. Aus unserer Erfahrung fehlen jedoch häufig Kenntnisse über und die Anwendung von Qualitätskontrollmaßnahmen in RNA-seq-Datensätzen. Hier stellen wir Rup (RNA-seq usability assessment pipeline) für die Qualitätskontrolle von RNA-seq-Massendaten für nachfolgende Genexpressionsanalysen vor, das eigenständig und für Nasslaborbiologen mit grundlegenden R-Kenntnissen leicht anwendbar ist. Rup hilft bei der Unterscheidung zwischen qualitativ hochwertigen Sequenzierungsdaten, die für nachgelagerte Genexpressionsexperimente geeignet sind, und solchen, die für eine allgemeine weitere Analyse ungeeignet sind. Rup umfasst Tests für mehrere häufig auftretende Probleme, wie z. B. unzureichende Lesezahlen oder Mapping, Identifizierung von Kontaminationen, Quantifizierung von rRNA-Fraktionen in den gesamten RNA-seq-Daten, Replikationsähnlichkeitstests, Verwendung von Echtdaten zur Demonstration und intuitive Visualisierung. Rup bietet eine Reihe von Werkzeugen, um experimentelle Mängel vor der standardisierten Transkriptomanalyse zu identifizieren und so die Datenqualität für einzelne Forscher und das Feld zu verbessern. Dies erhöht das Vertrauen in die Analyse von RNA-Seq-Massendaten und bildet eine Grundlage für zukünftige Richtlinien, die Mindestkriterien für die Qualitätskontrolle definieren, wodurch die Zuverlässigkeit und Transparenz der veröffentlichten RNA-Seq-Daten verbessert wird. Rup- und Testdaten sind unter https://github.com/oliverrupp/rup verfügbar.
Transkriptomik-Experimente (RNA-seq) untersuchen umfassend transkriptionelle Signaturen, die Phänotypen formen. Dieser Ansatz wurde in der Molekulargenetik von Pflanzen unersetzlich, um einzelne oder ko-exprimierte Gene und biologische Prozesse zu identifizieren, die beispielsweise an der Entwicklung, der Interaktion mit Pflanzenpathogenen oder der abiotischen Stressresistenz beteiligt sind1, 2 ,3. Jüngste Fortschritte in der RNA-seq-Technologie haben die Spezifität erhöht und den Nachweis verschiedener Isoformen und Varianten mit Einzelbasenauflösung ermöglicht, was die Identifizierung von Sequenzvariationen von größeren Indels bis hin zu Einzelnukleotid-Polymorphismen (SNPs) ermöglicht. Die durch RNA-seq erhaltenen Daten zeichnen sich durch einen breiten dynamischen Bereich aus, der den Nachweis sowohl von stark häufig als auch von niedrig exprimierten Transkripten erleichtert, und erfordern geeignete experimentelle Bedingungen für die Konsistenz. Darüber hinaus sind RNA-seq-Datensätze groß und daher rechenintensiv zu analysieren und zu speichern, was ein effizientes Datenmanagement und umfangreiche Rechenressourcen erfordert. RNA-seq erfordert bei jedem Schritt des Arbeitsablaufs einen qualitativ hochwertigen Input, da sich Schwachstellen in jeder Phase ausbreiten und die Ergebnisse beeinträchtigen können. Schlechte RNA-Integrität oder technische Probleme bei der Bibliotheksvorbereitung führen zu Verzerrungen und einer Verringerung der Genauigkeit. Die Parameter für die Erzielung hochwertiger Roh-Reads können je nach Next-Generation-Sequencing-Anlage (NGS) variieren.
Nach unseren Erfahrungen empfehlen wir, nur RNA mit einer RNA-Integritätszahl (RIN) über 7 als Input für die RNA-seq-Bibliothekserstellung zu verwenden, die auf eine weitgehend intakte mRNA-Struktur hinweist. Für eine erfolgreiche Sequenzierung sind etwa 2 μg Gesamt-RNA in einer Konzentration von 50–200 ng/μl für Standard-Bibliotheksvorbereitungsprotokolle erforderlich. Die RNA-Reinheit sollte durch ein OD260/280-Verhältnis zwischen 1,8 und 2,1 und ein OD260/230-Verhältnis von mehr als 1,5 mit einem Spektralphotometer bestätigt werden. Unzureichende Sequenzierungstiefe, niedrige Kartierungsraten oder eine Fehlausrichtung auf das Referenzgenom können die Genexpression und die Spleißquantifizierung weiter verzerren. Darüber hinaus kann ein unzureichendes Versuchsdesign, wie z. B. mangelnde Unterscheidung zwischen Proben verschiedener Gewebe oder Behandlungen und eine hohe Variabilität zwischen den Replikaten, zu Rauschen führen und die Reproduzierbarkeit verringern. Während viele Labore Transkriptome routinemäßig analysieren, werden strenge Qualitätskontrollen der ersten wesentlichen Analyseschritte oft nicht berichtet oder fehlen in Publikationen. Dies kann zu einer Überinterpretation der aus der Transkriptomanalyse abgeleiteten Ergebnisse und damit zu nicht reproduzierbaren Ergebnissen führen.
Hier bieten wir einen Workflow für die Qualitätskontrolle der ersten Schritte, die für qualitativ hochwertige Transkriptomanalysen von mRNA-Profilen zur Messung von transkriptionellen Veränderungen erforderlich sind. Unser Ziel ist es, Nasslaborbiologen mit begrenzten Kenntnissen in der Bioinformatik in die Lage zu versetzen, ihre primären Transkriptomik-Daten auszuwerten. Rup (Abbildung 1) ist für Forscher zugänglich, die mit grundlegenden Kenntnissen von R vertraut sind. Die Ausführung des hier vorgestellten Workflows wird den Forschern ein detailliertes Verständnis ihrer Primärdaten vermitteln, einschließlich ihrer potenziellen Einschränkungen für die nachfolgende Analyse. Unseres Wissens fehlt es bisher an einer praktischen Pipeline zur Auswertung von primären Transkriptomdaten in Kombination mit Richtlinien zur Unterscheidung von Daten mit hoher und niedriger Qualität.

Abbildung 1: Arbeitsablauf der RNA-seq in silico Qualitätskontrollpipeline. Die Eingabedateien für die Qualitätskontrolle werden aus RNA-seq-Daten abgeleitet, die durch die Sequenzierung von Pflanzenmaterial generiert werden, sowie aus öffentlich zugänglichen Datensätzen. Die bereitgestellte R-Pipeline bewertet die Sequenzqualität anhand von drei Hauptansätzen: Sequenzierungsqualität, Mapping-Qualität und Replikationsqualität. Es werden verschiedene Qualitätsmetriken berechnet, und statistische Ergebnisse werden mit gängigen R-Paketen wie ggplot2 und pheatmap (z. B. Balkendiagramme und Heatmaps) visualisiert. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.
Zuvor gemeldete Auswertungspipelines erforderten vorverarbeitete Daten (z. B. Leseausrichtung als .bam-Dateien), decken nicht alle Metriken ab oder werden nicht mehr gepflegt 4,5,6. Der Vorteil des hier vorgestellten Workflows liegt in seiner Vollständigkeit durch die Konsolidierung der häufigsten Probleme der Qualitätskontrolle in einer einzigen Pipeline. Darüber hinaus stellen wir Beispiele für qualitativ hochwertige Daten vor, die für alle nachgelagerten Analyseanwendungen geeignet sind, aber auch Beispiele für Daten von geringer Qualität und diskutieren deren spezifische Einschränkungen für die weitere Analyse. Mehrere bereits veröffentlichte Berichte beschreiben den Zweck von RNA-seq-Analysewerkzeugen und bieten vergleichende Bewertungen ihrer Leistung 7,8. Die Standards für die RNA-seq-Qualitätskontrolle und die methodische Berichterstattung sind jedoch nicht standardisiert und erschweren die Reproduzierbarkeit und die biologisch sinnvolle Interpretation von Transkriptomik-Experimenten.
Rup integriert hochwertige Standardwerkzeuge, Qualitätskontrollanalysen und die Visualisierung der Ergebnisse. Rup erfordert Rohsequenzierungs-Reads und ein annotiertes Genom als Eingabe und läuft unter Mac OS, Linux und "Windows Subsystem for Linux" (WSL) auf Windows-Systemen. Es integriert Werkzeuge für die Qualität der Sequenzierung und quantifiziert das Mapping von Reads auf ein Genom, beinhaltet die Messung des rRNA-Gehalts in den Proben und bietet Probenkorrelation zur Bewertung der Replikatkorrelation. Die Testdaten wurden mittels Lasermikrodissektion des zentralen Meristemgewebes von zwei verschiedenen Stadien der Pflanzenart Eschscholzia californica, einem Protokoll mit extrem niedrigem Input für die Bibliotheksvorbereitung, gewonnen und auf Novaseq 6000 sequenziert.
Rup kann als einzelnes Skript mit minimalen Eingabedateien ausgeführt werden. Die Pipeline, wie hier gezeigt, benötigt Illumina-Paired-End-RNA-Sequenzierungsdaten. Eine angepasste Pipeline für Single-End-Sequenzierung mit den gleichen Analyseschritten ist ebenfalls im GitHub-Repository hinterlegt. Es werden nur die Genomsequenz als Fasta-Datei, das Genmodell und die rRNA-Annotationen als gtf-Dateien und die Rohsequenzierungs-Reads als fastq.gz-Dateien benötigt. Stellen Sie sicher, dass die Genom- und Annotationsdateien als genome.fa, annotation.gtf und rRNA.gtf in einem Ordner bereitgestellt werden, der in der Variablen reference_folder angegeben ist. Wenn alle Sequenzierungsdateien als .fq.gz Dateien in der Variablen read_file_folder gespeichert sind, kann Rup wie folgt ausgeführt werden:
1. Vorbereitungen:
BiocManager::install(c("getopt", "ggplot2", "reshape2", "pheatmap", "fastqcr", "Rfastp", "Rsubread", "Rsamtools"))
# 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)
n_threads <- 8 # the number of available CPU cores
bamSortMemory <- "1024" # maximum memory for bam file sorting
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))
2. Qualitätsbewertung der Sequenzierung
HINWEIS: In diesem Schritt wird ein Balkendiagramm erstellt, das die Anzahl der Lesevorgänge vor und nach dem Trimmen jeder Probe anzeigt.
# 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)
# 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)
}
# 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)
# 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)

Abbildung 2: Sequenzierungs- und Trimmergebnisse. Lesen Sie die Zählungen vor (rot) und nach dem Trimmen (grün). Sample-s2_r1 hat bereits vor dem Trimmen eine niedrige Lesezahl, während beim Trimmen ein großer Teil der Sample-s2_r2-Reads entfernt wurde. Alle anderen Beispiele zeigen einen akzeptablen Verlust von Lesevorgängen durch das Kürzen. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.
3. Qualität der Kartierung
HINWEIS: In diesem Schritt werden Balkendiagramme berechnet, um einfach kartierte Reads (z. B. Transkripte von proteinkodierenden Genen) von mehrfach kartierten Reads (z. B. rRNA-Reads) und nicht kartierten Reads (z. B. Kontaminationen) zu unterscheiden.
# 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)
# 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"))
}
# 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)
# 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)
# 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"
# 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"
# 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)
# 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)

Abbildung 3: Lesen Sie die Übersicht über das Mapping und die Statistiken. Die meisten Proben zeigen eine hohe Anzahl zugeordneter Reads (linke Seite, braun) und eine geringe Anzahl von rRNA-Reads (rechte Seite, blau). Probe s2_r3 weist eine ungewöhnlich hohe Anzahl von Multi-Map-Reads auf (linke Seite, grün), was einer hohen rRNA-Lesezahl (rechte Seite) entspricht. Probe s2_r4 zeigt eine hohe Anzahl von nicht kartierten Reads in Kombination mit einer erwarteten rRNA-Read-Zahl, was auf eine Kontamination mit Reads eines anderen Organismus hindeutet. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.
4. Replizieren Sie die Qualität.
# 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")

Abbildung 4: Klassifizierung der Gen-Read-Zahl. Alle Gene in einer Probe werden basierend auf der Anzahl der zugewiesenen Reads in eine von fünf Kategorien eingeteilt. In Rot ist die Anzahl der Gene dargestellt, denen keine Reads zugewiesen wurden. Proben mit einer geringen Anzahl von insgesamt zugewiesenen Reads (s2_r1 bis s1_r4) weisen eine höhere Anzahl von Genen mit 10 bis 100 zugewiesenen Reads und eine niedrigere Anzahl von Genen mit mehr als 1000 Reads auf. Gene mit geringer Expression können in diesen Proben übersehen werden. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.
Die Pipeline ist vollständig als R-Skript implementiert und wurde auf den Betriebssystemen Linux und Mac OS getestet. Windows-Benutzer können das Windows-Subsystem für Linux (WSL) verwenden. Der Code und die Testdaten sind als GitHub-Repository verfügbar: https://github.com/oliverrupp/rup. Die Sequenzierungsdaten sind im Rahmen des ENA EBI-Projekts PRJEB96400 verfügbar.
Zehn Proben wurden künstlich aus zwei realen Proben hergestellt, um verschiedene Probleme zu veranschaulichen, die bei der Qualitätskontrolle von Bulk-RNA-seq mit Rup zur Analyse auftreten können. Die Probe s2_r1 wurde so konzipiert, dass sie eine niedrige Gesamtlesezahl aufweist, wobei die Probe s2_r2 einen großen Anteil an Lesevorgängen von geringer Qualität enthält, die durch den Trimmprozess verworfen werden sollen. Die Probe s2_r3 enthält einen großen Teil der rRNA-Reads, und die Probe enthält s2_r4 Kontaminanten-Reads, die nicht dem Referenzgenom zugeordnet werden konnten. Die Namen der Stichproben s1_r5 und s2_r5 wurden vertauscht, um die geringe Replikatkorrelation zu veranschaulichen.
Protokollsektion 2 ermöglicht die Identifizierung von Proben mit niedrigen Lesezahlen vor oder nach dem Trimmen (Abbildung 2). Das Balkendiagramm zeigt die niedrigere Lesezahl in der Stichprobe s1_r1 sowohl vor als auch nach dem Trimmen. Hier war die anfängliche Lesezahl niedrig. Die geringe Anzahl von Samples, die nach dem Trimmen s1_r2 wurden, deutet auf eine große Anzahl von Adaptersequenzen, Sequenzierungsfehlern, Primersequenzen und sequenzierten Poly-A/T-Dehnungen hin, die während des Trimmvorgangs entfernt wurden. Der Abbau der Eingangs-RNA kann auch die Anzahl der qualitativ hochwertigen Lesevorgänge verringern.
Protokollabschnitt 3 identifiziert Probleme bei Lesezuweisungen. Abbildung 3 zeigt die erhöhte Anzahl von Multi-Mapping-Reads in der Probe s2_r3 (grün) sowie die hohe Anzahl von rRNA-Reads. Die Kontamination der Probe s2_r4 zeigt sich an dem großen Anteil der Reads, die nicht dem Referenzgenom zugeordnet sind (rosa). Diese Reads korrelieren nicht mit rRNA, sondern mit Sequenzen aus einem Nicht-Zielorganismus. Abbildung 4 zeigt allgemeine Probleme, die mit einer geringen Anzahl zugewiesener Reads in Transkriptomen verbunden sind. In Proben mit einer geringen Rate an Reads, die eindeutig auf das Referenzgenom abgebildet wurden (s2_r1 bis r4), haben nur etwa ein Drittel der Gene mehr als 100 Reads, während in den anderen Proben etwa die Hälfte der Gene in diese Klassen fällt. Reads von Genen mit sehr geringer Expression können in den Proben s2_r1 r4 nicht gefunden werden, so dass eine vergleichende Expressionsanalyse mit diesen Proben sehr unzuverlässig ist und vermieden werden sollte.
Protokollabschnitt 4 kann verwendet werden, um Replikatausreißer zu identifizieren. Es wird erwartet, dass Replikate derselben Probe/desselben Zustands/desselben Gewebes untereinander eine höhere Korrelation aufweisen als Replikate aus anderen Proben/Bedingungen/Geweben. Abbildung 5 zeigt eine Korrelations-Heatmap der beiden Stichproben (S1 und S2) mit jeweils fünf Replikaten. Die Dendrogramme oben und an der Seite des Diagramms zeigen zwei Cluster mit fünf Replikaten in jedem Cluster. Der linke Cluster enthält vier Replikate von Probe 1 und ein Replikat von Probe 2 (s2_r5), der rechte Cluster enthält vier Replikate von Probe 2 und ein Replikat von Probe 1 (s1_r5). Wenn in diesem Fall die Beispielnamen s1_r5 und s2_r5 wieder zurückgetauscht werden, enthält jeder Cluster alle Replikate eines Beispiels, was auf einen Fehler bei der Replikatbeschriftung hinweisen kann. Andere Gründe dafür, dass Replikate nicht zusammenwachsen, könnten eine mangelnde Differenzierung zwischen den Proben/Bedingungen/Geweben oder die Clusterbildung von Replikaten nur aus Gründen der Sequenzierungsqualität sein. Letzteres kann auftreten, wenn alle Replikate mit außergewöhnlich niedrigen oder außergewöhnlich hohen Lesezahlen nach dem Trimmen und Zuordnen einen Cluster bilden.

Abbildung 5: Beispiel-Korrelations-Heatmap. Die Beispielkorrelations-Heatmap basiert auf log2 transformierten TPM-Werten und zeigt zwei unterschiedliche Cluster mit jeweils fünf Stichproben. Es wird erwartet, dass biologische/technische Replikate eine höhere Korrelation untereinander aufweisen als solche aus anderen Geweben/Behandlungen. Der linke Cluster enthält vier Replikate von Probe 1 und ein Replikat von Probe 2 (s2_r5), der rechte Cluster enthält vier Replikate von Probe 2 und ein Replikat von Probe 1 (s1_r5). Die einzelnen Replikate sollten auf möglichen Sample-Swapping oder Batch-Effekte überprüft werden, um ihre Clusterbildung in der Heatmap zu erklären. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.
Tabelle 1: Vergleich von RNA-seq-Qualitätsbewertungspipelines. Für eine detaillierte Analyse siehe Ergänzende Datei 1. Bitte klicken Sie hier, um diese Tabelle herunterzuladen.
Ergänzende Datei 1: Auswahl von Parametern und Referenzen. Die Analysen zeigen den Einfluss der Parameter- und Referenzauswahl auf die gesamten QC-Ergebnisse. Bitte klicken Sie hier, um diese Datei herunterzuladen.
Die Qualität differentieller Genexpressionsanalysen hängt stark von zwei Faktoren ab: der Anzahl der in jeder Probe sequenzierten undreplizierten Reads 16 und der Anzahl der Replikate pro Probe17,18. Hier stellen wir die benutzerfreundliche Rup-Pipeline vor, um die Anzahl der Reads zu bestimmen, die für die Quantifizierung der Genexpression in jedem Replikat geeignet sind. Verschiedene Metriken ermöglichen es den Forschern zu verstehen, warum Replikate niedrige zugewiesene Lesezahlen aufweisen, und Probleme in der Stichprobenkorrelation zu identifizieren. Obwohl Rup für die Qualitätsbewertung von pflanzlichen RNA-seq-Proben entwickelt wurde, ist es für andere eukaryotische Organismen gleichermaßen geeignet und erfordert keine weiteren Anpassungen für nicht-pflanzliche Proben (Supplemental File 1).
Der erste Schritt von Rup bestimmt die Gesamtanzahl der Lesevorgänge vor und nach dem Kürzen der Lesevorgänge. Die Anzahl der sequenzierten Reads bestimmt die Anzahl der nachweisbaren Gene, und wenn die Read-Zahl zu niedrig ist, bleiben viele differentiell exprimierte Gene unentdeckt. Ein großer Teil der Reads, die während des Qualitätstrimm- und Filterprozesses verworfen werden, könnte auf eine RNA-Degradation hinweisen, was zu einer Unterschätzung der Expression bei stark degradierten Genen führen könnte. Ob eine Probe für nachgelagerte Anwendungen verwendet werden kann, hängt jedoch vom Zielorganismus und dem Ziel der Studie ab. Um beispielsweise die Expression der meisten Pflanzengene zu erfassen, sind unserer Erfahrung nach 30 bis 50 Millionen Reads erforderlich, während für Pilze nur 10 Millionen ausreichen könnten. Daher wird Rup keine Qualitätsschwellen für den Ausschluss problematischer Stichproben definieren, sondern Metriken zur Identifizierung verschiedener Probleme bereitstellen, die auftreten können.
Im zweiten Schritt (Mapping-Qualität) wird die Genauigkeit und Zuverlässigkeit des Read-Alignments mit dem Referenzgenom bewertet und deren Eignung für nachgelagerte Analysen spezifiziert. Nicht alle sequenzierten Reads können für die Berechnung der Genhäufigkeit verwendet werden. Reads, die nicht mit dem Genom oder Transkriptom übereinstimmen, oder Reads, die mehreren Stellen im Genom zugeordnet sind, werden in den meisten Fällenignoriert 19 und tragen nicht zur Gesamtanzahl der Reads bei. Die Ergebnisse des zweiten Schritts der Pipeline können verwendet werden, um zu verstehen, warum Lesevorgänge bei der Berechnung der Abundanz nicht verwendet werden. Eine hohe Anzahl nicht kartierter Reads könnte auf eine Kontamination während der RNA-Extraktion (z. B. mit Pflanzenpathogenen oder Pflanzenfressern) oder ein unvollständiges Referenzgenom hinweisen. Unvollständige Genmodelle des Referenzgenoms können zu einer hohen Anzahl von "no feature"-Reads führen. Ein großer Teil der mehrfach zugeordneten Reads könnte durch eine hohe Anzahl von rRNA-Genen in der Bibliothek verursacht werden, was auf eine unzureichende rRNA-Entfernung während des Vorbereitungsprozesses der Sequenzierungsbibliothek hinweist (für den Fall, dass die mehrfach zugeordneten Reads wirklich von rRNAs abgeleitet sind, können sie weiter analysiert werden, indem der Pipeline eine rRNA-Annotationsdatei zur Verfügung gestellt wird, die dann automatisch die Anzahl der möglichen rRNA-Reads berechnet).
Die dritte, ebenso wichtige Qualitätsmetrik ist die Korrelation zwischen Replikaten derselben Stichprobe. Im Allgemeinen sollte die Korrelation zwischen Replikaten derselben Stichprobe höher sein als die Korrelation zwischen Replikaten verschiedener Stichproben. Eine geringe Korrelation zwischen den Replikaten könnte auf eine große biologische Variabilität zwischen den Replikaten, eine hohe Ähnlichkeit zwischen den Proben/Bedingungen/Geweben, andere Chargeneffekte oder sogar auf einen Probentausch oder eine falsche Markierung hinweisen. Rup berechnet paarweise Korrelationen zwischen allen Stichproben und erstellt eine geclusterte Heatmap der Stichprobenkorrelationen. Zusätzlich kann eine Hauptkomponentenanalyse (PCA) an den Proben berechnet werden, um mögliche Chargeneffekte zu identifizieren, die in weiteren nachgelagerten Analysen berücksichtigt werden müssen. Obwohl Batch-Effekte in nachgelagerten Analysen korrigiert werden können, kann es besser sein, problematische Proben zu entfernen, wenn die gemessene Differenz zu hoch ist.
Das Einsatzmaterial hat einen direkten Einfluss auf die Sequenzierungsqualität: Gewebeprobenahme, RNA-Extraktion und Bibliotheksvorbereitung sind entscheidende Schritte zur Minimierung von Problemen mit der RNA-seq-Qualität. Es sollten gleichbleibende Bedingungen aufrechterhalten werden, wenn möglich, zum Beispiel durch den Einsatz von Wachstumskammern. Schädlingsbekämpfungsmaßnahmen sollten rechtzeitig durchgeführt werden, da nur gesunde Personen für die Probenahme ausgewählt werden sollten. Es ist weiterhin ratsam, Proben am selben Tag und/oder zur gleichen Zeit zu entnehmen, um die Variation des zirkadianen Transkriptoms zu reduzieren. Es stehen zahlreiche RNA-Extraktionskits zur Verfügung, und die Auswahl eines geeigneten Kits für die Zielspezies kann die mRNA-Qualität verbessern. Die Protokolle zur Bibliotheksvorbereitung sollten Schritte zur Anreicherung von polyA+-RNA enthalten, um die rRNA-Fraktion zu minimieren.
Rup kann als erster Schritt der Qualitätskontrolle in der differentiellen Genexpressionsanalyse verwendet werden. Diese Qualitätskontrolle ist aus mehreren kritischen Gründen erforderlich, da sie die Qualität der Eingabedaten sicherstellt (identifiziert minderwertige Replikate/Proben, kann Qualitätsschwellenwerte für die Lesetiefe und die Mapping-Raten festlegen) und technische Probleme wie Sequenzierungsfehler, Batch-Effekte oder falsche Beschriftungen identifiziert. Wenn die RNA-Eingabe von ausreichender Qualität ist, kann Rup helfen, indem es minderwertige Reads kürzt oder falsch markierte Proben identifiziert. Die schlechte Qualität der Eingangs-RNA kann Rup jedoch nicht kompensieren. Bei minderwertiger RNA oder fehlerhaften Sequenzierungsdaten kann es erforderlich sein, Proben erneut zu entnehmen, das RNA-Extraktionsprotokoll anzupassen und/oder die Sequenzierung zu wiederholen. Gleiches gilt für Proben mit einem großen Anteil an nicht kartierten Reads, die möglicherweise durch Kontamination verursacht wurden. Aus unserer Erfahrung können jedoch nicht alle RNA-Extraktionen aufgrund mangelnder Materialverfügbarkeit einfach wiederholt werden. In diesem Fall sind sie für einige Downstream-Anwendungen ungeeignet: Proben mit einer geringen Anzahl von Single-Mapping-Reads sollten nicht in differentiellen Genexpressionsanalysen analysiert werden, enthalten aber dennoch Informationen für die Analyse des Vorhandenseins von Transkripten. In diesem Fall kann das Fehlen von Transkripten in Geweben/Behandlungen/Zuständen für die Analyse und Quantifizierung der Transkripthäufigkeit nicht berücksichtigt werden.
Rup enthält einige Einschränkungen. Zum Beispiel wird nicht direkt auf RNA-Abbau getestet, da vor der Vorbereitung und Sequenzierung der Bibliothek direkte Messungen der RNA-Integrität erforderlich sind. Ein Skript zum Identifizieren des RNA-Abbaus mithilfe der RSeQC-Module geneBodyCoverage.py und tin.py ist jedoch im Repository dieser Pipeline enthalten. Darüber hinaus testet Rup nicht auf GC-Gehalt und Verzerrung der Transkriptlänge. Die Größe des Referenzgenoms ist derzeit auf 4 GB begrenzt, und nach unserer Erfahrung überschätzt das Read-Mapping-Modul die Anzahl der mehrfach kartierten Reads in Polyploiden, wie in der Supplemental File am Beispiel eines Vergleichs des Read-Mappings zu einer haploiden vs. den diploiden Genomversionen von E. californica gezeigt wird, ein haploides Genom ist daher der bevorzugte Input für diese Pipeline. Die Rup-Ausgabequalität hängt weitgehend von der Qualität des Referenzgenoms und der Annotation ab, so dass eine unvollständige oder stark fragmentierte Genomsequenz zu einer Überschätzung der nicht kartierten Reads führen kann. Darüber hinaus führt eine unvollständige Annotation von Genmodellen zu einer Überschätzung von nicht zugeordneten Reads. Die Vollständigkeit und Duplizierung der Genomsequenz und der Genannotation kann mit Tools wie BUSCO20 abgeleitet werden.
Rup ist ein eigenständiges Werkzeug für alle anfänglichen Qualitätskontrollschritte, die für RNA-seq-Experimente unerlässlich sind. Im Gegensatz zu RSeQC und RNA-SeQC ist keine Vorverarbeitung der Rohsequenzierungs-Reads erforderlich. Die Qualitätsanalyse der Sequenzierung kann nicht mit RNA-SeQC durchgeführt werden, und Probenkorrelations-Heatmaps werden nicht mit RSeQC und RNA-SeQC berechnet (Tabelle 1). Darüber hinaus erfolgt die Normalisierungsausgabe in FPKM in RSeQC, und die Ausgabevisualisierung ist nicht als Standard für alle Module in RSeQC und RNA-SeQC implementiert. Daher werden mehrere Probleme der Qualitätskontrolle in den beiden alternativen Tools nicht getestet, darunter eine niedrige Lesezahl, ein hoher Anteil an getrimmten Lesevorgängen und die Identifizierung von Replikatausreißern. Ein Vergleich von Rup, RSeQC und RNA-SeQC ist in der Supplemental File 1 verfügbar. Darüber hinaus kann Rup komplementär zu RNA-SeQC oder RSeQC verwendet werden, z. B. können die von dieser Pipeline erzeugten sortierten BAM-Dateien als Input für diese Tools verwendet werden.
Derzeit ist Rup für eine kleine Anzahl von Samples optimiert, aber die zeitaufwändigsten Schritte, wie z. B. Quality Trimming und Read Mapping, können beispielsweise auf einem Computercluster oder einer Cloud-Infrastruktur vorberechnet werden. Für Genome, die größer als 4 Gb sind, ist dies obligatorisch, da der Rsubread-Aligner auf Genome beschränkt ist, die kleiner als 4 Gb sind. Die rRNA-Annotation erfolgt mit Barrnap, das die hochkonservierten rRNA-Gene identifiziert. Nach unserer Erfahrung erfordert die Annotation von rRNA-Genen keine Vollständigkeit, denn selbst wenn einige abnormale rRNA-Gene fehlen, ist ein grober Überblick über das Vorhandensein von rRNA im Datensatz für die Beurteilung der Datensatzqualität ausreichend. Zukünftige Versionen von Rup könnten auf eine andere Mapping-Methode umsteigen, wie z.B. das Pseudo-Alignment-Tool salmon21, um die Laufzeit zu verkürzen. Darüber hinaus könnten weitere Normalisierungs- und Korrekturmethoden, wie z. B. GC-Bias oder Längen-Bias-Korrektur, zu Rup hinzugefügt werden.
Zusammenfassend lässt sich sagen, dass Rup wesentliche Informationen liefert, um die Zuverlässigkeit und Reproduzierbarkeit der RNA-seq-Datenanalyse zu gewährleisten. Diese Pipeline berichtet umfassend über die wichtigsten RNA-seq-Qualitätsmetriken und erstellt Ausgabedateien für die direkte Verwendung in nachgelagerten Analysen. Es wurde als eigenständiges Tool für Forscher mit minimalen Bioinformatik-Kenntnissen entwickelt, um die Qualität ihrer primären Sequenzierungsdaten mit intuitiver Visualisierung zu bewerten.
Die Autoren haben keine Interessenkonflikte anzugeben.
Wir danken für die technische Unterstützung durch die Bioinformatics Core Facility an der Professur für Systembiologie der JLU Gießen sowie für die Bereitstellung von Rechenressourcen und die allgemeine Unterstützung durch das BiGi Service Center (BMFB-Stipendium 031A533) innerhalb der DE. NBI-Netzwerk. Die hier vorgestellte Arbeit wurde durch das Stipendium der Deutschen Forschungsgemeinschaft (DFG) BE2547/24-1 an A.B. gefördert, und wir danken auch der Justus-Liebig-Universität Gießen für die Unterstützung.
| Name | Unternehmen | Katalognummer | Kommentare |
|---|---|---|---|
| fastQCR | R | 0.1.3 | |
| FUJITSU ESPRIMO D958 Desktop | FUJITSU | die Pipeline wurde auf Ubuntu Linux 24.02 (32 GB RAM) entwickelt und getestet | |
| Getopt | R | 1.20.4 | |
| ggplot2 | R | 3.5.2 | |
| MacBook Pro | Apfel | die Pipeline wurde auf maxOS 15.4.1 (16 GB RAM) getestet | |
| Pheatmap | R | 1.0.13 | |
| R | 4.4.3 | ||
| Reshape2 | R | 1.4.4 | |
| rfastp | Bioleiter | 1.16.0 | |
| rsamtools | Bioleiter | 2.22.0 | |
| rsubread | Bioleiter | 2.20.0 |
Genehmigung beantragen, um den Text oder die Abbildungen dieses JoVE-Artikels zu verwenden
Genehmigung beantragen