Rup (RNA-seq Usability Assessment Pipeline) - Controllo di qualità per esperimenti di RNA-seq di massa negli eucarioti

Sommario

Questo protocollo consente il controllo di qualità iniziale per esperimenti di RNA-seq per biologi di laboratorio con esperienza limitata in bioinformatica.

Abstract

Gli approcci moderni nella scienza molecolare delle piante spesso richiedono esperimenti di RNA-seq di massa, ad esempio, per tracciare i cambiamenti globali nei trascrittomi dopo i trattamenti o per identificare i componenti chiave dei percorsi regolatori. Di conseguenza, diverse aree delle scienze vegetali si affidano a dati di RNA-seq di massa riproducibili e di alta qualità per il progresso scientifico. Tuttavia, in base alla nostra esperienza, la conoscenza e l'applicazione delle misure di controllo della qualità nei set di dati RNA-seq sono spesso carenti. Qui, introduciamo Rup (RNA-seq usability assessment pipeline) per il controllo di qualità di dati di RNA-seq di massa, per successive analisi di espressione genica, che è autonomo e facilmente applicabile per i biologi di laboratorio con conoscenze di base di R. Rup aiuta a discriminare tra dati di sequenziamento di alta qualità, adatti per esperimenti di espressione genica a valle, e quelli non adatti per ulteriori analisi generali. Rup include test per diversi problemi comunemente riscontrati come numeri di lettura o mappatura insufficienti, identificazione di contaminazioni, quantificazione di frazioni di rRNA nei dati totali RNA-seq, test di somiglianza replicata, utilizzando dati reali per la dimostrazione e offrendo una visualizzazione intuitiva. Rup fornisce una suite di strumenti per identificare le carenze sperimentali prima dell'analisi standardizzata del trascrittoma, migliorando così la qualità dei dati per i singoli ricercatori e per il campo. Ciò aumenta la fiducia nell'analisi dei dati di RNA-seq di massa e fornisce una base per le future linee guida che definiscono i criteri minimi di controllo della qualità, migliorando così l'affidabilità e la trasparenza dei dati RNA-seq pubblicati. I dati sui test e sui test sono disponibili all'https://github.com/oliverrupp/rup.

Introduzione

Gli esperimenti di trascrittomica (RNA-seq) interrogano in modo completo le firme trascrizionali che modellano i fenotipi. Questo approccio è diventato insostituibile nella genetica molecolare delle piante per identificare geni singoli o co-espressi e processi biologici coinvolti, ad esempio, nello sviluppo, nell'interazione tra patogeni delle piante o nella resistenza agli stress abiotici^{1, 2}, ³. I recenti progressi nella tecnologia RNA-seq hanno aumentato la specificità e hanno permesso il rilevamento di diverse isoforme e varianti a risoluzione a base singola, consentendo l'identificazione di variazioni di sequenza da indel più grandi a polimorfismi a singolo nucleotide (SNP). I dati ottenuti da RNA-seq sono caratterizzati da un ampio intervallo dinamico, che facilita la rilevazione di trascritti sia altamente abbondanti che a bassa espressione, e richiedono condizioni sperimentali appropriate per la coerenza. Inoltre, i set di dati RNA-seq sono di grandi dimensioni e quindi ad alta intensità di calcolo da analizzare e archiviare, il che comporta una gestione efficiente dei dati e ampie risorse di calcolo. L'RNA-seq richiede un input di alta qualità in ogni fase del flusso di lavoro, poiché i punti deboli in qualsiasi fase possono propagarsi e compromettere i risultati. Una scarsa integrità dell'RNA o problemi tecnici durante la preparazione della libreria porteranno a distorsioni e a una riduzione dell'accuratezza. I parametri per ottenere letture grezze di alta qualità possono variare tra le strutture di sequenziamento di nuova generazione (NGS).

