Quality assessment of long read data in multisample lrRNA-seq experiments using SQANTI-reads

SQANTI-reads can be used to evaluate long-read technology improvements

Long-read methods are rapidly improving, and both Nanopore and PacBio are updating their instruments, molecular protocols, and algorithms. Tools that can readily evaluate the broad impact on data quality of these technical improvements are highly needed to make informed decisions for downstream analyses. Our Drosophila data set contained the same set of samples processed on the MinION and PromethION instruments through several runs. In addition, the development of basecalling algorithms for ONT sequencing is a highly active and rapidly evolving field (Pagès-Gallego and de Ridder 2023; Diensthuber et al. 2024). We expect that as new basecallers emerge and algorithms are improved upon, users will want to evaluate the impact of different basecallers and algorithms on the quality of their experiments. This comparison can be done using SQANTI-reads. While Dorado is now the default basecaller for ONT, we demonstrate the utility of SQANTI-reads by comparing real-time Guppy basecalled reads and Dorado basecalled reads. We first evaluated whether Dorado effectively improved data quality without introducing biases and if data from several sequencing experiments could be merged. We compared the Guppy and Dorado basecallers using SQANTI-reads metrics. As anticipated, Dorado resulted in more reads with assignable barcodes, a higher number of mapped reads, more reads aligning to annotated genes, more reads aligning to annotated transcripts, and longer reads, without an increase in the proportion of reads with technical artifacts (Supplemental Fig. 1; Supplemental File 1). This confirms that Dorado improves basecalling accuracy without introducing unwanted biases and motivated the selection of Dorado basecalled reads for subsequent analyses.

In the Drosophila experiment, libraries were barcoded, pooled, and multiplexed across different MinION and PromethION runs, with a re-pooling step between the two machines. We used SQANTI-reads to compare the quality of the MinION and the PromethION runs, and to evaluate the consistency of the MinION and PromethION technology across TRs. The first MinION run (TR1) had higher percentages of reads classified as novel in catalog (NIC) and novel not in catalog (NNC) and with noncanonical junctions compared to the second and third MinION runs (TR2, TR3) and compared to the PromethION runs for the same libraries (Supplemental Fig. 2). NIC refers to reads with novel combinations of annotated junctions and NNC refers to reads with at least one unannotated junction. MinION runs TR2 and TR3 were similar in their quality metrics (described below) to the PromethION run of the same samples, and the TRs on the PromethION were similar. These results indicate that the technology performs consistently across instruments and runs. Based on these SQANTI-reads QC results, we aggregated data across TRs to further evaluate the quality of the lrRNA-seq experiment. Although TR1 had lower quality than the other TRs, the overall read numbers were low, and we decided to keep this TR in our evaluation of the samples.

SQANTI-reads metrics can be used to evaluate the global quality of the lrRNA-seq experiment

In a multisample lrRNA-seq experiment, all samples should be of similar quality. SQANTI3 uses the full-splice match (FSM) structural category to identify long-read sequences whose junctions are consistent with an annotated transcript model. However, for an lrRNA-seq experiment to accurately reflect the analyzed transcriptome, the reads should also capture the distribution of transcript lengths of the expressed transcriptome. The distribution of transcript lengths depends on many factors, including the species and tissue used as input. Consistent with the species transcript model annotations, we observe that Drosophila has less complex and shorter transcripts than humans (Supplemental Fig. 3). A data set with reads substantially shorter than the annotated, transcriptome but a high proportion of FSM indicates capture of short transcripts, while combining a high proportion of incomplete-splice match (ISM) may indicate RNA degradation. ISM refers to a read whose junction matches an annotated transcript but is missing junctions on the 5′ end, 3′ end, or both ends. We looked at these values for an initial assessment of the quality of the Drosophila experiment.

