Rup (RNA-seq Usability Assessment Pipeline) - Contrôle de la qualité pour les expériences de RNA-seq en vrac chez les eucaryotes

Résumé

Ce protocole permet un contrôle initial de la qualité des expériences de séquençage de l’ARN pour les biologistes de laboratoire ayant une expérience limitée en bioinformatique.

Résumé

Les approches modernes en science moléculaire des plantes nécessitent souvent des expériences de séquençage de l’ARN en vrac, par exemple, pour suivre les changements globaux des transcriptomes lors des traitements ou pour identifier les composants clés des voies de régulation. Par conséquent, divers domaines des sciences végétales s’appuient sur des données de séquençage de l’ARN en vrac de haute qualité et reproductibles pour le progrès scientifique. Cependant, d’après notre expérience, la connaissance et l’application des mesures de contrôle de la qualité dans les ensembles de données de séquençage de l’ARN font souvent défaut. Ici, nous présentons Rup (RNA-seq usability assessment pipeline) pour le contrôle de la qualité des données de séquençage d’ARN en vrac, pour les analyses ultérieures de l’expression génique, qui est autonome et facilement applicable pour les biologistes de laboratoire humide ayant une connaissance de base de R. Rup aide à faire la distinction entre les données de séquençage de haute qualité, adaptées aux expériences d’expression génique en aval, et celles qui ne conviennent pas à une analyse générale plus approfondie. Rup comprend des tests pour plusieurs problèmes couramment rencontrés tels que le nombre insuffisant de lectures ou de cartographie, l’identification des contaminations, la quantification des fractions d’ARNr dans les données totales de séquençage de l’ARN, les tests de similarité répliquée, l’utilisation de données réelles pour la démonstration et la visualisation intuitive. Rup fournit une suite d’outils pour identifier les lacunes expérimentales avant l’analyse standardisée du transcriptome, améliorant ainsi la qualité des données pour les chercheurs individuels et le terrain. Cela renforce la confiance dans l’analyse des données de RNA-seq en vrac et fournit une base pour les futures directives définissant des critères de contrôle de qualité minimum, améliorant ainsi la fiabilité et la transparence des données RNA-seq publiées. Les données sur le rup et les tests sont disponibles à l’adresse https://github.com/oliverrupp/rup.

Introduction

Les expériences de transcriptomique (RNA-seq) interrogent de manière exhaustive les signatures transcriptionnelles qui façonnent les phénotypes. Cette approche est devenue irremplaçable en génétique moléculaire végétale pour identifier les gènes uniques ou co-exprimés et les processus biologiques impliqués, par exemple, dans le développement, l’interaction entre les agents pathogènes des plantes ou la résistance au stress abiotique^{1, 2 ,3}. Les progrès récents de la technologie RNA-seq ont augmenté la spécificité et permis la détection de différentes isoformes et variantes à une résolution de base unique, permettant l’identification de variations de séquence allant de plus grands indels à des polymorphismes mononucléotidiques (SNP). Les données obtenues par RNA-seq sont caractérisées par une large gamme dynamique, facilitant la détection de transcrits très abondants et faiblement exprimés, et nécessitent des conditions expérimentales appropriées pour être cohérentes. De plus, les ensembles de données de séquençage de l’ARN sont volumineux et nécessitent donc beaucoup de calculs à analyser et à stocker, ce qui implique une gestion efficace des données et des ressources informatiques étendues. Le séquençage de l’ARN nécessite une entrée de haute qualité à chaque étape du flux de travail, car les faiblesses à chaque étape peuvent se propager et compromettre les résultats. Une mauvaise intégrité de l’ARN ou des problèmes techniques lors de la préparation de la bibliothèque entraîneront un biais et une réduction de la précision. Les paramètres permettant d’obtenir des lectures brutes de haute qualité peuvent varier d’une installation de séquençage de nouvelle génération (NGS).