In base alla nostra esperienza, si consiglia di utilizzare solo RNA con un numero di integrità dell'RNA (RIN) superiore a 7, che è indicativo di una struttura dell'mRNA in gran parte intatta, come input per la preparazione della libreria di RNA-seq. Per il successo del sequenziamento, sono necessari circa 2 μg di RNA totale a una concentrazione di 50-200 ng/μL per i protocolli standard di preparazione delle librerie. La purezza dell'RNA dovrebbe essere confermata da un rapporto OD_260/280 compreso tra 1,8 e 2,1 e da un rapporto OD_260/230 maggiore di 1,5 utilizzando uno spettrofotometro. Una profondità di sequenziamento insufficiente, bassi tassi di mappatura o disallineamento rispetto al genoma di riferimento possono distorcere ulteriormente l'espressione genica e la quantificazione dello splicing. Inoltre, un disegno sperimentale inadeguato, come la mancanza di discriminazione tra campioni di tessuti o trattamenti diversi e l'elevata variabilità tra le repliche, può introdurre rumore e diminuire la riproducibilità. Sebbene molti laboratori analizzino regolarmente i trascrittomi, i rigorosi controlli di qualità delle prime fasi essenziali dell'analisi spesso non vengono riportati o possono essere del tutto assenti dalle pubblicazioni. Ciò può portare a una sovrainterpretazione dei risultati derivati dall'analisi del trascrittoma e, di conseguenza, a risultati irriproducibili.

Qui, forniamo un flusso di lavoro per il controllo di qualità delle fasi iniziali necessarie per analisi del trascrittoma di alta qualità dei profili mRNA per misurare le alterazioni trascrizionali. Il nostro obiettivo è quello di consentire ai biologi wet-lab con conoscenze limitate in bioinformatica di valutare i loro dati trascrittomici primari. Rup (Figura 1) è accessibile ai ricercatori che hanno familiarità con le conoscenze di base di R. L'esecuzione del flusso di lavoro qui presentato fornirà ai ricercatori una comprensione dettagliata dei loro dati primari, comprese le loro potenziali limitazioni per l'analisi successiva. Per quanto ne sappiamo, finora manca una pipeline pratica di valutazione dei dati del trascrittoma primario in combinazione con linee guida per distinguere i dati di alta qualità da quelli di bassa qualità.

Figura 1: Flusso di lavoro della pipeline di controllo qualità RNA-seq in silico . I file di input per il controllo di qualità derivano dai dati RNA-seq generati dal sequenziamento del materiale vegetale, nonché da set di dati disponibili pubblicamente. La pipeline R fornita valuta la qualità della sequenza attraverso tre approcci principali: qualità del sequenziamento, qualità della mappatura e qualità della replica. Vengono calcolate varie metriche di qualità e i risultati statistici vengono visualizzati utilizzando pacchetti R comuni come ggplot2 e pheatmap (ad esempio, grafici a barre e mappe termiche). Clicca qui per visualizzare una versione più grande di questa figura.

