Fast and accurate differential transcript usage by testing equivalence class counts

Introduction

RNA sequencing with short-read sequencing technologies (RNA-seq) has been used for over a decade for exploring the transcriptome. While differential gene expression is one of the most widely used applications of this data, significantly higher resolution can be achieved by using the data to explore the multiple transcripts expressed from each gene locus. In particular, it has been shown that each gene can have multiple isoforms, sometimes with distinct functions, and the dominant transcript can be different across samples¹. Therefore, one important analysis task is to look for differential transcript usage (DTU) between samples.

DTU can be inferred through differential exon usage, where the proportions of RNA-Seq fragments aligning to each exon change relative to each other between biological groups. Anders et al.² showed that exon counts could be used to test for differential exon usage with a generalized linear model that accounts for biological variability. However, counting fragments across exons is not ideal because many fragments will align across multiple exons, making their assignment to an individual exon ambiguous. Moreover, individual exons often need to be partitioned into multiple disjoint counting bins when exon lengths differ between transcripts. Typically, there will be more counting bins than transcripts, resulting in lower power to detect differences between samples.

An alternative to using exon counts for testing DTU is to perform tests directly on estimated transcript abundances³. Recently, fast and accurate methods for quantifying gene expression at the transcript level have been developed^4,5. These methods use transcript annotations that include multiple known transcript sequences for each gene as a reference for the alignment. The expression levels of individual transcripts can be estimated from pseudo-aligned reads that are compatible with transcripts associated with a specific gene⁶. Transcript abundance estimates can be used as an alternative starting measure for DTU testing, which has been shown to perform comparably with state-of-the-art methods³. In addition, pseudo-alignment is significantly faster than methods that map to a genome. However, in the most comprehensive comparison using simulated data, exon-count based methods were shown to have slightly better performance compared with methods that first estimate transcript abundances³.

Conceptually, quantification by lightweight or ‘pseudo’ alignment begins by using a transcript annotation as a reference and then assigns each read as ‘compatible’ with one or more transcripts that are a close alignment to the read⁴. Because different transcripts of the same gene share large amounts of sequence, many reads are compatible with several transcripts. Reads are therefore assigned to an equivalence class, or transcript compatibility class, which reflects the combination of transcripts compatible with the read sequence (Figure 1). For the purposes of this work, we consider an equivalence class to be defined as in Bray et al.⁴, i.e. any fragments that are pseudo-aligned to the same set of transcripts are considered to be part of the same equivalence class. Figure 1 shows a toy example of a gene with three different transcripts. Depending on its sequence, a read can align to all three transcripts, only two of the transcripts or just one transcript. These different combinations result in four possible equivalence classes, containing read counts, for this gene.

Figure 1. The use of equivalence classes for detecting differential transcript usage (DTU) in a hypothetical gene.

The example shows a gene consisting of six exons (Ex1-6) and three transcripts (t_1-3) resulting in four equivalence classes (EC1-4). t₁ is predominantly expressed in condition 1 (S1), whereas t₃ is predominantly expressed in condition 2 (S2). The DTU is evident as a change in the relative counts for EC2, EC3 and EC4 between conditions. The pipelines for the three alternative methods for detecting DTU are shown: quantification of transcript expression followed by DTU testing, assignment of read counts to equivalence classes followed by testing of equivalence class counts (DECU) and assignment of read counts to exons followed by differential exon counts (DEU). Genes that are detected to have DECU or DEU are inferred to have DTU. The transcript quantification table in the left-most column is example data only, and is not based on real inference.

Recently, equivalence classes have been used for clustering single-cells^7,8 and Yi and colleagues have recently introduced direct differential testing on equivalence classes in a catch-all method to identify genes that display any transcript-level phenomena such as cancellation (isoform switching), domination (high abundance isoform(s) that mask transcript-level differences) and collapsing (multiple transcripts exhibiting small changes in the same direction)⁹. Here we focus on the case of isoform switching using methods originally designed for testing exon read counts. We evaluate the appropriateness of equivalence class read counts as an alternative choice for quantification compared to exon- and transcript-level quantification. We propose that DTU can be more accurately detected using equivalence class counts directly, rather than using these counts to first estimate individual transcript abundances before performing DTU. Soneson et al. applied a conceptually similar method with MISO¹⁰ by defining counting bins as combinations isoforms and counting according to isoform compatibility³. In our scenario, count-based DTU testing procedures such as DEXSeq are applied directly to equivalence classes generated from fast lightweight aligners, such as Salmon and Kallisto. DTU testing on equivalence class counts is not only fast but also bypasses inherent uncertainty in directly estimating transcript abundances before statistical testing.

We evaluate the performance of DTU testing on equivalence class read counts using real and simulated data, and show that the approach yields higher sensitivity and lower false discovery rates than estimating counts from transcript abundances, and performs faster with accuracies similar or better than counting across exons.

Results