First, we compared the number of reads and distributions of read lengths for all samples. The difference in sequencing depth (number of reads) between the two developmental stages was evident, as expected because of the additional PromethION run on the 3–8 day samples (Fig. 1A). For all samples, most reads were shorter than 1 kb, with <20% of them above the 1 kb threshold (Fig. 1B; Supplemental Figs. 4A, 5). When we evaluated all reads, 53% to 67% of the reads across samples were classified as FSM, 20% to 38% were labeled as ISM. For all samples, the proportion of NIC/NNC was <10% of the reads (Fig. 1C). For the reads that were above the 1 kb threshold, between 73% and 82% of the reads were FSM, with only 10% to 18% ISM (Fig. 1D). The decrease in proportion of ISM and increase in proportion of FSM when evaluating only the reads >1 kb suggests that the shorter reads represent incomplete transcript sequences, potentially due to mRNA degradation.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

SQANTI-reads analysis of Drosophila samples. (A) Number of mapped reads by experimental group labeled with read length. (B) Percentage of mapped reads >1 kb versus percentage of reads that are FSM. Dots represent the early stage (0–1 h after eclosion) and crosses indicate the adult stage (3–8 days old). The four genotypes are indicated with four different colors. (C) Percentage of reads mapping to genes in each SQANTI3 structural category (D) Percentage of reads mapping to genes in each SQANTI3-QC structural category for reads >1 kb (E) Number of genes detected with breakdown by the number of reads mapped to each gene. (F) Number of UJCs detected with breakdown by the number of reads associated with each UJC. (G) Percentage of UJCs detected with breakdown by the number of reads associated with each UJC. (H) Distribution of the percentage of FSM reads by gene across samples.

Reads that share all internal junctions are grouped together and annotated using a string made up of the junction locations, the unique junction chain (UJC) (Nanni et al. 2024). We examined both gene-level and UJC-level metrics. We found that, despite the large sequencing depth differences between developmental stages, the number of detected genes was only slightly lower in the 0–1 h samples (Fig. 1E). However, these genes were quantified with fewer reads (80% genes with <50 reads) than the 3–8 days samples, which had between 30% and 50% of genes with more than 100 reads (Fig. 1E). When evaluating UJC, we found that, while the number of UJC mirrored the sequencing depth pattern (Fig. 1F), with 3–8 days samples showing five times more UJC than 0–1 h samples, and a larger number of FSM and ISM UJC, there were many additional UJC detected by fewer than 10 reads, and usually by a single read (Fig. 1F,G). These UJC were most frequently NIC/NNC (Supplemental Fig. 4B,C). Downstream analyses would therefore need to address whether this represents novel low-expressed transcripts or technology errors. In contrast, the percentage of FSM reads between the two time points differed by less than 1× in all replicates (Fig. 1H). These results indicate that the higher sequencing depth of the 3–8 days samples does not change the number of detected genes or annotated transcripts (FSM). The higher read depth per gene/UJC suggests that more genes and transcripts will be able to be quantitatively evaluated in the 3–8 day samples compared to the 0–1 h samples.

In the Drosophila data, we noticed two samples (Sample 4 Rep 1 [RIL 12279, 0–1 h]—red arrow; Sample 6 Rep 3 [RIL 11255, 3–8 day]—teal arrow) that had the lowest percentage of FSM and highest percentage of ISM in the 0–1 h and 3–8 day groups, respectively (Fig. 1C). To determine whether these two samples were of overall lower quality than the rest, we examined their SQANTI-reads metrics. We found that Sample 4 Rep 1 had a lower proportion of FSM across all genes (Fig. 1H) and a higher proportion of genes quantified with only one read (Fig. 1F), while Sample 6 Rep 3 had a similar gene (Fig. 1H), UJC (Fig. 1G), and % FSM in genes (Fig. 1H) than other 3–8 day samples. We concluded that RIL 12279 rep 1 0–1 h is a lower-quality sample.

Altogether, this example shows that SQANTI-reads metrics can be used to compare samples and experimental conditions in a multisample experiment, detect outliers, and suggest points of attention for downstream data processing.

SQANTI-reads metrics can be used to identify systematic differences among samples

The previous example demonstrated that SQANTI-reads metrics are effective in assessing data set consistency. However, SQANTI-reads evaluates over 35 quality metrics, making it challenging to determine which features contribute to potential differences among samples. We include principal component analysis (PCA) analysis to identify which metrics are the most relevant for quality variability when there are differences among samples or between groups. The percentage of reads and UJCs in each structural category, percentage of artifact reads (RT-switching, noncanonical junctions, and intrapriming), percentage of junctions in each category, as well as length metrics, are included in the PCA.