Le pipeline di valutazione segnalate in precedenza richiedevano dati pre-elaborati (ad esempio, leggere l'allineamento come file .bam), non coprono tutte le metriche o non sono più mantenute ^4,5,6. Il vantaggio del flusso di lavoro qui presentato risiede nella sua completezza grazie al consolidamento dei problemi di controllo qualità più comuni in un'unica pipeline. Inoltre, forniamo esempi di dati di alta qualità adatti a tutte le applicazioni di analisi a valle, ma anche esempi di dati di bassa qualità e discutiamo i loro limiti specifici per ulteriori analisi. Diversi rapporti pubblicati in precedenza descrivono lo scopo degli strumenti di analisi dell'RNA-seq e forniscono valutazioni comparative delle loro prestazioni ^7,8. Tuttavia, gli standard di controllo della qualità e di segnalazione della metodologia RNA-seq non sono standardizzati e aggravano la riproducibilità e le interpretazioni biologicamente significative degli esperimenti di trascrittomica.

Rup integra strumenti standard di alta qualità, analisi del controllo qualità e visualizzazione dei risultati. Rup richiede letture di sequenziamento non elaborate e un genoma annotato come input ed è eseguito su Mac OS, Linux e "Windows Subsystem for Linux" (WSL) su sistemi Windows. Integra strumenti per la qualità del sequenziamento e quantifica la mappatura di lettura in un genoma, include la misurazione del contenuto di rRNA nei campioni e offre la correlazione dei campioni per valutare la correlazione della replica. I dati del test sono stati ottenuti utilizzando la microdissezione laser del tessuto meristematico centrale di due diversi stadi della specie vegetale Eschscholzia californica, un protocollo a bassissimo input per la preparazione di librerie, e sequenziati su Novaseq 6000.

Protocollo

Rup può essere eseguito come un singolo script con un numero minimo di file di input richiesti. La pipeline, come mostrato qui, richiede i dati di sequenziamento dell'RNA paired-end Illumina. Nel repository GitHub viene depositata anche una pipeline rettificata per il sequenziamento single-end con gli stessi passaggi di analisi. Sono necessari solo la sequenza del genoma come file fasta, il modello genico e le annotazioni dell'rRNA come file gtf e il sequenziamento grezzo come file fastq.gz. Assicurati che il genoma e i file di annotazione siano forniti come genome.fa, annotation.gtf e rRNA.gtf in una cartella specificata nella variabile reference_folder. Quando tutti i file di sequenziazione come file .fq.gz vengono salvati nella variabile read_file_folder, Rup può essere eseguito come segue:

1. Preparativi:

1. Installare tutti i pacchetti R necessari usando Bioconductor.

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

2. Imposta tutti i file necessari.
   1. Crea una cartella di origine che conterrà tutti i file necessari. Aggiungi la sequenza del genoma di riferimento nel file fasta chiamato "reference/genome.fa" alla cartella di origine. Fornire l'annotazione del modello genico come file gtf: "reference/annotation.gtf", Facoltativamente, fornire l'annotazione del gene rRNA come file gtf: "reference/rRNA.gtf".
   2. Aggiungere tutte le letture di sequenziamento come file fastq alla cartella "reads". Assicurarsi che tutti i nomi dei file fastq seguano lo stesso schema; _1.fastq.gz e _2.fastq.gz per le letture in avanti e all'indietro.

      # 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. Impostare il parametro in base alle risorse di calcolo.

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

4. Impostare il parametro in base al metodo di sequenziamento.
    

   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))

   
   
   NOTA: la lunghezza minima di lettura non deve essere inferiore a 25 bp, poiché letture più piccole potrebbero causare l'interruzione della mappatura. Valori più grandi si tradurranno in letture mappate in modo più univoco, ma anche in letture più scartate durante il taglio, a seconda della qualità della sequenziazione. Maggiori dettagli su questi parametri sono descritti negli integratori.

2. Valutazione della qualità del sequenziamento

NOTA: Questo passaggio creerà un grafico a barre che mostra il numero di letture prima e dopo il taglio di ciascun campione.

1. Eseguire fastqc⁹, che è uno strumento per il controllo generale della qualità dei dati di sequenza ad alto rendimento, che fornisce informazioni sul numero totale di letture sequenziate, la lunghezza di lettura, il contenuto medio di GC, la contaminazione dell'adattatore di sequenziamento e le sequenze sovrarappresentate (ad esempio, letture di rRNA). Analizza la qualità del sequenziamento in base ai rapporti sulla qualità di base e per sequenza. Eseguire fastqc sui file di sequenziazione non elaborati utilizzando il pacchetto R fastqcr¹⁰ e inserire tutti i file di input nella stessa directory (read_file_folder). Riepilogare i file di output da fastqc_folder in un oggetto R usando la funzione qc_aggregate.
   NOTA: I risultati verranno salvati nella directory fastqc_folder. È possibile accedere ai file .html creati da fastqc nella cartella di output per ulteriori statistiche.
    

   # 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. Con tutti i prefissi dei campioni memorizzati nel sample_prefixes vettoriale, iterare sui prefissi ed eseguire lo strumento di ritaglio della qualità fastp su tutti i campioni. Trim per rimuovere sequenze adattatrici, inneschi, basi di bassa qualità, sequenze poli-A/T e viene eseguito da Fastp¹¹, disponibile nel pacchetto R Rfastp¹². I file fastq tagliati verranno memorizzati nella cartella 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. Eseguire nuovamente fastqc sulle letture tagliate.
    

   # 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. Valutare il trimming raccogliendo il numero di letture prima e dopo il trimming per tutti i campioni. Crea un grafico a barre con i numeri letti utilizzando 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)

