Rup (RNA-seq Usability Assessment Pipeline) - Control de calidad para experimentos de RNA-seq a granel en eucariotas

Resumen

Este protocolo permite el control de calidad inicial de los experimentos de secuenciación de ARN para biólogos de laboratorio húmedo con experiencia limitada en bioinformática.

Resumen

Los enfoques modernos en la ciencia molecular de las plantas a menudo requieren experimentos de secuenciación de ARN a granel, por ejemplo, para rastrear los cambios globales en los transcriptomas durante los tratamientos o para identificar componentes clave de las vías reguladoras. En consecuencia, diversas áreas de las ciencias de las plantas dependen de datos de secuenciación de ARN a granel de alta calidad y reproducibles para el avance científico. Sin embargo, según nuestra experiencia, a menudo faltan conocimientos y aplicación de medidas de control de calidad en conjuntos de datos de secuenciación de ARN. Aquí, presentamos Rup (RNA-seq usability assessment pipeline) para el control de calidad de los datos de RNA-seq a granel, para análisis posteriores de expresión génica, que es independiente y fácilmente aplicable para biólogos de laboratorio húmedo con conocimientos básicos de R. Rup ayuda a discriminar entre datos de secuenciación de alta calidad, adecuados para experimentos de expresión génica posteriores, y aquellos que no son adecuados para análisis posteriores generales. Rup incluye pruebas para varios problemas comunes, como números de lectura o mapeo insuficientes, identificación de contaminaciones, cuantificación de fracciones de ARNr en los datos totales de secuenciación de ARN, pruebas de similitud de replicación, uso de datos reales para demostración y visualización intuitiva. Rup proporciona un conjunto de herramientas para identificar deficiencias experimentales antes del análisis estandarizado del transcriptoma, mejorando así la calidad de los datos para los investigadores individuales y el campo. Esto mejora la confianza en el análisis de datos de secuenciación de ARN a granel y proporciona una base para futuras directrices que definan criterios mínimos de control de calidad, mejorando así la fiabilidad y transparencia de los datos de secuenciación de ARN publicados. Los datos de Rup y prueba están disponibles en https://github.com/oliverrupp/rup.

Introducción

Los experimentos de transcriptómica (RNA-seq) interrogan exhaustivamente las firmas transcripcionales que dan forma a los fenotipos. Este enfoque se volvió insustituible en la genética molecular de las plantas para identificar genes únicos o coexpresados y procesos biológicos involucrados, por ejemplo, en el desarrollo, la interacción con patógenos de plantas o la resistencia al estrés abiótico^{1, 2 ,3}. Los avances recientes en la tecnología de secuenciación de ARN han aumentado la especificidad y han permitido la detección de diferentes isoformas y variantes con una resolución de base única, lo que permite la identificación de variaciones de secuencia desde indeles más grandes hasta polimorfismos de un solo nucleótido (SNP). Los datos obtenidos por RNA-seq se caracterizan por un amplio rango dinámico, lo que facilita la detección de transcripciones altamente abundantes y de baja expresión, y requieren condiciones experimentales apropiadas para su consistencia. Además, los conjuntos de datos de secuenciación de ARN son grandes y, por lo tanto, intensivos en computación para analizar y almacenar, lo que implica una gestión de datos eficiente y amplios recursos informáticos. RNA-seq requiere una entrada de alta calidad en cada paso del flujo de trabajo, ya que las debilidades en cualquier etapa pueden propagarse y comprometer los resultados. La mala integridad del ARN o los problemas técnicos durante la preparación de la biblioteca conducirán a sesgos y una reducción de la precisión. Los parámetros para obtener lecturas sin procesar de alta calidad pueden variar entre las instalaciones de secuenciación de próxima generación (NGS).