We applied SQANTI-reads PCA analysis of quality features to investigate differences in read quality among various LRS methods used in the Long-read RNA-seq Genome Annotation Assessment Project (LRGASP) challenge (Pardo-Palacios et al. 2024b), focusing on the WTC11 data set. We evaluated triplicate transcriptome measurements of the WTC11 human cell line, analyzed using three technologies: cDNA PacBio Sequel II (cDNA PacBio), cDNA Oxford Nanopore MinION (cDNA ONT), and direct RNA Oxford Nanopore MinION (dRNA ONT). The analysis revealed that WTC11 samples clustered based on the long-read technology applied (Fig. 2A). Specifically, PC1, which explains 56% of the variance, distinguished cDNA ONT samples from those generated by the other two technologies, while PC2, accounting for 35% of the variance, highlighted differences between dRNA ONT and cDNA PacBio. To further explore these differences, we examined the loadings for each principal component. Quality features with the highest positive loadings in PC1 included the number of reads, the percentage of reads, and the proportion of UJCs in the NNC category, while features with high negative loadings included Intergenic and Genic Genomic reads. Several junction-related variables also exhibited high absolute loadings on PC1 (Fig. 2B). SQANTI-reads plots confirmed these structural category differences between cDNA ONT samples and other library preparations. cDNA ONT had both the highest proportion of NNC reads and UJCs (Fig. 2C,D) and also had the lowest proportion of intergenic reads (Fig. 2C). Other differences in sequencing throughput and junction characteristics were also confirmed (Supplemental Fig. 6).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

SQANTI-reads PCA analysis of LRGASP WTC11 samples. (A) PCA using SQANTI-reads quality features. The percentage of variance explained by each principal component (PC) is labeled on each axis in brackets. (B) Top 10 Loadings for PC1 and PC2. (C) Distribution of reads in SQANTI3 structural categories. (D) Distribution of UJCs in structural categories. (E) Distribution of read lengths for all mapped reads.

Upon examining the feature loadings for PC2, we found that variables with high contributions included several metrics related to read length (Fig. 2B). Consequently, we evaluated the SQANTI-reads “Lengths of All Mapped Reads” plot for this experiment. Indeed, we observed that the cDNA PacBio method produced a significantly higher proportion of reads between 1–2 kb, 2–3 kb, and >3 kb, as suggested by their negative loadings, compared to the dRNA ONT method, which predominantly generated reads shorter than 1 kb (Fig. 2E). Similarly, the percentage of reads assigned as ISM with high positive values was higher in dRNA ONT samples (Fig. 2C).

In conclusion, we showed that the SQANTI-reads PCA analysis is an effective tool for uncovering significant read quality differences between LRS methods. The technological differences explored here are based on a benchmarking experiment that has been superseded by new technology and the results here should not be interpreted as an evaluation of current technological capabilities. Nonetheless, we have shown that SQANTI-reads can be used to evaluate LRS methods and for revealing technological biases.

SQANTI-reads identifies potentially under-annotated genes