Figura 2: Risultati del sequenziamento e del taglio. Leggere i conteggi prima (rosso) e dopo il taglio (verde). Il s2_r1 campione ha già un basso numero di letture prima del taglio, mentre il taglio ha rimosso una grande frazione delle letture del s2_r2 campione. Tutti gli altri campioni mostrano una perdita accettabile di letture dovute al taglio. Clicca qui per visualizzare una versione più grande di questa figura.

3. Qualità della mappatura

NOTA: Questo passaggio calcola i grafici a barre per discriminare le letture a mappatura singola (ad esempio, trascrizioni da geni codificanti proteine) da letture multi-mappate (ad esempio, letture di rRNA) e letture non mappate (ad esempio, contaminazioni).

1. Utilizzare il pacchetto Rsubread¹⁴ per mappare le letture al genoma di riferimento. Innanzitutto, costruisci un indice del file fasta del genoma di riferimento (genome_fasta_file).
   NOTA: Questo passaggio viene eseguito una volta per ogni genoma di riferimento.
    

   # 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. Itera su tutti i campioni e allinea le letture al genoma di riferimento utilizzando la funzione align() del pacchetto Rsubread. Questa funzione crea file .bam nella cartella di output.
    

   # 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. Conta le letture per ogni gene utilizzando la funzione featureCounts() dal pacchetto Rsubread. Assicurarsi che il file di annotazione (annotation_file) sia in formato GTF. Il conteggio legge solo con una singola corrispondenza con il genoma.
    

   # 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. Contare le reads mappate sui geni dell'rRNA per valutare il contenuto di rRNA nei campioni. Consenti il conteggio delle letture multimappate qui. Specificare i loci del gene rRNA in un file GTF (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. Raccogli le statistiche di assegnazione delle letture create da featureCounts. Le statistiche includono il numero di allineamenti per ciascun campione nelle diverse categorie (ad esempio, Assegnato, Non mappato, Multimappato).
    

   # 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. Raccogli le statistiche dell'assegnazione dell'rRNA.
    

   # 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. Traccia le statistiche di mappatura delle letture come un grafico a barre.
    

   # 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. Classificare i geni in gruppi in base al numero di letture assegnate e tracciare i risultati come un grafico a barre.
    

   # 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)

Figura 3: Leggi la panoramica e le statistiche della mappatura. La maggior parte dei campioni mostra un numero elevato di letture assegnate (lato sinistro, marrone) e un numero basso di letture di rRNA (lato destro, blu). Il campione s2_r3 ha una quantità insolitamente elevata di letture multi-mappate (lato sinistro, verde) corrispondente a un numero elevato di letture di rRNA (lato destro). Il campione s2_r4 mostra un numero elevato di letture non mappate in combinazione con un numero di letture di rRNA attese, indicativo di contaminazione con letture di un altro organismo. Clicca qui per visualizzare una versione più grande di questa figura.

4. Replica la qualità.

1. Assicurarsi che le repliche dello stesso campione o tessuto mostrino una correlazione più elevata tra loro rispetto alle repliche di campioni o tessuti diversi. Tracciare una mappa termica raggruppata delle correlazioni replicate per visualizzare i dati utilizzando i conteggi di lettura non elaborati o i conteggi normalizzati TPM (Transcripts Per Million)¹⁵ .
    

   # 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")

Figura 4: Classificazione del conteggio delle letture geniche. Tutti i geni in un campione sono classificati in una delle cinque categorie in base al numero di letture assegnate. In rosso sono mostrati il numero di geni a cui non sono assegnate alcuna lettura. I campioni con un numero basso di letture totali assegnate (da s2_r1 a s1_r4) hanno un numero maggiore di geni con 10-100 letture assegnate e un numero inferiore di geni con più di 1000 letture. I geni con bassa espressione possono non essere presenti in questi campioni. Clicca qui per visualizzare una versione più grande di questa figura.

Risultati

La pipeline è completamente implementata come script R ed è stata testata nei sistemi operativi Linux e Mac OS. Gli utenti Windows possono utilizzare il sottosistema Windows per Linux (WSL). Il codice e i dati di test sono disponibili come repository GitHub: https://github.com/oliverrupp/rup. I dati di sequenziamento sono disponibili nell'ambito del progetto EBI ENA PRJEB96400.