Según nuestra experiencia, recomendamos usar solo ARN con un número de integridad de ARN (RIN) superior a 7, lo que es indicativo de una estructura de ARNm en gran parte intacta, como entrada para la preparación de la biblioteca de secuenciación de ARN. Para una secuenciación exitosa, se requieren aproximadamente 2 μg de ARN total a una concentración de 50-200 ng/μL para los protocolos de preparación de bibliotecas estándar. La pureza del ARN debe confirmarse mediante una relación DO_260/280 entre 1,8 y 2,1 y una relación DO_260/230 superior a 1,5 utilizando un espectrofotómetro. La profundidad de secuenciación insuficiente, las bajas tasas de mapeo o la desalineación con el genoma de referencia pueden distorsionar aún más la expresión génica y la cuantificación del empalme. Además, un diseño experimental inadecuado, como la falta de discriminación entre muestras de diferentes tejidos o tratamientos y la alta variabilidad entre réplicas, puede introducir ruido y disminuir la reproducibilidad. Si bien muchos laboratorios analizan los transcriptomas de forma rutinaria, los estrictos controles de calidad de los primeros pasos esenciales del análisis a menudo no se informan o pueden estar completamente ausentes de las publicaciones. Esto puede llevar a una sobreinterpretación de los resultados derivados del análisis del transcriptoma y, en consecuencia, a resultados irreproducibles.

Aquí, proporcionamos un flujo de trabajo para el control de calidad de los pasos iniciales necesarios para los análisis de transcriptoma de alta calidad de los perfiles de ARNm para medir las alteraciones transcripcionales. Nuestro objetivo es permitir que los biólogos de laboratorio húmedo con conocimientos limitados en bioinformática evalúen sus datos transcriptómicos primarios. Rup (Figura 1) es accesible para los investigadores que están familiarizados con el conocimiento básico de R. La ejecución del flujo de trabajo presentado aquí proporcionará a los investigadores una comprensión detallada de sus datos primarios, incluidas sus posibles limitaciones para el análisis posterior. Hasta donde sabemos, hasta ahora falta una línea práctica de evaluación de datos de transcriptoma primario en combinación con pautas para distinguir datos de alta y baja calidad.

Figura 1: Flujo de trabajo de la canalización de control de calidad RNA-seq in silico . Los archivos de entrada para el control de calidad se derivan de los datos de secuenciación de ARN generados por la secuenciación del material vegetal, así como de conjuntos de datos disponibles públicamente. La canalización de R proporcionada evalúa la calidad de la secuencia a través de tres enfoques principales: calidad de secuenciación, calidad de mapeo y calidad de replicación. Se calculan varias métricas de calidad y los resultados estadísticos se visualizan utilizando paquetes comunes de R como ggplot2 y pheatmap (por ejemplo, gráficos de barras y mapas de calor). Haga clic aquí para ver una versión más grande de esta figura.

Las canalizaciones de evaluación notificadas anteriormente requerían datos preprocesados (por ejemplo, alineación de lectura como archivos .bam), no cubren todas las métricas o ya no se mantienen ^4,5,6. La ventaja del flujo de trabajo presentado aquí radica en su exhaustividad al consolidar los problemas de control de calidad más comunes en una sola tubería. Además, proporcionamos ejemplos de datos de alta calidad adecuados para todas las aplicaciones de análisis posteriores, pero también ejemplos de datos de baja calidad, y discutimos sus limitaciones específicas para un análisis posterior. Varios informes publicados anteriormente describen el propósito de las herramientas de análisis de secuenciación de ARN y proporcionan evaluaciones comparativas de su rendimiento ^7,8. Sin embargo, los estándares de control de calidad y notificación metodológica de RNA-seq no están estandarizados y agravan la reproducibilidad y las interpretaciones biológicamente significativas de los experimentos transcriptómicos.

Rup integra herramientas estándar de alta calidad, análisis de control de calidad y visualización de los resultados. Rup requiere lecturas de secuenciación sin procesar y un genoma anotado como entrada y se ejecuta en Mac OS, Linux y "Windows Subsystem for Linux" (WSL) en sistemas Windows. Integra herramientas para la calidad de la secuenciación y cuantifica el mapeo de lectura a un genoma, incluye la medición del contenido de ARNr en las muestras y ofrece correlación de muestras para evaluar la correlación de replicación. Los datos de la prueba se obtuvieron mediante microdisección láser del tejido meristemático central de dos etapas diferentes de la especie de planta Eschscholzia californica, un protocolo de entrada ultrabaja para la preparación de bibliotecas, y se secuenciaron en Novaseq 6000.