Long-read data often contain a large number of sequences that cannot be exactly matched to existing annotations. In many cases, these UJCs belong to annotated genes and are identified by only a few reads, as illustrated by the SQANTI-reads analysis in Figures 1 and 2. This suggests they could be either low-expressed transcripts or technological artifacts. However, in some instances, a high proportion of reads in a gene may correspond to the same novel UJC, indicating the possibility of a previously unannotated transcript that warrants closer examination. SQANTI-reads includes a customizable decision tree to identify such cases (see Methods; Fig. 3A). Basically, the tool identifies genes with a high (R > threshold) number of reads, with novel UJC containing a large fraction (Q > threshold) of the gene's splice sites and capturing a high proportion (P > threshold) of the reads. We applied this approach to the WTC11 PacBio data using default parameters (R = 100, P = 20, and Q = 80). In addition, we classify expressed genes as well annotated or underannotated using the parameters P and R. The logic for defining expressed genes as well annotated or underannotated is described in Figure 3A. The annotation category for all expressed genes (default: number of reads in gene [R] > 100) is provided in the gene_classification.csv file.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Evaluation of WTC11 PacBio samples for under-annotated genes. (A) Decision tree for classifying genes as well-annotated or under-annotated and classifying transcripts as putative novel transcripts. (B) Number of genes by their annotation status according to SQANTI-reads parameters. (C) The coverage (percent of total reads) versus length (percentage of maximum junctions) for all UJCs in under-annotated genes with well-covered candidate transcripts. UJCs that meet the thresholds for putative novel transcripts are colored in orange. (D,E) UpSet plot of putative novel transcripts with the most expressed annotated transcript in that gene for well-annotated genes with well-covered FSMs (D) and with an FSM detected but without <20% of the total reads in the gene (E). Upper and middle panels show the distribution of the number and proportion of nucleotides different between the putative novel transcript and most expressed FSM.

From the set of expressed genes, we identified 8556 well-annotated genes, 88% of which have a well-covered annotated transcript (>20% of total gene coverage) (Fig. 3B; Supplemental Fig. 7). We also identified 101 genes for which there are no reads with an FSM match to an annotated transcript. Of these, 54% have a well-covered UJC. For all expressed genes with a well-covered unannotated UJC, we identified 424 that contained most of the observed junctions in that gene (>80%), and we label these putative novel transcripts (Fig. 3C; Supplemental Fig. 7). Of these, 316 were NIC and 108 were NNC (Supplemental Fig. 7). The SQANTI-reads output for putative novel transcripts is included in the putative_underannotation.csv file (Table 1).

For genes with at least one putative novel transcript (R > 100, P > 20, and Q > 80) and an annotated transcript (FSM), we selected the FSM with the highest proportion of reads. We then compared the structure of the putative novel transcripts to the most expressed annotated transcript using TranD (Nanni et al. 2024). For genes with an annotated transcript that is relatively highly expressed (>20% of the reads for that gene), 103 putative candidate transcripts differed from the annotated transcript by donor/acceptor variation, suggesting a possible alternative splice site. In addition, 10 putative candidate transcripts had an extra exon, 15 a skipped exon, and 9 with both missing and skipped exons relative to the most expressed annotated UJC (Fig. 3D). For the genes where the annotated transcript represented <20% of the reads in that gene, the putative novel transcript differed from the annotated transcript by an alternative exon in 147 cases (33 extra exons, 86 skipped exons, and 43 with both an extra and skipped exon) (Fig. 3E). Details for this analysis are provided in the Supplemental Methods.

This analysis shows that SQANTI-reads can readily identify under-annotated genes and flag putative novel transcript models that contain interpretable alternative exonic patterns that deserve further attention.

SQANTI-reads metrics for donors/acceptors identify noisy splicing and potentially novel splice sites

SQANTI-reads calculates the mean, standard deviation, and coefficient of variation (CV) for all expressed annotated donors/acceptors using the absolute distances of each read donor and acceptor from the nearest annotated donors and acceptors (Fig. 4A). A CV > 0 indicates variability in the donor/acceptor, with higher CV values indicating more variability. Variability around a splice junction may be due to weak splicing (Wang and Marín 2006) or to technology errors and mapping accuracy, for example, due to junction ambiguity (Li 2018). We evaluated these metrics on the WTC11 data set for reference junctions with at least 10 reads. We found similar patterns in the variability (CV > 0) in donors and acceptors (Fig. 4B). All three technologies identify donors and acceptors with variability around the splice site (CV > 0). Donors/acceptors with CV > 0 consistently across the three technologies are highly suspicious of “noisy” splicing or a weak splice site and may warrant follow-up (Supplemental Fig. 8A,B). The SQANTI-reads output file cv.csv identifies the donors/acceptors with CV > 0, making it straightforward to follow up on particular locations with tools such as the Integrative Genomics Viewer (IGV) (Robinson et al. 2011).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Variation in donors/acceptors. Metrics are only calculated for annotated donor/acceptors with a minimum threshold of reads (10 by default). (A) Metrics for the classification of donor/acceptor variation. When all reads align to the annotated donor/acceptor this is classified as a reference match (Ref Match). When all reads align to the same donor/acceptor location, but this is not the annotated position this is classified as CV = 0. When reads align in multiple positions in proximity to an annotated donor/acceptor this is classified as CV > 0. (B) Classification of the number (left) and percentage (right) of detected acceptors faceted by technology. (C) Classification of the number (left) and percentage (right) of detected donors faceted by technology.