Dieci campioni sono stati creati artificialmente da due campioni reali per esemplificare i diversi problemi che possono essere incontrati durante il controllo di qualità dell'RNA-seq di massa utilizzando Rup per l'analisi. Il s2_r1 campione è stato progettato per mostrare un basso numero totale di letture, il campione s2_r2 contiene una grande frazione di letture di bassa qualità da scartare dal processo di rifinitura. Il s2_r3 del campione include un'ampia frazione di letture di rRNA e il s2_r4 del campione include letture di contaminanti che non sono riusciti a mappare il genoma di riferimento. I nomi dei campioni s1_r5 e s2_r5 sono stati scambiati per illustrare una bassa correlazione di replica.

La sezione 2 del protocollo consente l'identificazione di campioni con numeri di lettura bassi prima o dopo il taglio (Figura 2). Il grafico a barre mostra il numero di lettura inferiore nel campione s1_r1 sia prima che dopo il taglio. In questo caso, il numero di letture iniziali era basso. Il basso numero di s1_r2 di lettura del campione dopo il trimming suggerisce una grande quantità di sequenze adattatrici, errori di sequenziamento, sequenze di primer, tratti poli-A/T sequenziati che sono stati rimossi durante il processo di trimming. La degradazione dell'RNA in ingresso può anche ridurre il numero di letture di alta qualità.

La sezione 3 del protocollo identifica i problemi nelle assegnazioni di lettura. La Figura 3 illustra l'elevato numero di letture multi-mappate nel campione di s2_r3 (verde), nonché l'elevato numero di letture di rRNA. La contaminazione del s2_r4 del campione è evidente dall'ampia frazione di reads non mappate sul genoma di riferimento (rosa). Queste letture non sono correlate all'rRNA ma alle sequenze di un organismo non bersaglio. La Figura 4 mostra i problemi generali associati a un basso numero di letture assegnate nei trascrittomi. Nei campioni con un basso tasso di reads mappati in modo univoco al genoma di riferimento (da s2_r1 a r4) solo circa un terzo dei geni ha più di 100 reads assegnate, mentre negli altri campioni, circa la metà dei geni rientra in queste classi. Le letture di geni con espressione molto bassa potrebbero non essere trovate nei campioni s2_r1 a r4 e, di conseguenza, l'analisi comparativa dell'espressione con questi campioni sarà altamente inaffidabile e dovrebbe essere evitata.

La sezione 4 del protocollo può essere utilizzata per identificare i valori anomali replicati. Ci si aspetta che le repliche dello stesso campione/condizione/tessuto mostrino una correlazione più elevata tra loro rispetto alle repliche di altri campioni/condizioni/tessuti. La Figura 5 mostra una heatmap di correlazione dei due campioni (S1 e S2) con cinque repliche ciascuno. I dendrogrammi in alto e a lato del grafico mostrano due cluster con cinque repliche in ciascun cluster. Il cluster di sinistra contiene quattro repliche del campione 1 e una replica del campione 2 (s2_r5), il cluster di destra contiene quattro repliche del campione 2 e una replica del campione 1 (s1_r5). In questo caso, quando i nomi dei campioni s1_r5 e s2_r5 vengono nuovamente scambiati, ogni cluster contiene tutte le repliche di un campione, il che potrebbe indicare un errore di etichettatura della replica. Altri motivi per cui le repliche non si raggruppano insieme potrebbero essere la mancanza di differenziazione tra i campioni/condizioni/tessuti, o il raggruppamento delle repliche solo per motivi di qualità del sequenziamento. Quest'ultimo può verificarsi quando tutte le repliche con un numero di letture eccezionalmente basso o quelle con un numero di letture eccezionalmente alto dopo il taglio e la mappatura formano un cluster.