Protocolo

Rup se puede ejecutar como un solo script con un mínimo de archivos de entrada requeridos. La canalización, como se muestra aquí, requiere datos de secuenciación de ARN de extremos emparejados de Illumina. También se deposita en el repositorio de GitHub una canalización ajustada para la secuenciación de un solo extremo con los mismos pasos de análisis. Solo se necesitan la secuencia del genoma como un archivo fasta, el modelo genético y las anotaciones de ARNr como archivos gtf, y las lecturas de secuenciación sin procesar como archivos fastq.gz. Asegúrese de que los archivos de genoma y anotación se proporcionan como genome.fa, annotation.gtf y rRNA.gtf en una carpeta especificada en la variable reference_folder. Cuando todos los archivos de secuenciación como archivos .fq.gz se guardan en la variable read_file_folder, Rup se puede ejecutar de la siguiente manera:

1. Preparativos:

1. Instale todos los paquetes de R necesarios con Bioconductor.

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

2. Configure todos los archivos necesarios.
   1. Cree una carpeta de origen que contenga todos los archivos necesarios. Agregue la secuencia del genoma de referencia en un archivo fasta llamado "reference/genome.fa" a la carpeta de origen. Proporcione la anotación del modelo genético como un archivo gtf: "reference/annotation.gtf", Opcionalmente, proporcione la anotación del gen rRNA como archivo gtf: "reference/rRNA.gtf".
   2. Agregue todas las lecturas de secuenciación como archivos fastq a la carpeta "lecturas". Asegúrese de que todos los nombres de archivo fastq sigan el mismo patrón; _1.fastq.gz y _2.fastq.gz para las lecturas hacia adelante y hacia atrás.

      # 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. Establezca el parámetro de acuerdo con los recursos informáticos.

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

4. Establezca el parámetro según el método de secuenciación.
    

   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 longitud mínima de lectura no debe ser inferior a 25 pb, ya que las lecturas más pequeñas pueden hacer que se rompa la asignación. Los valores más grandes darán como resultado lecturas asignadas de forma más única, pero también más lecturas descartadas durante el recorte, dependiendo de la calidad de la secuencia. Más detalles sobre estos parámetros se describen en los suplementos.

2. Evaluación de la calidad de la secuenciación

NOTA: Este paso creará un gráfico de barras que muestra el número de lecturas antes y después de recortar cada muestra.

1. Ejecute fastqc⁹, que es una herramienta para el control de calidad general de datos de secuencias de alto rendimiento, que proporciona información sobre el número total de lecturas secuenciadas, la longitud de lectura, el contenido medio de GC, la contaminación del adaptador de secuenciación y las secuencias sobrerrepresentadas (por ejemplo, lecturas de ARNr). Analice la calidad de la secuenciación según los informes de calidad base y secuencial. Ejecute fastqc en los archivos de secuenciación sin procesar usando el paquete de R fastqcr¹⁰ y coloque todos los archivos de entrada en el mismo directorio (read_file_folder). Resuma los archivos de salida de fastqc_folder en un objeto R mediante la función qc_aggregate.
   NOTA: Los resultados se guardarán en el directorio fastqc_folder. Se puede acceder a los archivos .html creados por fastqc en la carpeta de salida para obtener estadísticas adicionales.
    

   # 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 todos los prefijos de muestra almacenados en el sample_prefixes vectorial, itere sobre los prefijos y ejecute la herramienta de recorte de calidad fastp en todas las muestras. Recorte para eliminar secuencias adaptadoras, cebadores, bases de baja calidad, secuencias poli-A/T y es ejecutado por Fastp¹¹, que está disponible en el paquete R Rfastp¹². Los archivos fastq recortados se almacenarán en la carpeta 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. Ejecute fastqc nuevamente en las lecturas recortadas.
    

   # 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. Evalúe el recorte recopilando el número de lecturas antes y después del recorte para todas las muestras. Cree un gráfico de barras con los números leídos usando 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: Resultados de secuenciación y recorte. Lea los recuentos antes (rojo) y después de recortar (verde). La muestra s2_r1 ya tiene un recuento de lecturas bajo antes del recorte, mientras que el recorte eliminó una gran fracción de las lecturas de s2_r2 muestra. Todas las demás muestras muestran una pérdida aceptable de lecturas por el recorte. Haga clic aquí para ver una versión más grande de esta figura.