D’après notre expérience, nous recommandons d’utiliser uniquement de l’ARN dont l’indice d’intégrité de l’ARN (RIN) est supérieur à 7, ce qui indique une structure d’ARNm largement intacte, comme entrée pour la préparation de la banque de séquençage d’ARN. Pour un séquençage réussi, environ 2 μg d’ARN total à une concentration de 50 à 200 ng/μL sont nécessaires pour les protocoles standard de préparation de banques. La pureté de l’ARN doit être confirmée par un rapport DO_260/280 compris entre 1,8 et 2,1 et un rapport DO_260/230 supérieur à 1,5 à l’aide d’un spectrophotomètre. Une profondeur de séquençage insuffisante, de faibles taux de cartographie ou un désalignement par rapport au génome de référence peuvent fausser davantage l’expression des gènes et la quantification de l’épissage. De plus, une conception expérimentale inadéquate, telle qu’un manque de discrimination entre les échantillons de différents tissus ou traitements et une grande variabilité entre les répétitions, peut introduire du bruit et diminuer la reproductibilité. Bien que de nombreux laboratoires analysent régulièrement les transcriptomes, il arrive souvent que les contrôles de qualité rigoureux des premières étapes d’analyse essentielles ne soient pas signalés ou soient totalement absents des publications. Cela peut conduire à une surinterprétation des résultats dérivés de l’analyse du transcriptome et, par conséquent, à des résultats non reproductibles.

Ici, nous fournissons un flux de travail pour le contrôle de la qualité des étapes initiales requises pour des analyses transcriptologiques de haute qualité des profils d’ARNm afin de mesurer les altérations transcriptionnelles. Notre objectif est de permettre aux biologistes de laboratoire ayant des connaissances limitées en bioinformatique d’évaluer leurs données transcriptomiques primaires. Rup (Figure 1) est accessible aux chercheurs qui connaissent les connaissances de base de R. L’exécution du flux de travail présenté ici fournira aux chercheurs une compréhension détaillée de leurs données primaires, y compris leurs limites potentielles pour une analyse ultérieure. À notre connaissance, il n’existe pas jusqu’à présent de pipeline pratique d’évaluation des données de transcriptome primaire, associé à des lignes directrices permettant de distinguer les données de haute qualité des données de faible qualité.

Figure 1 : Flux de travail du pipeline de contrôle qualité in silico RNA-seq. Les fichiers d’entrée pour le contrôle de la qualité sont dérivés des données de séquençage de l’ARN générées par le séquençage du matériel végétal, ainsi que des ensembles de données accessibles au public. Le pipeline R fourni évalue la qualité de la séquence à travers trois approches principales : la qualité du séquençage, la qualité de la cartographie et la qualité de la réplication. Diverses mesures de qualité sont calculées et les résultats statistiques sont visualisés à l’aide de packages R courants tels que ggplot2 et pheatmap (par exemple, des graphiques à barres et des cartes thermiques). Veuillez cliquer ici pour voir une version agrandie de cette figure.

Les pipelines d’évaluation précédemment signalés nécessitaient des données prétraitées (par exemple, l’alignement de lecture sous forme de fichiers .bam), ne couvrent pas toutes les métriques ou ne sont plus gérés ^4,5,6. L’avantage du flux de travail présenté ici réside dans son exhaustivité en regroupant les problèmes de contrôle qualité les plus courants dans un seul pipeline. De plus, nous fournissons des exemples de données de haute qualité adaptées à toutes les applications d’analyse en aval, mais aussi des exemples de données de faible qualité, et discutons de leurs limites spécifiques pour une analyse plus approfondie. Plusieurs rapports publiés précédemment décrivent l’objectif des outils d’analyse RNA-seq et fournissent des évaluations comparatives de leurs performances ^7,8. Cependant, les normes de contrôle de la qualité du séquençage de l’ARN et de rapport méthodologique ne sont pas normalisées et aggravent la reproductibilité et les interprétations biologiquement significatives des expériences de transcriptomique.

Rup intègre des outils standard de haute qualité, l’analyse du contrôle de la qualité et la visualisation des résultats. Rup nécessite des lectures de séquençage brutes et un génome annoté en entrée et s’exécute sur Mac OS, Linux et le « sous-système Windows pour Linux » (WSL) sur les systèmes Windows. Il intègre des outils pour la qualité du séquençage et quantifie la cartographie de lecture d’un génome, inclut la mesure du contenu en ARNr dans les échantillons et propose une corrélation d’échantillons pour évaluer la corrélation des réplications. Les données de test ont été obtenues à l’aide de la microdissection laser du tissu du méristème central de deux stades différents de l’espèce végétale Eschscholzia californica, un protocole d’entrée ultra-faible pour la préparation de la banque, et séquencées sur Novaseq 6000.