Figura 5: Mappa termica di correlazione del campione. La mappa termica di correlazione del campione si basa sui valori TPM trasformati del log₂ e mostra due cluster distinti di cinque campioni ciascuno. Si prevede che le repliche biologiche/tecniche mostrino una correlazione più elevata tra loro rispetto a quelle di altri tessuti/trattamenti. Il cluster di sinistra contiene quattro repliche del campione 1 e una replica del campione 2 (s2_r5), il cluster di destra contiene quattro repliche del campione 2 e una replica del campione 1 (s1_r5). Le singole repliche devono essere controllate per verificare la presenza di possibili scambi di campioni o effetti batch per spiegare il loro raggruppamento nella heatmap. Clicca qui per visualizzare una versione più grande di questa figura.

Tabella 1: Confronto delle pipeline di valutazione della qualità dell'RNA-seq. Per un'analisi dettagliata, vedere File supplementare 1. Clicca qui per scaricare questa tabella.

File supplementare 1: Selezione di parametri e riferimenti. Le analisi mostrano l'influenza della selezione dei parametri e dei riferimenti sui risultati complessivi del controllo qualità. Clicca qui per scaricare questo file.

Discussione

La qualità delle analisi di espressione genica differenziale dipende fortemente da due fattori: il numero di reads sequenziate in ciascun campione e replicato¹⁶ e il numero di repliche per campione^17,18. Qui, presentiamo la pipeline Rup di facile utilizzo per determinare il numero di reads adatte alla quantificazione dell'espressione genica in ciascuna replica. Diverse metriche consentono ai ricercatori di capire perché le repliche mostrano un basso numero di letture assegnate e di identificare i problemi nella correlazione del campione. Sebbene Rup sia stato sviluppato per la valutazione della qualità di campioni di RNA-seq vegetali, è ugualmente adatto per altri organismi eucarioti, non richiedendo ulteriori aggiustamenti per campioni non vegetali (File supplementare 1).

Il primo passaggio di Rup determina il conteggio totale delle letture prima e dopo il taglio delle letture. Il numero di letture sequenziate determina il numero di geni rilevabili e, se il numero di letture è troppo basso, molti geni espressi in modo differenziale rimarranno non rilevati. Una grande frazione di reads scartate durante il processo di trimming e filtraggio della qualità potrebbe indicare una degradazione dell'RNA, che potrebbe portare a una sottostima dell'espressione di geni altamente degradati. Tuttavia, la possibilità di utilizzare un campione per applicazioni a valle dipende dall'organismo bersaglio e dall'obiettivo dello studio. Ad esempio, per catturare l'espressione della maggior parte dei geni delle piante, nella nostra esperienza, sono necessari da 30 a 50 milioni di read, mentre per i funghi potrebbero essere sufficienti solo 10 milioni. Pertanto, Rup non definirà soglie di qualità per l'esclusione di campioni problematici, ma piuttosto fornirà metriche per identificare i diversi problemi che possono essere incontrati.

La seconda fase (mappatura della qualità) valuta l'accuratezza e l'affidabilità dell'allineamento delle letture al genoma di riferimento, specificandone l'idoneità per le analisi a valle. Non tutte le letture sequenziate possono essere utilizzate per il calcolo dell'abbondanza genica; Le letture che non si allineano al genoma o al trascrittoma o le letture che mappano più posizioni sul genoma vengono ignorate nella maggior parte dei casi¹⁹ e non contribuiscono al conteggio complessivo delle letture. I risultati del secondo passaggio della pipeline possono essere utilizzati per comprendere il motivo per cui le letture non vengono utilizzate nel calcolo dell'abbondanza. Un numero elevato di letture non mappate potrebbe indicare una contaminazione durante l'estrazione dell'RNA (ad esempio, con patogeni delle piante o erbivori) o un genoma di riferimento incompleto. Modelli genetici incompleti del genoma di riferimento possono portare a un numero elevato di letture "senza funzionalità". Un'ampia frazione di reads multi-mappate potrebbe essere causata da un numero elevato di geni rRNA nella libreria che indicano una rimozione insufficiente di rRNA durante il processo di preparazione della libreria di sequenziamento (nel caso in cui le reads multi-mappate siano veramente derivate da rRNA, possono essere analizzate ulteriormente fornendo un file di annotazione rRNA alla pipeline, che quindi calcola automaticamente il numero di possibili reads di rRNA).