A reference match junction indicates that the splice signal is strong (Wang and Marín 2006; Dent et al. 2021). Results for cDNA PacBio and dRNA ONT were similar, with both showing a higher number and proportion of reference match donors/acceptors compared to cDNA ONT (Fig. 4B,C), despite these technologies detecting similar numbers of FSM UJCs (Fig. 2D). We compared the FSM UJCs identified by the three methods (Supplemental Fig. 9A). Most of the FSM UJCs were detected by all three technologies (n = 19,690), with a similar number detected by cDNA PacBio only (n = 8110) and cDNA ONT (n = 9360) only. We hypothesized that the difference in the number of reference match donors/acceptors was potentially due to longer transcripts with more junctions being detected in cDNA PacBio compared to shorter transcripts with fewer junctions in dRNA ONT and cDNA ONT. For the FSMs detected only in one technology, we plotted the distribution of the number of junctions and confirmed that the cDNA PacBio FSM transcripts had a larger number of junctions compared to dRNA ONT and cDNA ONT (Supplemental Fig. 10). This agrees with cDNA PacBio showing longer reads than both ONT technologies in the WTC11 data set (Fig. 2E).

Splice junctions may differ from annotated sites due to the presence of novel donors/acceptors. The category CV = 0 identifies donors/acceptors with no variability but differing from the annotated donor/acceptor, representing strong candidates for bona fide alternative splice sites. We evaluated the donors and acceptors with CV = 0 in all the technologies (Supplemental Fig. 8E,F). We identified 51 donors and 34 acceptors with CV = 0 in all three technologies. Of these, 76% of donors and 65% of acceptors are within 12 nt of an annotated donor/acceptor, indicating a potential misannotation of the splice site (Supplemental Fig. 11). These donors and acceptors with CV = 0 across all technologies were detected by a minimum of 10 reads and, in some cases, could be detected with >100 reads (Supplemental Fig. 12). Supplemental Figure 13 shows an example of one of the reference donors with CV = 0 detected across all three technologies. This donor is in the H2AZ1 gene at position 99,948,814 on Chromosome 4. All the reads associated with this donor map to position 99,948,811, which is 3 nt away from the annotated donor position (Supplemental Fig. 13). Donors and acceptors with CV = 0 detected consistently across samples and/or technologies indicate potential robust detection of alternative splice sites.

Memory usage and runtime of SQANTI-reads

Computational efficiency is a critical factor for users when choosing a quality control tool. To evaluate memory usage and running time of SQANTI-reads, data sets containing five different number of reads were generated, as well as five replicates for each data set. After independently processing these data sets through SQANTI3, they were tested five times in 15 combinations to analyze memory usage and runtime for: a single sample of each size, three samples of each size, and five samples of each size.

The performance evaluation of SQANTI-reads revealed that memory consumption remained consistent at 2 GB of RAM across all tested conditions, regardless of the number of reads or the number of samples analyzed. Runtime, however, increased with both the number of reads per sample and the total number of samples processed simultaneously (Supplemental Fig. 14).

When processing a single sample, runtimes ranged from 100 sec for data sets with 10,000 reads to ∼3 min for data sets with 1,000,000 reads. For multisample runs, processing three samples of 1,000,000 reads each required ∼4 min, while processing five samples of the same size took about 5 min. These findings highlight the scalability of SQANTI-reads, with a predictable increase in runtime proportional to the data set size and sample number. Overall, SQANTI-reads demonstrated computational efficiency with minimal memory requirements, making it suitable for large-scale transcriptomic analyses.