Protocole

Rup peut être exécuté comme un seul script avec un minimum de fichiers d’entrée requis. Le pipeline, comme illustré ici, nécessite des données de séquençage de l’ARN à extrémité appariée Illumina. Un pipeline ajusté pour le séquençage à une extrémité avec les mêmes étapes d’analyse est également déposé dans le référentiel GitHub. Seules la séquence du génome sous forme de fichier fasta, le modèle de gène et les annotations d’ARNr sous forme de fichiers gtf, et les lectures de séquençage brutes sous forme de fichiers fastq.gz sont nécessaires. Assurez-vous que les fichiers de génome et d’annotation sont fournis sous la forme genome.fa, annotation.gtf et rRNA.gtf dans un dossier spécifié dans la variable reference_folder. Lorsque tous les fichiers de séquençage en tant que fichiers .fq.gz sont enregistrés dans la variable read_file_folder, Rup peut être exécuté comme suit :

1. Préparatifs :

1. Installez tous les packages R nécessaires à l’aide de Bioconductor.

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

2. Configurez tous les fichiers nécessaires.
   1. Créez un dossier source qui contiendra tous les fichiers nécessaires. Ajoutez la séquence génomique de référence dans le fichier fasta appelé « reference/genome.fa » au dossier source. Fournissez l’annotation du modèle de gène sous la forme d’un fichier gtf : « reference/annotation.gtf », si vous le souhaitez, fournissez l’annotation du gène de l’ARNr sous la forme d’un fichier gtf : « reference/rRNA.gtf ».
   2. Ajoutez toutes les lectures de séquençage sous forme de fichiers fastq dans le dossier « reads ». Assurez-vous que tous les noms de fichiers fastq suivent le même modèle ; _1.fastq.gz et _2.fastq.gz pour les lectures aller et arrière.

      # 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. Définissez le paramètre en fonction des ressources de calcul.

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

4. Réglez le paramètre en fonction de la méthode de séquençage.
    

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

   
   
   REMARQUE : La longueur de lecture minimale ne doit pas être inférieure à 25 pb, car des lectures plus petites peuvent entraîner l’interruption du mappage. Des valeurs plus élevées se traduiront par des lectures mappées de manière plus unique, mais aussi par un plus grand nombre de lectures ignorées lors du découpage, en fonction de la qualité du séquençage. Plus de détails sur ces paramètres sont décrits dans les suppléments.

2. Évaluation de la qualité du séquençage

REMARQUE : Cette étape créera un graphique à barres indiquant le nombre de lectures avant et après le découpage de chaque échantillon.

1. Exécutez fastqc⁹, qui est un outil de contrôle de qualité général des données de séquence à haut débit, fournissant des informations sur le nombre total de lectures séquencées, la longueur de lecture, le contenu GC moyen, la contamination de l’adaptateur de séquençage et les séquences surreprésentées (par exemple, les lectures d’ARNr). Analysez la qualité du séquençage selon les rapports de qualité de base et par séquence. Exécutez fastqc sur les fichiers de séquençage bruts à l’aide du package R fastqcr¹⁰ et placez tous les fichiers d’entrée dans le même répertoire (read_file_folder). Synthétisez les fichiers de sortie de fastqc_folder en un objet R à l’aide de la fonction qc_aggregate.
   REMARQUE : Les résultats seront enregistrés dans le répertoire fastqc_folder. Les fichiers .html créés par fastqc dans le dossier de sortie sont accessibles pour des statistiques supplémentaires.
    

   # 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. Avec tous les préfixes d’échantillon stockés dans le sample_prefixes vectoriel, itérez sur les préfixes et exécutez l’outil de découpage de qualité fastp sur tous les échantillons. Le découpage pour supprimer les séquences d’adaptateur, les amorces, les bases de mauvaise qualité, les séquences poly-A/T est exécuté par Fastp¹¹, qui est disponible dans le package R Rfastp¹². Les fichiers fastq tronqués seront stockés dans le dossier 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. Exécutez à nouveau fastqc sur les lectures tronquées.
    

   # 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. Évaluez le découpage en recueillant le nombre de lectures avant et après le parage pour tous les échantillons. Créez un graphique à barres avec les nombres lus à l’aide de 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)