La terza metrica di qualità, altrettanto importante, è la correlazione tra le repliche dello stesso campione. In generale, la correlazione tra repliche di uno stesso campione dovrebbe essere superiore alla correlazione tra repliche di campioni diversi. Una bassa correlazione tra le repliche potrebbe indicare un'ampia variabilità biologica tra le repliche, un'elevata somiglianza tra campioni/condizioni/tessuti, altri effetti batch o persino lo scambio di campioni o l'etichettatura errata. Rup calcola le correlazioni a coppie tra tutti i campioni e produce una mappa di calore raggruppata delle correlazioni dei campioni. Inoltre, è possibile calcolare un'analisi dei componenti principali (PCA) sui campioni per identificare i possibili effetti batch che devono essere riconosciuti in ulteriori analisi a valle. Sebbene gli effetti dei lotti possano essere corretti nelle analisi a valle, potrebbe essere meglio rimuovere i campioni problematici se la differenza misurata è troppo elevata.

Il materiale in ingresso ha un impatto diretto sulla qualità del sequenziamento: il campionamento dei tessuti, l'estrazione dell'RNA e la preparazione della libreria sono passaggi fondamentali per ridurre al minimo i problemi di qualità dell'RNA-seq. Le condizioni costanti dovrebbero essere mantenute, se possibile, ad esempio utilizzando camere di crescita. Le misure di controllo dei parassiti devono essere eseguite in modo tempestivo, poiché solo gli individui sani devono essere selezionati per il campionamento. Si consiglia inoltre di raccogliere i campioni lo stesso giorno e/o alla stessa ora per ridurre la variazione del trascrittoma circadiano. Sono disponibili numerosi kit per l'estrazione dell'RNA e la selezione di un kit appropriato per le specie target può migliorare la qualità dell'mRNA. I protocolli di preparazione della libreria dovrebbero includere passaggi per l'arricchimento dell'RNA polyA+ per ridurre al minimo la frazione di rRNA.

Rup può essere utilizzato come fase iniziale del controllo di qualità nell'analisi dell'espressione genica differenziale. Questo controllo di qualità è necessario per diversi motivi critici, in quanto garantisce la qualità dei dati di input (identifica repliche/campioni di bassa qualità, può impostare soglie di qualità per la profondità di lettura e le velocità di mappatura) e identifica problemi tecnici come errori di sequenziamento, effetti batch o etichettatura errata. Se l'input dell'RNA è di qualità sufficiente, Rup può aiutare tagliando le letture di bassa qualità o identificando i campioni etichettati in modo errato. Tuttavia, Rup non è in grado di compensare la scarsa qualità dell'RNA in ingresso. In caso di RNA di bassa qualità o dati di sequenziamento errati, potrebbe essere necessario raccogliere nuovamente i campioni, regolare il protocollo di estrazione dell'RNA e/o ripetere il sequenziamento. Lo stesso vale per i campioni con un'ampia frazione di letture non mappate che potrebbero essere state causate da contaminazione. Tuttavia, in base alla nostra esperienza, non tutte le estrazioni di RNA possono essere semplicemente ripetute a causa della mancanza di disponibilità di materiale. In questo caso, non sono adatti per alcune applicazioni a valle: i campioni con un basso numero di reads a mappatura singola non devono essere analizzati nelle analisi di espressione genica differenziale, ma contengono comunque informazioni per l'analisi della presenza del trascritto. In questo caso, l'assenza di trascritti in tessuti/trattamenti/condizioni non può essere considerata per l'analisi e la quantificazione dell'abbondanza di trascritti.