3. Calidad del mapeo

NOTA: Este paso calcula gráficos de barras para discriminar lecturas de mapeo único (por ejemplo, transcripciones de genes que codifican proteínas) de lecturas de mapeo múltiple (por ejemplo, lecturas de ARNr) y lecturas no mapeadas (por ejemplo, contaminaciones).

1. Utilice el paquete Rsubread¹⁴ para asignar las lecturas al genoma de referencia. Primero, construya un índice del archivo fasta del genoma de referencia (genome_fasta_file).
   NOTA: Este paso se realiza una vez para cada genoma de referencia.
    

   # 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. Iterar sobre todas las muestras y alinear las lecturas con el genoma de referencia utilizando la función align() del paquete Rsubread. Esta función crea archivos .bam en la carpeta de salida.
    

   # 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. Cuente las lecturas de cada gen usando la función featureCounts() del paquete Rsubread. Asegúrese de que el archivo de anotación (annotation_file) esté en formato GTF. Count solo lee con una sola coincidencia con el 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. Cuente las lecturas que se asignan a los genes de ARNr para evaluar el contenido de ARNr en las muestras. Permita que se cuenten aquí las lecturas multimapeadas. Especifique los loci del gen rRNA en un archivo 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. Recopile las estadísticas de asignación de lectura creadas por featureCounts. Las estadísticas incluyen el número de alineaciones para cada muestra en las diferentes categorías (por ejemplo, Asignado, Sin mapear, Multimapeado).
    

   # 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. Recopile las estadísticas de la asignación de ARNr.
    

   # 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. Trace las estadísticas de asignación de lectura como un gráfico de barras.
    

   # 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. Clasifique los genes en grupos según el número de lecturas que se les asignen y grafique los resultados como un gráfico de barras.
    

   # 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: Lea la descripción general y las estadísticas de mapeo. La mayoría de las muestras muestran un alto número de lecturas asignadas (lado izquierdo, marrón) y un bajo número de lecturas de ARNr (lado derecho, azul). La muestra s2_r3 tiene una cantidad inusualmente alta de lecturas de múltiples mapas (lado izquierdo, verde) correspondiente a un número alto de lectura de ARNr (lado derecho). La muestra s2_r4 muestra un alto número de lecturas no mapeadas en combinación con un número de lectura de ARNr esperado, lo que sugiere contaminación con lecturas de otro organismo. Haga clic aquí para ver una versión más grande de esta figura.

4. Replicar la calidad.