Figure 2 : Résultats du séquençage et du découpage. La lecture compte avant (rouge) et après la coupe (vert). Le s2_r1 d’échantillon a déjà un faible nombre de lectures avant le découpage, tandis que le découpage a supprimé une grande partie des lectures de s2_r2 d’échantillon. Tous les autres échantillons montrent une perte de lecture acceptable lors du rognage. Veuillez cliquer ici pour voir une version agrandie de cette figure.

3. Qualité de la cartographie

REMARQUE : Cette étape calcule des diagrammes à barres pour distinguer les lectures à cartographie unique (par exemple, les transcrits de gènes codant pour des protéines) des lectures à plusieurs cartes (par exemple, les lectures d’ARNr) et des lectures non cartographiées (par exemple, les contaminations).

1. Utilisez le package Rsubread¹⁴ pour mapper les lectures au génome de référence. Tout d’abord, construisez un index du fichier fasta du génome de référence (genome_fasta_file).
   REMARQUE : Cette étape est effectuée une fois pour chaque génome de référence.
    

   # 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. Itérez sur tous les échantillons et alignez les lectures sur le génome de référence à l’aide de la fonction align() du package Rsubread. Cette fonction crée des fichiers .bam dans le dossier de sortie.
    

   # 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. Comptez les lectures pour chaque gène à l’aide de la fonction featureCounts() du package Rsubread. Assurez-vous que le fichier d’annotation (annotation_file) est au format GTF. Le comptage n’est lu qu’avec une seule correspondance avec le génome.
    

   # 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. Comptez les lectures correspondant aux gènes de l’ARNr pour évaluer le contenu de l’ARNr dans les échantillons. Permet de compter les lectures multimappées ici. Spécifiez les loci du gène de l’ARNr dans un fichier 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. Collectez les statistiques d’attribution de lecture créées par featureCounts. Les statistiques incluent le nombre d’alignements pour chaque échantillon dans les différentes catégories (par exemple, Attribué, Non Mappé, MultiMappé).
    

   # 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. Collecter les statistiques d’attribution de l’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. Tracez les statistiques de mappage de lecture sous forme de graphique à barres.
    

   # 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. Classez les gènes en groupes en fonction du nombre de lectures qui leur sont attribuées et tracez les résultats sous forme de graphique à barres.
    

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

Figure 3 : Lire la vue d’ensemble et les statistiques de la cartographie. La plupart des échantillons montrent un nombre élevé de lectures attribuées (côté gauche, marron) et un faible nombre de lectures d’ARNr (côté droit, bleu). L’échantillon s2_r3 présente un nombre inhabituellement élevé de lectures multi-cartographiées (côté gauche, vert) correspondant à un nombre élevé de lectures d’ARNr (côté droit). L’échantillon s2_r4 montre un nombre élevé de lectures non cartographiées en combinaison avec un nombre attendu de lectures d’ARNr, suggérant une contamination par des lectures provenant d’un autre organisme. Veuillez cliquer ici pour voir une version agrandie de cette figure.

4. Reproduire la qualité.