Rup include alcune limitazioni. Ad esempio, non verifica direttamente la degradazione dell'RNA poiché sono necessarie misurazioni dirette dell'integrità dell'RNA prima della preparazione e del sequenziamento della libreria. Tuttavia, uno script per identificare la degradazione dell'RNA utilizzando i moduli RSeQC geneBodyCoverage.py e tin.py è incluso nel repository di questa pipeline. Inoltre, Rup non verifica il contenuto GC e la distorsione della lunghezza della trascrizione. La dimensione del genoma di riferimento è attualmente limitata a 4 Gb e, secondo la nostra esperienza, il modulo di mappatura delle letture sovrastima le letture multi-mappate nei poliploidi, come esemplificato nel file supplementare in un confronto tra la mappatura delle letture con le versioni del genoma aploide rispetto a quelle diploidi di E. californica, un genoma aploide è quindi l'input preferito per questa pipeline. La qualità dell'output Rup dipende in gran parte dalla qualità del genoma di riferimento e dall'annotazione, in modo tale che una sequenza genomica incompleta o altamente frammentata potrebbe portare a una sovrastima delle letture non mappate. Inoltre, l'annotazione incompleta del modello genico porta a una sovrastima delle letture non assegnate. La completezza e la duplicazione della sequenza del genoma e dell'annotazione genica possono essere dedotte con strumenti come BUSCO²⁰.

Rup è uno strumento autonomo per tutte le fasi iniziali del controllo di qualità, essenziali negli esperimenti di RNA-seq. A differenza di RSeQC e RNA-SeQC, non è richiesta la pre-elaborazione delle letture di sequenziamento grezze. L'analisi della qualità del sequenziamento non può essere eseguita con RNA-SeQC e le mappe di calore della correlazione dei campioni non vengono calcolate da RSeQC e RNA-SeQC (Tabella 1). Inoltre, l'output di normalizzazione è in FPKM in RSeQC e la visualizzazione dell'output non è implementata come standard per tutti i moduli in RSeQC e RNA-SeQC. Pertanto, diversi problemi di controllo della qualità non vengono testati nei due strumenti alternativi, che includono un basso numero di letture, un'alta frazione di letture tagliate e l'identificazione di valori anomali replicati. Un confronto tra Rup, RSeQC e RNA-SeQC è disponibile nel File Supplementare 1. Inoltre, Rup può essere utilizzato in modo complementare a RNA-SeQC o RSeQC, ad esempio, i file BAM ordinati prodotti da questa pipeline possono essere utilizzati come input per questi strumenti.

Attualmente, Rup è ottimizzato per un numero limitato di campioni, ma i passaggi più dispendiosi in termini di tempo, come il taglio di qualità e la mappatura di lettura, possono essere precalcolati, ad esempio, su un cluster di computer o su un'infrastruttura cloud. Per i genomi più grandi di 4 Gb, questo è obbligatorio, poiché l'allineatore Rsubread è limitato ai genomi più piccoli di 4 Gb. L'annotazione dell'rRNA viene eseguita con barrnap, che identifica i geni rRNA altamente conservati. Secondo la nostra esperienza, l'annotazione del gene dell'rRNA non richiede completezza, perché anche se mancano alcuni geni anomali dell'rRNA, una panoramica grossolana della presenza di rRNA nel set di dati è sufficiente per la valutazione della qualità del set di dati. Le versioni future di Rup potrebbero passare a un metodo di mappatura diverso, come lo strumento di pseudo-allineamento salmon²¹, per ridurre il tempo di esecuzione. Inoltre, a Rup potrebbero essere aggiunti altri metodi di normalizzazione e correzione, come la distorsione GC o la correzione della distorsione della lunghezza.

In sintesi, Rup fornisce informazioni essenziali per garantire l'affidabilità e la riproducibilità dell'analisi dei dati RNA-seq. Questa pipeline riporta in modo completo le principali metriche di qualità dell'RNA-seq e produce file di output per l'uso diretto nelle analisi a valle. È stato progettato come strumento autonomo per i ricercatori con conoscenze bioinformatiche minime per valutare la qualità dei dati di sequenziamento primario con una visualizzazione intuitiva.

Materiali

Elenco dei materiali utilizzati in questo articolo

Nome	Azienda	Numero di catalogo	Commenti
fastQCR	R	0.1.3
FUJITSU ESPRIMO D958 Desktop	FUJITSU		la pipeline è stata sviluppata e testata su Ubuntu Linux 24.02 (32 Gb RAM)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	Mela		la pipeline è stata testata su maxOS 15.4.1 (16 Gb RAM)
Pheatmap	R	1.0.13
R		4.4.3
rimodellazione2	R	1.4.4
RFASP	Bioconduttore	1.16.0
rsamtools	Bioconduttore	2.22.0
rsubread	Bioconduttore	2.20.0