Here we propose an alternative pipeline for performing DTU and evaluate its performance using simulated and real datasets²³. The method we propose is to first perform alignment with a lightweight aligner and extract equivalence class (EC or transcript compatibility) counts. These EC counts are assigned to genes using the annotation of the transcripts matching to the EC. Next, each gene is tested for DTU between conditions using a count based statistical testing method where exon counts are replaced with EC counts (Figure 1). Significant genes can then be interpreted to have a difference between the relative abundance of transcripts of that gene between conditional groups. In evaluating the EC approach, we used Salmon for pseudo-alignment and DEXSeq for differential testing. We then compared DTU results against the alternative quantification and counting approaches, also using DEXSeq for testing (see Methods). It should be noted that we are not attempting to evaluate the statistical testing method (DEXSeq) in relation to other methods, as this has been done previously in several papers^3,9.

The datasets we used to evaluate performance were simulated data from human and drosophila from Soneson et al.³ and biological data from Bottomly et al.¹¹. Each of the Soneson datasets consisted of two sample groups, each with three replicates, where 1000 genes were randomly selected to have DTU such that the expression levels of the two most abundant transcripts were switched. The Bottomly dataset contains 10 and 11 replicates each from two mouse strains that were used to call truth and then were subsampled to three replicates in the testing scenarios.

Fewer equivalence classes are expressed than exons

The number of counting bins used for DTU detection has an impact on sensitivity. More bins leads to lower average counts per bin and therefore lower statistical power per bin and more multiple testing correction. We therefore examined the number of ECs, transcripts and exons present in each dataset. Although the theoretical number of ECs from a set of transcripts can be calculated from the annotation and has the potential to be large, not all combinations of transcripts exist or are expressed. The number of equivalence classes calculated from pseudo-alignment depends on the experimental data as only ECs with reads assigned to them are reported. We compared the number of transcripts and exons in the three datasets (with at least one read) to the number of ECs. In both the simulated human and drosophila datasets, as well as in the Bottomly mouse data, the number of ECs is greater than the number of transcripts, but substantially fewer than the number of exons, indicating that there might be more power for testing DTU using ECs, compared to exon counts (Figure 2a).

Figure 2. The number of counting bins and variance between replicates.

(a) The number of transcripts, equivalence classes and exons per gene, where each feature has at least one associated read. (b) The density of the log₂ of the variance of counts over the mean for each feature (calculated per condition).

Performance of equivalence classes for DTU detection

Several methods were previously tested on the simulated data from Soneson et al.³; DEXSeq’s default counting pipeline and featureCounts were shown to perform best. We recalculated exon counts using DEXSeq’s counting pipeline (as recommended by Soneson et al., we excluded region of genes that overlapped on the same strand in the input annotation) and ran Salmon⁵ to obtain both transcript abundance estimates and equivalence class counts. All other comparison results were obtained from Soneson et al.³. For the simulated datasets, we found that ECs had the highest sensitivity in both the drosophila and human datasets (Figure 3a) with a TPR of 0.697 and 0.739 respectively (FDR < 0.05). However, ECs also had a slightly higher FDR than exon-counting methods.

Figure 3. The performance of the equivalence class method for differential transcript usage.

(a) The equivalence class method compared to other state-of-the-art methods on simulated data described in Soneson et al.³. (b) The ability of the equivalence class, transcript and exon-based methods to recreate the results of a full comparison (10 vs. 11) of the Bottomly data, using only a (randomly selected) subset of samples (3 vs. 3) across 20 iterations. The union of all genes called as significant across all three methods is used to calculate the FDR, and the intersect (genes called by all three methods) is used for the TPR. Full results (union, intersect and each method’s individual truth set) is shown in Supplementary Figure 3.

We next tested the performance of the EC method on a biological dataset from Bottomly et al. We tested the complete RNA-seq dataset (10 vs. 11) for DTU using DEXseq on counts generated from transcript abundance estimates, exons and ECs. To calculate the FDR, we considered the set of 'true' DTU genes to be the union of all genes called significant (FDR < 0.05) across the three methods. To calculate the TPR, the intersect of genes called by all methods was used. Supplementary Figure 2²¹ shows the number of significant genes and overlap between all three methods. ECs called the highest number of genes with significant DTU (1485 genes, in contrast to the 748 and 391 genes called significant by the transcript and exon-based methods respectively). Similar to the FDR experiments described in Pimentel et al.¹², we randomly selected three samples per condition and performed DTU using all three methods and repeated this for 20 iterations. Figure 3b shows the results. EC-based testing performed the best, with a mean FDR of 0.305 across all iterations (compared to a mean FDR of 0.569 and 0.373 for the transcript and exon-based methods respectively). The mean TPR was also slightly higher for ECs at 0.544, compared to exons at 0.539 and 0.36 for the transcript-based method. Results for all three combinations of the ‘truth gene’ sets (union, intersect and individual) are shown in Supplementary Figure 3²¹. The EC-based method had consistently lower FDR, which is also illustrated by the rank-order plot (Supplementary Figure 4²¹), showing the number of false positives present in the top 500 FDR-ranked genes. In terms of the TPR, ECs performed better than transcripts, but worse than exons when using the union of all methods as the truth set. In the Bottomly analysis, Salmon was used as a representative method for transcript abundance estimation. We also performed the analysis with Kallisto, which gave results consistent with Salmon (Supplementary Figure 5²¹).