1. Asegúrese de que las réplicas de la misma muestra o tejido muestren una mayor correlación entre sí que con las réplicas de diferentes muestras o tejidos. Trace un mapa de calor agrupado de las correlaciones replicadas para visualizar los datos utilizando los recuentos de lectura sin procesar o los recuentos normalizados de transcripciones por millón (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")

Figura 4: Clasificación del recuento de lectura de genes. Todos los genes de una muestra se clasifican en una de cinco categorías según el número de lecturas asignadas. En rojo se muestra el número de genes sin lecturas asignadas. Las muestras con un número bajo de lecturas asignadas totales (s2_r1 a s1_r4) tienen un mayor número de genes con 10 a 100 lecturas asignadas y un menor número de genes con más de 1000 lecturas. Los genes con baja expresión pueden pasar desapercibidos en estas muestras. Haga clic aquí para ver una versión más grande de esta figura.

Resultados

La canalización está completamente implementada como un script R y se ha probado en los sistemas operativos Linux y Mac OS. Los usuarios de Windows pueden usar el Subsistema de Windows para Linux (WSL). El código y los datos de prueba están disponibles como repositorio de GitHub: https://github.com/oliverrupp/rup. Los datos de secuenciación están disponibles en el marco del proyecto ENA EBI PRJEB96400.

Se crearon diez muestras artificialmente a partir de dos muestras reales para ejemplificar diversos problemas que pueden encontrarse durante el control de calidad de la secuenciación de ARN a granel utilizando Rup para el análisis. La muestra s2_r1 fue diseñada para exhibir un número de lectura total bajo, s2_r2 de muestra contiene una gran fracción de lecturas de baja calidad para ser descartadas por el proceso de recorte. La s2_r3 de muestra incluye una gran fracción de lecturas de ARNr y s2_r4 de muestra incluye lecturas de contaminantes que no se asignaron al genoma de referencia. Los nombres de las muestras s1_r5 y s2_r5 se intercambiaron para ilustrar una baja correlación de replicación.

La sección 2 del protocolo permite la identificación de muestras con números de lectura bajos antes o después del recorte (Figura 2). El gráfico de barras muestra el número de lectura más bajo en la muestra s1_r1 antes y después del recorte. Aquí, el número de lectura inicial fue bajo. El bajo número de lecturas de s1_r2 de muestra después del recorte sugiere una gran cantidad de secuencias de adaptadores, errores de secuenciación, secuencias de cebadores, tramos de poli-A/T secuenciados que se eliminaron durante el proceso de recorte. La degradación del ARN de entrada también puede reducir el número de lecturas de alta calidad.

La sección 3 del protocolo identifica problemas en las asignaciones de lectura. La Figura 3 muestra el elevado número de lecturas multimapeadas en la muestra s2_r3 (verde), así como el alto número de lecturas de ARNr. La contaminación de la s2_r4 de la muestra es evidente por la gran fracción de lecturas no asignadas al genoma de referencia (rosa). Estas lecturas no se correlacionan con el ARNr, sino con secuencias de un organismo no objetivo. La Figura 4 muestra problemas generales asociados con un bajo número de lecturas asignadas en transcriptomas. En muestras con una baja tasa de lecturas asignadas únicamente al genoma de referencia (s2_r1 a r4), solo alrededor de un tercio de los genes tienen más de 100 lecturas asignadas, mientras que en las otras muestras, aproximadamente la mitad de los genes caen en estas clases. Es posible que no se encuentren lecturas de genes con muy baja expresión en las muestras s2_r1 a r4 y, en consecuencia, el análisis comparativo de expresión con estas muestras será muy poco confiable y debe evitarse.

La sección 4 del protocolo se puede utilizar para identificar valores atípicos replicados. Se espera que las réplicas de la misma muestra/condición/tejido muestren una mayor correlación entre sí que las réplicas de otras muestras/condiciones/tejidos. La Figura 5 muestra un mapa de calor de correlación de las dos muestras (S1 y S2) con cinco réplicas cada una. Los dendrogramas en la parte superior y lateral del gráfico muestran dos grupos con cinco réplicas en cada grupo. El clúster izquierdo contiene cuatro réplicas de la muestra 1 y una réplica de la muestra 2 (s2_r5), el clúster derecho contiene cuatro réplicas de la muestra 2 y una réplica de la muestra 1 (s1_r5). En este caso, cuando los nombres de muestra s1_r5 y s2_r5 se vuelven a intercambiar, cada clúster contiene todas las réplicas de una muestra, lo que puede indicar un error de etiquetado de réplica. Otras razones por las que las réplicas no se agrupan podrían ser la falta de diferenciación entre las muestras/condiciones/tejidos, o la agrupación de réplicas solo por razones de calidad de secuenciación. Esto último puede ocurrir cuando todas las réplicas con recuentos de lectura excepcionalmente bajos o excepcionalmente altos después del recorte y el mapeo forman un grupo.

Figura 5: Mapa de calor de correlación de muestra. El mapa de calor de correlación de muestra se basa en los valores de TPM transformados de log₂ y muestra dos clústeres distintos de cinco muestras cada uno. Se espera que las réplicas biológicas/técnicas muestren una mayor correlación entre sí que las de otros tejidos/tratamientos. El clúster izquierdo contiene cuatro réplicas de la muestra 1 y una réplica de la muestra 2 (s2_r5), el clúster derecho contiene cuatro réplicas de la muestra 2 y una réplica de la muestra 1 (s1_r5). Las réplicas individuales deben verificarse para detectar posibles intercambios de muestras o efectos por lotes para explicar su agrupación en el mapa de calor. Haga clic aquí para ver una versión más grande de esta figura.

Tabla 1: Comparación de tuberías de evaluación de la calidad de RNA-seq. Para un análisis detallado, consulte el Archivo complementario 1. Haga clic aquí para descargar esta tabla.

Archivo complementario 1: Selección de parámetros y referencias. Los análisis muestran la influencia de la selección de parámetros y referencias en los resultados generales del control de calidad. Haga clic aquí para descargar este archivo.

Discusión

La calidad de los análisis de expresión génica diferencial depende en gran medida de dos factores: el número de lecturas secuenciadas en cada muestra y replicada¹⁶ y el número de réplicas por muestra^17,18. Aquí, presentamos la tubería Rup fácil de usar para determinar el número de lecturas adecuadas para la cuantificación de la expresión génica en cada réplica. Diferentes métricas permiten a los investigadores comprender por qué las réplicas muestran números de lectura asignados bajos e identificar problemas en la correlación de la muestra. Aunque Rup se desarrolló para la evaluación de la calidad de las muestras de secuenciación de ARN de las plantas, es igualmente adecuado para otros organismos eucariotas, ya que no requiere más ajustes para muestras que no son de plantas (Archivo complementario 1).

El primer paso de Rup determina el recuento total de lecturas antes y después del recorte de lecturas. El número de lecturas secuenciadas determina el número de genes detectables, y si el número de lecturas es demasiado bajo, muchos genes expresados diferencialmente permanecerán sin detectar. Una gran fracción de las lecturas descartadas durante el proceso de recorte y filtrado de calidad podría indicar la degradación del ARN, lo que podría conducir a una subestimación de la expresión de genes altamente degradados. Sin embargo, si una muestra se puede utilizar para aplicaciones posteriores depende del organismo objetivo y del objetivo del estudio. Por ejemplo, para capturar la expresión de la mayoría de los genes de las plantas, en nuestra experiencia, se requieren de 30 a 50 millones de lecturas, mientras que para los hongos, solo 10 millones podrían ser suficientes. Por lo tanto, Rup no definirá umbrales de calidad para la exclusión de muestras problemáticas, sino que proporcionará métricas para identificar diversos problemas que puedan encontrarse.

El segundo paso (calidad de mapeo) evalúa la precisión y confiabilidad de la alineación de lectura con el genoma de referencia, especificando su idoneidad para análisis posteriores. No todas las lecturas secuenciadas se pueden usar para el cálculo de la abundancia de genes; Las lecturas que no se alinean con el genoma o el transcriptoma o las lecturas que se asignan a múltiples ubicaciones en el genoma se ignoran en la mayoría de los casos¹⁹ y no contribuyen a los recuentos generales de lecturas. Los resultados del segundo paso de la canalización se pueden usar para comprender por qué las lecturas no se utilizan en el cálculo de abundancia. Un alto número de lecturas no mapeadas podría indicar contaminación durante la extracción de ARN (por ejemplo, con patógenos de plantas o herbívoros) o un genoma de referencia incompleto. Los modelos genéticos incompletos del genoma de referencia pueden dar lugar a un alto número de lecturas "sin característica". Una gran fracción de las lecturas de múltiples mapas podría deberse a un alto número de genes de ARNr en la biblioteca que indican una eliminación insuficiente de ARNr durante el proceso de preparación de la biblioteca de secuenciación (en el caso de que las lecturas de múltiples mapas se deriven realmente de los ARNr, se pueden analizar más a fondo proporcionando un archivo de anotación de ARNr a la canalización, que luego calcula automáticamente el número de posibles lecturas de ARNr).

La tercera métrica de calidad, igualmente importante, es la correlación entre réplicas de la misma muestra. Generalmente, la correlación entre réplicas de la misma muestra debe ser mayor que la correlación entre réplicas de diferentes muestras. Una baja correlación entre las réplicas podría indicar una gran variabilidad biológica entre las réplicas, una alta similitud entre las muestras/condiciones/tejidos, otros efectos de lote o incluso intercambio de muestras o etiquetado incorrecto. Rup calcula correlaciones por pares entre todas las muestras y produce un mapa de calor agrupado de las correlaciones de la muestra. Además, se puede calcular un análisis de componentes principales (PCA) en las muestras para identificar posibles efectos de lotes que deben reconocerse en análisis posteriores. Aunque los efectos por lotes se pueden corregir en los análisis posteriores, podría ser mejor eliminar las muestras problemáticas si la diferencia medida es demasiado alta.

El material de entrada tiene un impacto directo en la calidad de la secuenciación: el muestreo de tejidos, la extracción de ARN y la preparación de bibliotecas son pasos críticos para minimizar los problemas de calidad de la secuenciación de ARN. Se deben mantener condiciones consistentes, si es posible, por ejemplo, mediante el uso de cámaras de crecimiento. Las medidas de control de plagas deben llevarse a cabo de manera oportuna, ya que solo se deben seleccionar individuos sanos para el muestreo. Además, es aconsejable recolectar muestras el mismo día y / o al mismo tiempo para reducir la variación del transcriptoma circadiano. Hay numerosos kits de extracción de ARN disponibles, y la selección de un kit apropiado para la especie objetivo puede mejorar la calidad del ARNm. Los protocolos de preparación de la biblioteca deben incluir pasos para el enriquecimiento de ARN poliA+ para minimizar la fracción de ARNr.

Rup se puede utilizar como un paso inicial de control de calidad en el análisis de expresión génica diferencial. Este control de calidad es necesario por varias razones críticas, ya que garantiza la calidad de los datos de entrada (identifica réplicas/muestras de baja calidad, puede establecer umbrales de calidad para la profundidad de lectura y las tasas de mapeo) e identifica problemas técnicos como errores de secuenciación, efectos de lote o etiquetado incorrecto. Si la entrada de ARN es de calidad suficiente, Rup puede ayudar recortando lecturas de baja calidad o identificando muestras mal etiquetadas. Sin embargo, Rup no puede compensar la mala calidad del ARN de entrada. En casos de ARN de baja calidad o datos de secuenciación erróneos, puede ser necesario volver a recolectar muestras, ajustar el protocolo de extracción de ARN y/o repetir la secuenciación. Lo mismo se aplica a las muestras con una gran fracción de lecturas no mapeadas que pueden haber sido causadas por la contaminación. Sin embargo, según nuestra experiencia, no todas las extracciones de ARN pueden repetirse simplemente debido a la falta de disponibilidad de material. En este caso, no son adecuados para algunas aplicaciones posteriores: las muestras con un bajo número de lecturas de un solo mapa no deben analizarse en análisis de expresión génica diferencial, pero aún contienen información para el análisis de la presencia de transcripciones. En este caso, la ausencia de transcripciones en tejidos/tratamientos/condiciones no puede considerarse para el análisis y cuantificación de la abundancia de transcripciones.

Rup incluye algunas limitaciones. Por ejemplo, no prueba directamente la degradación del ARN, ya que se requieren mediciones directas de la integridad del ARN antes de la preparación y secuenciación de la biblioteca. Sin embargo, en el repositorio de esta canalización se incluye un script para identificar la degradación del ARN mediante los módulos RSeQC geneBodyCoverage.py y tin.py. Además, Rup no prueba el contenido de GC y el sesgo de longitud de la transcripción. El tamaño del genoma de referencia está actualmente limitado a 4 Gb y, según nuestra experiencia, el módulo de mapeo de lectura sobreestima las lecturas de múltiples mapas en poliploides, como se ejemplifica en el archivo complementario en una comparación del mapeo de lectura con un haploide frente a las versiones del genoma diploide de E. californica, un genoma haploide es, por lo tanto, la entrada preferida para esta canalización. La calidad de salida de Rup depende en gran medida de la calidad del genoma de referencia y la anotación, de modo que una secuencia del genoma incompleta o altamente fragmentada podría conducir a una sobreestimación de lecturas no mapeadas. Además, la anotación incompleta del modelo genético conduce a una sobreestimación de lecturas no asignadas. La integridad y duplicación de la secuencia del genoma y la anotación de genes se pueden inferir con herramientas como BUSCO²⁰.

Rup es una herramienta independiente para todos los pasos iniciales de control de calidad esenciales en los experimentos de secuenciación de ARN. A diferencia de RSeQC y RNA-SeQC, no se requiere preprocesamiento de las lecturas de secuenciación sin procesar. El análisis de calidad de secuenciación no se puede realizar con RNA-SeQC, y los mapas de calor de correlación de muestras no se calculan con RSeQC y RNA-SeQC (Tabla 1). Además, la salida de normalización está en FPKM en RSeQC y la visualización de salida no se implementa como estándar para todos los módulos en RSeQC y RNA-SeQC. Por lo tanto, varios problemas de control de calidad no se prueban en las dos herramientas alternativas, que incluyen un número de lectura bajo, una fracción alta de lecturas recortadas y la identificación de valores atípicos replicados. Una comparación de Rup, RSeQC y RNA-SeQC está disponible en el Archivo complementario 1. Además, Rup se puede usar de forma complementaria a RNA-SeQC o RSeQC, por ejemplo, los archivos BAM ordenados producidos por esta tubería se pueden usar como entrada para estas herramientas.

Actualmente, Rup está optimizado para un pequeño número de muestras, pero los pasos que consumen más tiempo, como el recorte de calidad y el mapeo de lectura, se pueden precalcular, por ejemplo, en un clúster de computadoras o una infraestructura en la nube. Para genomas mayores de 4 Gb, esto es obligatorio, ya que el alineador Rsubread está limitado a genomas menores de 4 Gb. La anotación de ARNr se realiza con barrnap, que identifica los genes de ARNr altamente conservados. Según nuestra experiencia, la anotación de genes de ARNr no requiere exhaustividad, porque incluso si faltan algunos genes de ARNr anormales, una descripción general de la presencia de ARNr en el conjunto de datos es suficiente para la evaluación de la calidad del conjunto de datos. Las versiones futuras de Rup podrían cambiar a un método de mapeo diferente, como la herramienta de pseudoalineación salmon²¹, para disminuir el tiempo de ejecución. Además, se podrían agregar más métodos de normalización y corrección, como el sesgo GC o la corrección del sesgo de longitud, a Rup.

En resumen, Rup proporciona información esencial para garantizar la fiabilidad y reproducibilidad del análisis de datos de RNA-seq. Esta canalización informa de forma exhaustiva de las principales métricas de calidad de RNA-seq y produce archivos de salida para su uso directo en los análisis posteriores. Fue diseñado como una herramienta independiente para que los investigadores con conocimientos mínimos de bioinformática evalúen la calidad de sus datos de secuenciación primaria con una visualización intuitiva.

Materiales

Lista de materiales utilizados en este artículo

Nombre	Empresa	Número de catálogo	Comentarios
fastqcr	R	0.1.3
FUJITSU ESPRIMO D958 Sobremesa	FUJITSU		la cadena se desarrolló y probó en Ubuntu Linux 24.02 (32 Gb de RAM)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro & nbsp;	Manzana		la tubería se probó en maxOS 15.4.1 (16 Gb de RAM)
pheatmap	R	1.0.13
R		4.4.3
Reshape2	R	1.4.4
RFASP	Bioconductor	1.16.0
rsamtools	Bioconductor	2.22.0
rsubread	Bioconductor	2.20.0