1. Assurez-vous que les répétitions du même échantillon ou tissu présentent une corrélation plus élevée entre elles qu’avec les répétitions provenant d’échantillons ou de tissus différents. Tracez une carte thermique en cluster des corrélations de réplication pour visualiser les données à l’aide du nombre de lectures brutes ou du nombre normalisé de transcriptions par million (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")

Figure 4 : Classification du nombre de gènes lus. Tous les gènes d’un échantillon sont classés dans l’une des cinq catégories suivantes en fonction du nombre de lectures attribuées. Le nombre de gènes sans aucune lecture attribuée est indiqué en rouge. Les échantillons avec un faible nombre total de lectures attribuées (s2_r1 à s1_r4) ont un nombre plus élevé de gènes avec 10 à 100 lectures attribuées et un nombre inférieur de gènes avec plus de 1000 lectures. Les gènes à faible expression peuvent passer inaperçus dans ces échantillons. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Résultats

Le pipeline est entièrement implémenté sous la forme d’un script R et a été testé sur les systèmes d’exploitation Linux et Mac OS. Les utilisateurs de Windows peuvent utiliser le sous-système Windows pour Linux (WSL). Le code et les données de test sont disponibles sous la forme d’un dépôt GitHub : https://github.com/oliverrupp/rup. Les données de séquençage sont disponibles dans le cadre du projet ENA EBI PRJEB96400.

Dix échantillons ont été créés artificiellement à partir de deux échantillons réels pour illustrer divers problèmes qui peuvent être rencontrés lors du contrôle de la qualité du séquençage de l’ARN en vrac à l’aide de Rup pour l’analyse. L’échantillon s2_r1 a été conçu pour présenter un faible nombre total de lectures, s2_r2 échantillon contient une grande fraction de lectures de mauvaise qualité à éliminer par le processus de découpage. La s2_r3 de l’échantillon comprend une grande partie des lectures d’ARNr et la s2_r4 de l’échantillon comprend des lectures de contaminants qui n’ont pas réussi à correspondre au génome de référence. Les noms des échantillons s1_r5 et s2_r5 ont été échangés pour illustrer la faible corrélation de réplication.

La section 2 du protocole permet d’identifier les échantillons dont le nombre de lectures est faible avant ou après le parage (figure 2). Le graphique à barres montre le nombre de lectures le plus bas dans l’échantillon s1_r1 avant et après le découpage. Ici, le nombre initial de lectures était faible. Le faible nombre de lectures d’échantillons s1_r2 après le découpage suggère un grand nombre de séquences d’adaptateur, d’erreurs de séquençage, de séquences d’amorces, d’étirements poly-A/T séquencés qui ont été supprimés pendant le processus de découpage. La dégradation de l’ARN d’entrée peut également réduire le nombre de lectures de haute qualité.

La section 3 du protocole identifie les problèmes dans les affectations de lecture. La figure 3 illustre le nombre élevé de lectures multi-cartographiées dans l’échantillon s2_r3 (vert), ainsi que le nombre élevé de lectures d’ARNr. La contamination de l’échantillon s2_r4 est apparente par la grande fraction de lectures non cartographiées sur le génome de référence (rose). Ces lectures ne sont pas corrélées à l’ARNr, mais aux séquences d’un organisme non cible. La figure 4 montre les problèmes généraux associés au faible nombre de lectures attribuées dans les transcriptomes. Dans les échantillons avec un faible taux de lectures cartographiées uniquement au génome de référence (s2_r1 à r4), seulement environ un tiers des gènes ont plus de 100 lectures attribuées alors que dans les autres échantillons, environ la moitié des gènes entrent dans ces classes. Il se peut que l’on ne trouve pas de lectures de gènes à très faible expression dans les échantillons s2_r1 à r4 et, par conséquent, l’analyse comparative de l’expression avec ces échantillons sera très peu fiable et devrait être évitée.

La section 4 du protocole peut être utilisée pour identifier les valeurs aberrantes répliquées. On s’attend à ce que les répétitions d’un même échantillon, d’une même affection ou d’un même tissu présentent une corrélation plus élevée entre elles qu’avec les répétitions d’autres échantillons, affections ou tissus. La figure 5 montre une carte thermique de corrélation des deux échantillons (S1 et S2) avec cinq répétitions chacun. Les dendrogrammes en haut et sur le côté du graphique montrent deux grappes avec cinq répétitions dans chaque grappe. Le cluster de gauche contient quatre répliques de l’échantillon 1 et une réplique de l’échantillon 2 (s2_r5), le cluster de droite contient quatre répliques de l’échantillon 2 et une réplique de l’échantillon 1 (s1_r5). Dans ce cas, lorsque les noms d’échantillons s1_r5 et que s2_r5 sont à nouveau échangés, chaque cluster contient toutes les répétitions d’un échantillon, ce qui peut indiquer une erreur d’étiquetage de réplication. D’autres raisons pour lesquelles les répétitions ne se regroupent pas pourraient être un manque de différenciation entre les échantillons/conditions/tissus, ou le regroupement des réplicats pour des raisons de qualité de séquençage uniquement. Ce dernier peut se produire lorsque toutes les réplications avec un nombre de lectures exceptionnellement faible ou ceux avec un nombre de lectures exceptionnellement élevé après l’élagage et la cartographie forment un cluster.

Figure 5 : Carte thermique de corrélation d’échantillons. La carte thermique de corrélation d’échantillons est basée sur les valeurs TPM transformées log₂ et montre deux groupes distincts de cinq échantillons chacun. On s’attend à ce que les réplications biologiques/techniques présentent une corrélation plus élevée entre elles que celles provenant d’autres tissus/traitements. Le cluster de gauche contient quatre répliques de l’échantillon 1 et une réplique de l’échantillon 2 (s2_r5), le cluster de droite contient quatre répliques de l’échantillon 2 et une réplique de l’échantillon 1 (s1_r5). Les répétitions uniques doivent être vérifiées pour détecter un éventuel échange d’échantillons ou des effets de lot afin d’expliquer leur regroupement dans la carte thermique. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Tableau 1 : Comparaison des pipelines d’évaluation de la qualité du séquençage de l’ARN. Pour une analyse détaillée, voir le dossier supplémentaire 1. Veuillez cliquer ici pour télécharger ce tableau.

Fichier supplémentaire 1 : Sélection des paramètres et des références. Les analyses montrent l’influence de la sélection des paramètres et des références sur les résultats globaux du CQ. Veuillez cliquer ici pour télécharger ce fichier.

Discussion

La qualité des analyses d’expression génique différentielle dépend fortement de deux facteurs : le nombre de lectures séquencées dans chaque échantillon et répétition¹⁶ et le nombre de réplicats par échantillon^17,18. Ici, nous présentons le pipeline Rup convivial pour déterminer le nombre de lectures appropriées pour la quantification de l’expression génique dans chaque répétition. Différentes mesures permettent aux chercheurs de comprendre pourquoi les réplicats présentent de faibles nombres de lectures attribuées et d’identifier les problèmes de corrélation des échantillons. Bien que Rup ait été mis au point pour l’évaluation de la qualité des échantillons de séquençage de l’ARN des plantes, il convient également à d’autres organismes eucaryotes, ne nécessitant pas d’autres ajustements pour les échantillons non végétaux (Fichier supplémentaire 1).

La première étape de Rup détermine le nombre total de lectures avant et après le découpage des lectures. Le nombre de lectures séquencées détermine le nombre de gènes détectables, et si le nombre de lectures est trop faible, de nombreux gènes exprimés de manière différentielle ne seront pas détectés. Une grande partie des lectures éliminées au cours du processus de rognage et de filtrage de la qualité pourrait indiquer une dégradation de l’ARN, ce qui pourrait conduire à une sous-estimation de l’expression des gènes fortement dégradés. Cependant, la possibilité d’utiliser un échantillon pour des applications en aval dépend de l’organisme cible et de l’objectif de l’étude. Par exemple, pour capturer l’expression de la plupart des gènes végétaux, d’après notre expérience, 30 à 50 millions de lectures sont nécessaires, alors que pour les champignons, seulement 10 millions pourraient suffire. Ainsi, Rup ne définira pas de seuils de qualité pour l’exclusion d’échantillons problématiques, mais fournira plutôt des mesures permettant d’identifier divers problèmes pouvant être rencontrés.

La deuxième étape (qualité cartographique) évalue la précision et la fiabilité de l’alignement des lectures sur le génome de référence, en spécifiant leur aptitude aux analyses en aval. Toutes les lectures séquencées ne peuvent pas être utilisées pour le calcul de l’abondance des gènes ; Les lectures qui ne s’alignent pas sur le génome ou le transcriptome ou les lectures qui correspondent à plusieurs emplacements sur le génome sont ignorées dans la plupart des cas¹⁹ et ne contribuent pas au nombre total de lectures. Les résultats de la deuxième étape du pipeline peuvent être utilisés pour comprendre pourquoi les lectures ne sont pas utilisées dans le calcul de l’abondance. Un nombre élevé de lectures non cartographiées pourrait indiquer une contamination lors de l’extraction de l’ARN (p. ex. par des agents pathogènes des plantes ou des herbivores) ou un génome de référence incomplet. Des modèles de gènes incomplets du génome de référence peuvent entraîner un nombre élevé de lectures « sans caractéristique ». Une grande partie des lectures multi-mappées pourrait être causée par un nombre élevé de gènes d’ARNr dans la banque, indiquant une élimination insuffisante de l’ARNr pendant le processus de préparation de la bibliothèque de séquençage (dans le cas où les lectures multi-mappées sont vraiment dérivées d’ARNr, elles peuvent être analysées plus en détail en fournissant un fichier d’annotation d’ARNr au pipeline, qui calcule ensuite automatiquement le nombre de lectures possibles d’ARNr).

La troisième mesure de qualité, tout aussi importante, est la corrélation entre les répétitions d’un même échantillon. En général, la corrélation entre les répétitions d’un même échantillon doit être plus élevée que la corrélation entre les répétitions d’échantillons différents. Une faible corrélation entre les répétitions pourrait indiquer une grande variabilité biologique entre les répétitions, une grande similitude entre les échantillons/conditions/tissus, d’autres effets de lot, ou même un échange d’échantillons ou un étiquetage erroné. Rup calcule les corrélations par paires entre tous les échantillons et produit une carte thermique en grappes des corrélations d’échantillons. De plus, une analyse en composantes principales (ACP) peut être calculée sur les échantillons afin d’identifier les effets de lot possibles qui doivent être pris en compte dans d’autres analyses en aval. Bien que les effets de lot puissent être corrigés dans les analyses en aval, il peut être préférable de retirer les échantillons problématiques si la différence mesurée est trop élevée.

Le matériel d’entrée a un impact direct sur la qualité du séquençage : l’échantillonnage des tissus, l’extraction de l’ARN et la préparation de la banque sont des étapes essentielles pour minimiser les problèmes de qualité du séquençage de l’ARN. Des conditions constantes doivent être maintenues, si possible, par exemple en utilisant des chambres de culture. Les mesures de lutte antiparasitaire doivent être prises en temps opportun, car seuls les individus en bonne santé doivent être sélectionnés pour l’échantillonnage. Il est en outre conseillé de prélever des échantillons le même jour et/ou à la même heure pour réduire la variation du transcriptome circadien. De nombreux kits d’extraction d’ARN sont disponibles, et la sélection d’un kit approprié pour l’espèce cible peut améliorer la qualité de l’ARNm. Les protocoles de préparation de la banque doivent inclure des étapes d’enrichissement de l’ARN polyA+ afin de minimiser la fraction d’ARNr.

Rup peut être utilisé comme étape initiale de contrôle de la qualité dans l’analyse différentielle de l’expression génique. Ce contrôle qualité est nécessaire pour plusieurs raisons critiques, car il garantit la qualité des données d’entrée (identifie les répétitions/échantillons de mauvaise qualité, peut définir des seuils de qualité pour la profondeur de lecture et les taux de mappage) et identifie les problèmes techniques tels que les erreurs de séquençage, les effets de lot ou les erreurs d’étiquetage. Si l’entrée d’ARN est de qualité suffisante, Rup peut aider en coupant les lectures de mauvaise qualité ou en identifiant les échantillons mal étiquetés. Cependant, Rup ne peut pas compenser la mauvaise qualité de l’ARN d’entrée. Dans le cas d’ARN de mauvaise qualité ou de données de séquençage erronées, il peut être nécessaire de collecter à nouveau des échantillons, d’ajuster le protocole d’extraction de l’ARN et/ou de répéter le séquençage. Il en va de même pour les échantillons avec une grande fraction de lectures non cartographiées qui peuvent avoir été causées par une contamination. Cependant, d’après notre expérience, toutes les extractions d’ARN ne peuvent pas être simplement répétées en raison d’un manque de disponibilité du matériel. Dans ce cas, ils ne sont pas adaptés à certaines applications en aval : les échantillons avec un faible nombre de lectures mono-cartographiées ne doivent pas être analysés dans des analyses différentielles d’expression génique, mais ils contiennent toujours des informations pour l’analyse de la présence de transcrits. Dans ce cas, l’absence de transcrits dans les tissus/traitements/affections ne peut pas être prise en compte pour l’analyse et la quantification de l’abondance des transcrits.

Rup comprend certaines limitations. Par exemple, il ne teste pas directement la dégradation de l’ARN, car des mesures directes de l’intégrité de l’ARN sont nécessaires avant la préparation et le séquençage de la banque. Cependant, un script permettant d’identifier la dégradation de l’ARN à l’aide des modules RSeQC geneBodyCoverage.py et tin.py est inclus dans le référentiel de ce pipeline. De plus, Rup ne teste pas le biais de contenu GC et de longueur de transcrit. La taille de référence du génome est actuellement limitée à 4 Gb, et selon notre expérience, le module de cartographie de lecture surestime les lectures multi-cartographiées chez les polyploïdes, comme l’illustre le fichier supplémentaire dans une comparaison de la cartographie de lecture à une version haploïde par rapport aux versions diploïdes du génome d’E. californica, un génome haploïde est donc l’entrée préférée pour ce pipeline. La qualité de la sortie du Rup dépend en grande partie de la qualité du génome de référence et de l’annotation, de sorte qu’une séquence génomique incomplète ou très fragmentée pourrait conduire à une surestimation des lectures non cartographiées. De plus, l’annotation incomplète d’un modèle de gène conduit à une surestimation des lectures non attribuées. L’exhaustivité et la duplication de la séquence génomique et de l’annotation du gène peuvent être déduites à l’aide d’outils tels que BUSCO²⁰.

Rup est un outil autonome pour toutes les étapes initiales de contrôle de la qualité essentielles dans les expériences de séquençage de l’ARN. Contrairement au RSeQC et au RNA-SeQC, le prétraitement des lectures de séquençage brutes n’est pas nécessaire. L’analyse de la qualité du séquençage ne peut pas être effectuée avec l’ARN-SeQC, et les cartes thermiques de corrélation d’échantillons ne sont pas calculées par le RSeQC et l’ARN-SeQC (Tableau 1). De plus, la sortie de normalisation est en FPKM dans RSeQC et la visualisation de sortie n’est pas implémentée en tant que norme pour tous les modules en RSeQC et RNA-SeQC. Ainsi, plusieurs problèmes de contrôle de la qualité ne sont pas testés dans les deux outils alternatifs, notamment le faible nombre de lectures, la fraction élevée de lectures tronquées et l’identification des valeurs aberrantes répliquées. Une comparaison de Rup, RSeQC et RNA-SeQC est disponible dans le fichier supplémentaire 1. De plus, Rup peut être utilisé en complément de RNA-SeQC ou RSeQC, par exemple, les fichiers BAM triés produits par ce pipeline peuvent être utilisés comme entrée pour ces outils.

Actuellement, Rup est optimisé pour un petit nombre d’échantillons, mais les étapes les plus chronophages, comme le découpage de qualité et le mappage de lecture, peuvent être précalculées, par exemple, sur un cluster d’ordinateurs ou une infrastructure cloud. Pour les génomes de plus de 4 Gb, cela est obligatoire, car l’aligneur Rsubread est limité aux génomes de moins de 4 Gb. L’annotation de l’ARNr se fait avec barrnap, qui identifie les gènes de l’ARNr hautement conservés. D’après notre expérience, l’annotation du gène de l’ARNr n’a pas besoin d’être exhaustive, car même si certains gènes anormaux de l’ARNr sont manquants, un aperçu approximatif de la présence de l’ARNr dans l’ensemble de données est suffisant pour l’évaluation de la qualité de l’ensemble de données. Les futures versions de Rup pourraient passer à une méthode de mappage différente, comme l’outil de pseudo-alignement salmon²¹, pour réduire le temps d’exécution. De plus, d’autres méthodes de normalisation et de correction, comme le biais GC ou la correction du biais de longueur, pourraient être ajoutées à Rup.

En résumé, Rup fournit des informations essentielles pour assurer la fiabilité et la reproductibilité de l’analyse des données de séquençage de l’ARN. Ce pipeline rapporte de manière exhaustive les principales mesures de qualité du RNA-seq et produit des fichiers de sortie pour une utilisation directe dans les analyses en aval. Il a été conçu comme un outil autonome permettant aux chercheurs ayant des connaissances minimales en bioinformatique d’évaluer la qualité de leurs données de séquençage primaire avec une visualisation intuitive.

Matériaux

Liste des matériaux utilisés dans cet article

Nom	Entreprise	Numéro de catalogue	Commentaires
fastqcr	R	0.1.3
FUJITSU ESPRIMO D958 Desktop	FUJITSU		le pipeline a été développé et testé sur Ubuntu Linux 24.02 (32 Go de RAM)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro & nbsp ;	Pomme		le pipeline a été testé sur maxOS 15.4.1 (16 Go de RAM)
Pheatmap	R	1.0.13
R		4.4.3
Reshape2	R	1.4.4
RFASP	Bioconducteur	1.16.0
rsamtools	Bioconducteur	2.22.0
rsubread	Bioconducteur	2.20.0