Discussion

DTU detection has previously been approached by either testing for changes to the read counts across exons or changes in the relative abundance of transcripts. These approaches are intuitive but are not necessarily optimal for short read data analysis. In particular, individual exons are not necessarily the optimal unit of isoform quantification as there are often many more exons than transcripts. In addition, transcript quantification can be difficult because read assignment is ambiguous. Fortunately, transcript quantification methods generate equivalence class counts as a forestep to estimating abundances. We propose that equivalence classes are the optimal unit for performing count based differential testing. Equivalence class counts benefit from the advantages of both exon and transcript counts: they can be generated quickly through pseudo-alignment, there are fewer expressed than exons, and they retain the low variance between replicates seen in exon counts compared to transcripts abundances.

Here we evaluated the use of equivalence classes as the counting unit for differential transcript usage. We used two simulated datasets from drosophila and human and one biological dataset from mouse. Our results suggest that equivalence class counts provide equal or better accuracy in DTU detection compared to exon counts or estimated transcript abundances. We also found the analysis was quick to run and we provide code to convert pseudo alignments into gene level EC annotations.

The ECs used in our evaluation are defined using only the set of transcripts for which reads are compatible. Extensions to this model have been proposed that incorporate read-level information, such as fragment length, to more accurately calculate the probability of a read arising from a given transcript¹³. Although, we do not consider probability-based equivalence classes in this work, incorporating this information for DTU deserves exploration in future work. In addition, EC counts may be calculated from full read alignment rather than pseudo-alignment^14,15, which has the potential to improve accuracy further. In this work, we limited our investigation to comparing the best counting metric preceding DTU statistical testing, using DEXSeq as a representative method. Evaluation of statistical testing methods for DTU is outside the scope of this manuscript and would require further work.

One limitation of using equivalence classes is in the interpretation of the results. Although we can detect DTU at the gene-level, it is not simple to determine which isoforms have changed abundance without further work. We propose that superTranscripts¹⁶, which are a method for visualising the transcriptome, could be used for interpretation. Alternatively, transcript abundances, which are generated together with ECs, can still be used to provide insight into the isoform switching.

Finally, in this work, we have focused on differential transcript usage, but EC counts have the potential to be useful in a range of other expression analysis. EC counts have already been applied to areas such as clustering and dimensionality reduction⁷, gene-level differential expression⁹, single-cell transcriptomics^7,8 and fusion detection²⁰. We foresee that equivalence classes could serve as a base unit of measurement in many other types of analyses.

Methods

We detected sequence content bias in the Bottomly RNA-seq data using FastQC v0.11.4, and therefore performed trimming using Trimmomatic¹⁴ 0.35, using recommended parameters. The simulated Soneson data was not trimmed.

To obtain transcript abundance counts, Salmon⁵ v0.13.0 (development version) was run on the drosophila, human and Bottomly datasets in quant mode using default parameters. To obtain EC counts, the --dumpEq argument was used, as well as the --skipQuant to skip the quantification step. Kallisto⁴ 0.43.0 was run in pseudo mode with the --batch argument to run all samples simultaneously. Fragment length and standard deviation were estimated from all reads of a single sample from the Bottomly data (SRR099223). Equivalence classes were then matched between samples and compiled into a matrix using the python scripts (create_salmon_ec_count_matrix.py and create_kallisto_ec_count_matrix.py), available on GitHub and archived on Zenodo²². Equivalence classes mapping to more than a single gene were removed. No other filtering was performed on any of the data types.

To perform the exon-based counts, raw reads were first aligned using STAR¹⁵ v2.5.2a, then the DEXSeq-count annotation was prepared excluding overlapping exon-parts, from different genes, on the same strand (--aggregate=’no’). DEXSeq-count was then run using default parameters. The same genome and transcriptome references for drosophila and human were used as in Soneson et al.³, with only protein-coding transcripts considered for the Salmon index. For the Bottomly data, we used the NCBIM37 mm9 mouse genome and Ensembl release 67 transcriptome. Non-protein-coding transcripts were filtered out, as with the Soneson transcriptome reference. DEXSeq v1.26 was used to run all DTU analyses.

An earlier version of this article can be found on bioRxiv (DOI: https://doi.org/10.1101/501106).

Data availability

Software availability

Pipeline used to reproduce the quantification data generated in this paper:https://github.com/Oshlack/ec-dtu-pipe.

Archived source code at time of publication:https://doi.org/10.5281/zenodo.2567596²³.

Source code to run the analyses and generate the paper figures:https://github.com/Oshlack/ec-dtu-paper.

Archived source code at time of publication:https://doi.org/10.5281/zenodo.2561549²².

License: MIT license.