satuRn: Scalable analysis of differential transcript usage for bulk and single-cell RNA-sequencing applications

satuRn model

As input, satuRn requires a matrix of transcript-level expression counts (or, analogously, a matrix of exon-level or equivalence class-level counts), which may be obtained either through pseudo-alignment using kallisto¹ or Salmon,² or by classical alignment-based tools followed by transcript-level quantification (e.g. STAR^33,34 and RSEM³⁵). Let Y_gti denote the observed expression value for a given transcript t = 1, … , T_g of gene g = 1, ..., G in cell or sample i = 1, … , n. The total expression of gene g in sample i can then be expressed as

i.e., by taking the sum of expression values for all T_g transcripts belonging to gene g in sample i. The usage of transcript t in sample or cell i can then be estimated as

Next, we adopt a quasi-binomial (QB) generalized linear modelling (GLM) strategy to model transcript usage. As opposed to canonical maximum likelihood models, this quasi-likelihood modelling strategy only requires the specification of the first two moments of the response distribution, i.e., the mean and the variance. We define the mean of the QB model as

In this notation, π_gti is the expected probability of observing transcript t within the pool of transcripts (1, … , T_g) belonging to gene g in sample i and, as such, corresponds to its expected usage for that sample. We model π_gti using a logit link function, where β_t is a p x 1 column vector of regression parameters modelling the association between the average usage and the covariates for transcript t. Finally, ${\mathit{X}}_i^T$ is a row in the n x p design matrix X that corresponds with the covariate pattern of sample i, with p the number of parameters of the mean model, i.e., the length of vector β_t.

The variance of the QB model can be described as

with $Y_{g.i}{\mathrm{\pi }}_{\mathit{gti}}\left(1-{\mathrm{\pi }}_{\mathit{gti}}\right)$ the canonical variance of the binomial distribution and ϕ_gt a transcript-specific overdispersion parameter to describe additional variance in the data with respect to the binomial variance. We adopt the empirical Bayes procedure from Smyth et al.,²³ as implemented in the squeezeVar function of the limma Bioconductor R package, to stabilize the estimates of ϕ_gt by borrowing information across transcripts, which is adopted in the default edgeR quasi-likelihood workflow for bulk RNA-seq data.^22,36 Note that stabilizing the dispersion estimation is particularly useful in datasets with a small sample size.

Taken together, the quasi-binomial thus allows us to model the log-odds of drawing a particular transcript t from the pool of transcripts from the corresponding gene g across samples. The intercept also has an interpretation of a log-odds and the remaining mean model parameters are log-odds ratios, which may thus be interpreted in terms of differential transcript usage. We adopt t-tests that are computed based on the log-odds ratio estimates of the QB model and the posterior variance, as obtained from the empirical Bayes procedure. P-values are computed assuming a t-distribution under the null hypothesis with posterior degrees of freedom calculated as the sum of the residual degrees of freedom and the prior degrees of freedom from the empirical Bayes procedure.

For bulk analyses, the implementation of satuRn as described above provides a high performance and a good control of the FDR. However, for single-cell datasets we observed that our inference is too liberal (see Extended data²⁵ figure S14), which could suggest that the theoretical null, the t-distribution, is no longer valid. Indeed, in large-scale inference settings, failure of the theoretical null distribution is often observed. Efron³⁷ (Chapter 6) describes four reasons why the theoretical null distribution may fail; failed mathematical assumptions, correlation across features (transcript expression), correlation across subjects (samples or cells), and unobserved confounders in observational studies. To avoid these issues, Efron proposes to exploit the massive parallel data structure of omics datasets to empirically estimate the null distribution of the test statistics.³⁷ To this end, Efron converts the test statistic to z-scores, which should follow a standard normal distribution under the theoretical null, and then proposes to approximate the empirical null distribution with a normal distribution with unknown mean (${\mu }^{\ast }$) and standard deviation (${\sigma }^{\ast }$), which can be estimated by maximum likelihood on a subset of the test statistics near zero.

As such, we first convert the two-sided p-values to z-scores according to

with Φ the cumulative distribution function for the standard normal distribution, p_gt the original two-sided p-value indicating the statistical significance of differential usage of transcript t from gene g between the conditions of interest, sign(S_gt) the sign of the t-test statistic S_gt and Z_gt the resulting z-score. Next, we adopt the maximum likelihood procedure, implemented in the locfdr function of the locfdr R package from CRAN,³⁸ to estimate the mean ${\mu }^{\ast }$ and standard deviation ${\sigma }^{\ast }$of the empirical null distribution. Based on these estimates, we recompute the z-scores and corresponding p-values as follows

Finally, the resulting (empirical) p-values are corrected for multiple testing with the FDR method of Benjamini and Hochberg.³² As opposed to the original p-values that were calculated based on the theoretical null distribution for the t-statistics, we found that this procedure allows for a better FDR control in single-cell applications.

DTU tools literature

Table 1 provides a brief description of the DTU methods that were included in the performance benchmarks of this paper. For more details, we refer to the Extended data²⁵ and the respective original publications. Note that all methods were run in R3.6.1 using default settings.

Filtering

We adopted two different strategies for filtering transcripts in each of the RNA-seq datasets in the performance benchmarks.

The first filtering strategy uses the filterByExpr function implemented in edgeR.³⁹ This filtering strategy only retains transcripts that have at least an expression level of min.count counts-per-million (CPM, calculated as the number of read counts divided by the total number of reads in the dataset and multiplied by one million) in at least n samples or cells. In addition, the sum of the CPM of the transcript across all cells or samples must be at least min.total.count. For the bulk RNA-seq datasets, we use the default settings (min.count = 10, n = min(10, 0.7*n_s) and min.total.count = 10), with n_s the number of samples or cells in the smallest group. For the scRNA-seq datasets, the settings are adjusted to; min.count = 1 (as requiring a transcript to be expressed in all single-cells is a stringent criterium), n = 0.5*n_s and min.total.count = 0. In addition, if only one transcript of a gene passes this filtering criterion, it is omitted from the analysis, as DTU analysis is meaningless when only one transcript is retained. As such, we specifically set the parameters to generate a very lenient filtering criterium.

The second filtering strategy uses the dmFilter function implemented in DRIMSeq.²¹ This filter is more stringent and specifically designed for DTU analysis. The filtering process can be thought of as proceeding in three steps. The first step requires the transcripts to have a count of at least 10 in at least n_s samples. The second filtering step requires the transcript to make up at least 10% of the total count of its corresponding gene in at least n_s samples or cells. The third filtering step removes all transcripts for which the corresponding gene has a count below 10 in any of the samples or cells in the dataset. Again, if only one transcript of a gene passes this filtering criterion, it is omitted from the analysis.

Abundance metrics

We additionally compared the results between two quantification strategies: counts and scaledTPM abundances. Counts are the transcript-level abundances estimated by kallisto¹ or Salmon², i.e., the number of reads mapping to each transcript. scaledTPM normalized abundances are obtained by first computing transcript-per-million (TPM) abundances by (1) dividing the transcript counts by the length of the corresponding transcript in kilobases, (2) dividing by a cell/sample specific scaling factor, which is computed by taking the sum across transcripts of all length-corrected abundances from (1) and (3) multiplying by one million. Finally, the TPM abundances are scaled per cell/sample, so that the sum of the TPM abundances in the cell/sample equals the original number of counts for that cell/sample.

Bulk simulation study

To evaluate the performance of the different DTU analysis methods, we first adopt three simulated bulk RNA-seq datasets from previous publications: the simulated dataset from Love et al.¹⁸ (dataset 1) and both the Drosophila melanogaster (dataset 2) and Homo sapiens (dataset 3) simulation studies from Van den Berge et al..⁴⁰ All three datasets were generated based on parameter values obtained from real RNA-seq samples, to mimic real RNA-seq data as close as possible.

Notably, there is a subtle difference in how DTU is introduced between the two simulation frameworks. For dataset 1 from Love et al.,¹⁸ the origin of DTU is twofold: On the one hand, DTU was specifically introduced by swapping the transcript-per-million (TPM) abundances between two expressed isoforms. On the other hand, DTU was also obtained as a consequence of introducing DTE, where a single expressed isoform was induced to be differentially expressed at a certain log fold change, which leads to DTU if this transcript belongs to a gene expressing multiple isoforms. For datasets 2 and 3 from Van den Berge et al.,⁴⁰ there is only one source of DTU. The number of differentially used transcripts within a gene was sampled ranging from a minimum of two up to a random number drawn from a binomial distribution with size equal to the number of transcripts and success probability 1/3, thus allowing for differential usage of more than two transcripts within the same gene. For the selected transcripts, the percentage of usage withing the corresponding gene is first computed on the TPM scale. Next, DTU was introduced by swapping these percentages between transcripts. Subsequently, the swapped percentages are multiplied by the original gene-level read count per kilobase (RPK) to obtain transcript-level RPK. Finally, TPM and counts are computed based on these transcript-level RPK. Additionally, dataset 1 uses Salmon² (version 1.1.0) for estimating transcript-level abundances, whereas datasets 2 and 3 were quantified with kallisto¹ (version 0.46.0).

Real bulk RNA-seq dataset evaluation

We evaluate the performance of the different DTU methods on real bulk RNA-seq data, by subsampling a homogeneous set of samples from the large bulk RNA-seq dataset available from the Genotype-Tissue Expression (GTEx) consortium⁴¹ release version 8. Nine datasets were generated non-parametrically. More specifically, we first selected samples from adrenal gland tissue that were extracted with the RNA extraction method “RNA Extraction from Paxgene-derived Lysate Plate Based”. From the remaining samples we subsampled nine datasets, comprising three repeats for each of three sample sizes; 5 versus 5, 20 versus 20 and 50 versus 50 samples. Next, DTU is artificially in 15% of the genes. The number of differentially used transcripts within a gene was determined with the strategy of Van den Berge et al.⁴⁰ For the selected transcripts, the counts of transcripts within the same genes are swapped in one group of samples or cells, inducing DTU signal between the groups. The GTEx data was quantified with RSEM³⁵ version 1.3.0.

Real single-cell RNA-seq datasets evaluation

We evaluate the performance of the different DTU methods on real scRNA-seq datasets. These scRNA-seq datasets were generated non-parametrically by subsampling a homogeneous set of cells from three real scRNA-seq datasets,^26,29,42 after which DTU is artificially introduced as discussed above (see Methods, Real bulk RNA-seq dataset evaluation paragraph). The datasets from Chen et al.²⁹ and Tasic et al.⁴² were quantified using Salmon v1.1.0, whereas the dataset from Darmanis et al.²⁶ was quantified with Salmon v0.8.2 as obtained from Soneson et al.⁴³

For the dataset of Chen et al.,²⁹ which was used to construct Figures 4 and the Extended data Figure S7,²⁵ we selected a homogeneous population of cells by considering only the EpiStem cells of female mice, resulting in a dataset of 120 cells. From this homogeneous population of cells, we then subsampled six datasets, comprising three repeats for each of two sample sizes: 20 versus 20 and 50 versus 50 cells. Next, DTU was artificially introduced by swapping transcript counts between groups of cells. Finally, we adopted either edgeR or DRIMSeq for filtering.

The other two scRNA-seq datasets were generated analogously. For the dataset of Tasic et al.,⁴² which was used to construct Figure S8 (Extended data²⁵), we selected a homogeneous population of cells by considering only the Lamp5 cells in the anterior lateral motor cortex of mice without any eye conditions, resulting in a dataset of 897 cells. After introducing DTU by swapping transcript counts, we randomly subsampled 20, 75 or 200 cells from each group. For the dataset of Darmanis et al.,²⁶ which was used to construct Extended data²⁵ Figure S9, we selected the immune cells that clustered together in tSNE cluster 8 of the original publication, resulting in a dataset of 248 cells. After introducing DTU, we randomly subsampled 20, 50 or 100 cells from each group.

Case study differential gene expression analysis

We perform a differential gene expression (DGE) analysis on a subset of the Tasic single-cell dataset,⁴² i.e. between different the cell types originating from the ALM and VISp regions of the glutamatergic L5 IT subclass. We use the quasi-likelihood method of edgeR⁴⁴ to model the gene expression profiles and additionally adopt the edgeR glmTreat function to test differential expression against a log2-fold change threshold (log2-fold change = 1). Statistical significance was evaluated at the 5% FDR level.

Performance assessment

We assess the performance of different DTU methods on a bulk simulation dataset with scatterplots of the true positive rate (TPR) versus the false discovery rate (FDR), according to the following definitions:

where FN, FP and TP denote the numbers of false negatives, false positives and true positives, respectively. The FDR is the expected value of the false discovery proportion (FDP), which is the ratio of the number of false positives to the total number of positives. The FDR-TPR curves are constructed using the Bioconductor R package ICOBRA version 1.14.0.⁴⁵

Scalability benchmark

The scalability benchmark was run on subsets of the Chen scRNA-seq dataset,²⁹ which contains 617 cells in total. For the scalability benchmark with respect to the number of cells in the dataset, we randomly subsample a certain number of cells (8, 16, 32, 64, 128 or 256 cells per group) from the dataset (without introducing DTU or selecting specific homogeneous cell populations). Next, we filter this subsample using the edgeR-based filtering criterion. This was done to remove very lowly abundant transcripts, which may otherwise cause problems in the parameter estimation procedure. From the remaining transcripts, we randomly subsampled to a total of 30,000 transcripts before running the DTU analysis. To allow for a scalability benchmark of BANDITS, which scales poorly to the number of transcripts (Figure 5B), we also generated a dataset with only 1,000 transcripts (see Extended data figure S1²⁵). For the scalability benchmark with respect to the number of transcripts, we randomly sampled two groups of 16 cells from the dataset. After applying the edgeR-based filter, we sampled eight distinct numbers of transcripts: 1,000, 2,000, 5,000, 10,000, 15,000, 20,000, 25,000, 30,000 and 35,000 prior to the DTU analysis.

We additionally perform a scalability benchmark on the GTEx⁴¹ bulk RNA-seq dataset. Using the same strategy as described above for the scRNA-seq scalability benchmark, we here assess scalability with respect to the number of samples (8, 16, 32 or 64 samples per group) while keeping the number of transcripts fixed at 30,000 and scalability with respect to the number of transcripts (same range as above) while keeping the number of samples fixed at sixteen. Results are displayed in the Extended data²⁵ figures S20 and S21.

All scalability benchmarks were run on a single core of a virtual machine with an Intel(R) Xeon(R) CPU E5-2420 v2 (2.20GHz, Speed: 2200 MHz) processor and 30GB RAM.

Performance on simulated bulk RNA-seq datasets

To evaluate the performance of satuRn, we adopt three simulated bulk RNA-seq datasets from previous publications. Dataset 1 was obtained from Love et al.¹⁸ and contains two groups of 12 samples each, which we subsample without replacement to evaluate 3vs3, 6vs6 and 10vs10 two-group comparisons. Datasets 2 and 3 are the Drosophila melanogaster and Homo sapiens simulation studies from Van den Berge et al.⁴⁰ and Soneson et al.,¹⁷ which both contain two groups of five samples each. In brief, all datasets were constructed by generating sequencing reads based on parameters that are estimated from real bulk RNA-seq data. DTU between groups of samples was artificially introduced in the data, prior to the quantification of expression using either Salmon² (dataset 1) or kallisto¹ (datasets 2 and 3). Notably, there are some methodological differences between the simulation framework of dataset 1 and that of datasets 2 and 3 with respect to the read generation and the simulation of DTU signal (see Methods). In terms of transcript filtering, we adopt two different strategies as implemented by edgeR⁴⁴ and DRIMSeq,²¹ which correspond to a lenient and more stringent filtering, respectively (see Methods).

The result of the performance evaluation of satuRn with respect to other DTU methods on the three simulated bulk datasets is displayed in Figure 2. Figure 2A shows the average performance over three 6 versus 6 subsamples for dataset 1, after filtering with edgeR. Figures 2B and 2C display the performance on datasets 2 and 3 after edgeR filtering, respectively. In all three datasets, satuRn outperforms NBSplice, edgeR diffsplice and limma diffsplice. Intriguingly, the performance of DRIMSeq varies strongly between the three datasets. This discrepancy may be explained by the different strategies for generating reads and introducing DTU between dataset 1 on the one hand, and datasets 2 and 3 on the other hand (see Methods). We furthermore find the performance of satuRn is on par with the best performing tools from the literature, DEXSeq and DoubleExpSeq. In addition, both satuRn and DoubleExpSeq provide a stringent control of the FDR, while DEXSeq and DRIMSeq are often too liberal, as reported previously.¹⁸

Figure 2. Performance evaluation of satuRn on three simulated bulk RNA-Seq datasets.

Each curve visualizes the performance of each method by displaying the sensitivity of the method (true positive rate, TPR) with respect to the false discovery rate (FDR). The three circles on each curve represent working points when the FDR level is set at nominal levels of 1%, 5% and 10%, respectively. The circles are filled if the empirical FDR is equal or below the imposed FDR threshold. The performance of satuRn is on par with the best tools from the literature, DEXSeq and DoubleExpSeq, for all datasets. In addition, our method consistently controls the FDR close to its imposed nominal FDR threshold.

We also evaluated the effects of sample size and different filtering criteria on the performance of the different DTU methods (see Extended data²⁵ figures S2, S13). Neither sample size nor filtering criterion had a profound impact on the ranking of the performances of the different DTU methods; satuRn, DEXSeq and DoubleExpSeq remain the best performing methods overall. In addition, we studied the impact of using either transcript count estimates or normalized abundance estimates (scaledTPM, see Methods) as input data for the DTU algorithms. We observed a slightly higher performance in all datasets when providing count estimates, except for Dataset 1 from Love et al.¹⁸ Love et al.¹⁸ introduced DTU either by swapping TPM abundances or by introducing a fold change in transcript expression in one group of samples, again on TPM abundances. However, when back-transforming the swapped TPM abundances to count abundances, the transcript lengths were not swapped. Hence, the simulation strategy can induce differences in the total count for a gene, which in turn may alter the usage of all (target and off-target) transcripts in that gene when counts are used as an input, which gives rise to many false positives in off-target transcripts of DTU genes. When introducing DTU by swapping counts (real bulk and single-cell evaluations), we could expect spurious signal in analyses on the TPM scale for the same reason. However, the difference between the transcript count and TPM analyses was much less pronounced, which was also the case for the Van den Berge et al. simulation approach. Therefore, all performance evaluations in the body of this publication were generated with count estimates as input data, except for Figure 2, panel A. For a full overview on the effects of sample size, filtering criteria and data input type, we refer to Figures S2 and S13 of the Extended data.²⁵

Performance on a real bulk RNA-seq dataset

While simulation studies are common for evaluating the performance of DE analysis methods, there is currently no consensus on the simulation strategy that best mimics real (sc)RNA-seq data. In addition, simulation frameworks typically generate data according to parametric assumptions on the data-generating mechanism, thus potentially favoring DE methods that adopt similar distributional assumptions in their statistical model.⁴⁰ An alternative procedure is to non-parametrically modify a real dataset. Here, we obtained different subsamples from the large bulk RNA-seq dataset available from the Genotype-Tissue Expression (GTEx) consortium,⁴¹ generating nine datasets in total, i.e. three repeats for each of three sample sizes; 5 versus 5, 20 versus 20 and 50 versus 50 samples. We then artificially introduced DTU in these data by swapping transcript counts between groups of samples (see Methods for details). Again, we adopt two different filtering strategies as implemented by edgeR⁴⁴ and DRIMSeq²¹ (see Methods).

The results of the performance evaluation of satuRn on the real bulk datasets upon edgeR filtering is displayed in Figure 3. In agreement with the results obtained from the simulated bulk RNA-seq study, we observe that the performance of satuRn is on par with DEXSeq and DoubleExpSeq. Again, satuRn provides a conservative FDR control. While the FDR control of DoubleExpSeq is good overall, it appears to become too liberal with increasing sample size. In this evaluation, DRIMSeq performs poorly, in contrast to simulated bulk RNA-seq datasets 2 and 3, but in line with the performance evaluation on the simulated bulk RNA-seq dataset 1. Note that DEXSeq, DRIMSeq and NBSplice were omitted from the analysis of the largest dataset (50 versus 50 samples), as these methods do not scale to such large datasets (Figure 1). Adopting the DRIMSeq-based filtering did not have a qualitative impact on the performance (Extended data Figure S6²⁵).

Figure 3. Performance evaluation of satuRn on a real bulk RNA-Seq dataset.

Each curve visualizes the performance of each method by displaying the sensitivity of the method (true positive rate, TPR) with respect to the false discovery rate (FDR). The three circles on each curve represent working points when the FDR level is set at nominal levels of 1%, 5% and 10%, respectively. The circles are filled if the empirical FDR is equal or below the imposed FDR threshold. The performance of satuRn is on par with the best tools from the literature, DEXSeq and DoubleExpSeq. In addition, satuRn consistently controls the FDR close to its imposed nominal FDR threshold, while DoubleExpSeq becomes more liberal with increasing sample sizes. Note that DEXSeq, DRIMSeq and NBSplice were omitted from the larger comparison, as these methods do not scale to large datasets (Figure 1).

Performance on real single-cell data

Finally, we evaluate the performance of satuRn on single-cell RNA-seq data. As with the real bulk analysis, the single-cell datasets were generated by subsetting from three different real scRNA-seq datasets^26,29,42 (see Methods). Again, we subsampled three repeats of different sample sizes, artificially introduced DTU with the swapping strategy and applied either the edgeR- or DRIMSeq-based filtering criterium (see Methods for details).

By subsampling the Chen et al.²⁹ dataset, we generated three repeats of two sample sizes, i.e. 20 versus 20 and 50 versus 50 cells. The results of the performance evaluation of satuRn on this dataset upon edgeR filtering is displayed in Figure 4. The performance of satuRn is slightly better than that of the best tool from the literature, DoubleExpSeq. As compared to the evaluations on bulk data, we observe a performance drop for DEXSeq relative to satuRn and DoubleExpSeq. This, in combination with its poor scalability (Figure 1), greatly compromises the use of DEXSeq for the analysis of scRNA-seq data. satuRn again provides a stringent control of the FDR, while the inference of DoubleExpSeq is too liberal, again becoming more problematic for larger sample sizes. Adopting the DRIMSeq filter did not have a qualitative impact on the performances (Extended data figure S7²⁵). The results of the performance evaluations on the other two scRNA-seq datasets^26,42 are in strong agreement with the results displayed here, with satuRn performing at least on par with DoubleExpSeq and satuRn additionally controlling the FDR around the nominal level (Extended data²⁵ figures S8 and S9).

Figure 4. Performance evaluation of satuRn on a real scRNA-Seq dataset.

Each curve visualizes the performance of each method by displaying the sensitivity of the method (true positive rate, TPR) with respect to the false discovery rate (FDR). The three circles on each curve represent working points when the FDR level is set at nominal levels of 1%, 5% and 10%, respectively. The circles are filled if the empirical FDR is equal or below the imposed FDR threshold. The performance of satuRn is on par with the best tools from the literature, DEXSeq and DoubleExpSeq. In addition, our method provides a stringent control of the FDR, while DoubleExpSeq becomes more liberal with increasing sample sizes (see also Extended data figure S7²⁵). Note that DEXSeq and DRIMSeq were omitted from the two largest comparisons, as these methods do not scale to large datasets (Figure 1). NBSplice was omitted from all comparisons, as it does not converge on datasets with many zeros, such as scRNA-Seq datasets.

Notably, we found that the theoretical null distribution of the test statistics from satuRn failed to provide good FDR control in single-cell analyses (Extended data figure S14²⁵). To obtain proper p-values with satuRn in single-cell applications, we therefore estimate the null distribution of the test statistic empirically (see Methods, satuRn paragraph). Note that the performances in the bulk RNA-seq benchmarks are not affected by adopting the empirical null distribution, because no deviations of the theoretical null distribution occur in bulk data. However, the empirical null resulted in much improved FDR control in scRNA-seq datasets (Extended data²⁵ figure S14). We therefore adopt the empirical null estimation as the default setting in satuRn. As such, all satuRn results in this publication are relying on the empirical null strategy. As a final remark, we likewise attempted to improve the FDR control of DoubleExpSeq. However, in all analyses with DoubleExpSeq we observed a large spike of p-values equal to 1, which poses a problem when estimating the empirical null distribution (Extended data²⁵ figure S15). Therefore, this strategy could not be used to improve the FDR control of DoubleExpSeq.

In these benchmark single-cell datasets, DTU signal was introduced non-parametrically by swapping transcript counts (see Methods), thus introducing DTU signal of different magnitude in the different transcripts. In Extended data²⁵ Figures S16-S19, we stratified the FDR-TPR curves on the difference in the observed average transcript usage between the two groups of cells. In all benchmark datasets, satuRn and DoubleExpSeq are more successful in detecting small differences as compared to the other methods. Indeed, while all methods succeed in detecting large differences in absolute transcript usage, the performance gain of satuRn and DoubleExpSeq in the results without stratification can be attributed to picking up DTU with smaller effect sizes.

Scalability benchmark

We performed a computational benchmark of satuRn to investigate its scalability with respect to the number of samples/cells and the number of transcripts in an RNA-seq dataset. All scalability benchmarks were run on a single core of a Linux machine with an Intel(R) Xeon(R) CPU E5-2420 v2 (2.20GHz, Speed: 2200 MHz) processor and 30GB RAM. The results are displayed in Figure 5.

Figure 5. Scalability evaluation on scRNA-Seq data. A: Runtime with respect to the number of cells in a scRNA-Seq dataset.

Left panel: DRIMSeq and especially DEXSeq scale poorly with the number of cells in the dataset. Right panel: Detailed plot of the fastest methods. satuRn scales linearly with increasing numbers of cells, with a slope that is comparable to that of limma diffsplice. As such, satuRn can perform a DTU analysis on a dataset with two groups of 256 cells each and 30,000 transcripts in less than three minutes. For all sample sizes, the number of transcripts in the datasets were set at 30,000. Note that NBSplice was not included in this analysis as it fails to converge on datasets with a large proportion of zero counts. B: Runtime with respect to the number of transcripts in a scRNA-Seq dataset. Left panel: DEXSeq, DRIMSeq and especially BANDITS scale poorly to the number of transcripts in the dataset. Right panel: Detailed plot of the fastest methods. satuRn scales linearly with increasing numbers of transcripts, but with a steeper slope than edgeR diffsplice, DoubleExpSeq and limma diffsplice. The number of cells in the dataset was set fixed to two groups of 16 cells. All scalability benchmarks were run on a single core.

Figure 5A displays the scalability with respect to the number of cells in the dataset, (different subsets of the single-cell dataset of Chen et al.²⁹), while keeping the number of transcripts in the dataset fixed at 30,000. From the left panel, it is clear that DRIMSeq and especially DEXSeq scale very poorly with the number of cells in the dataset, which was already shown in Figure 1B. In the right panel, we focus on the four remaining methods. satuRn scales linearly with increasing numbers of cells, with a slope comparable to limma diffsplice. As such, satuRn can perform a DTU analysis on a dataset with two groups of 256 cells each and 30,000 transcripts in less than three minutes. Note that BANDITS³⁰ was not included in this benchmark, as it does not scale to datasets with this many transcripts. For a performance and scalability evaluation of BANDITS on datasets with a lower number of transcripts, we refer to Extended data figure S1.²⁵ NBSplice was also omitted as it fails to converge on datasets with a large proportion of zero counts.

Figure 5B shows the scalability with respect to the number of transcripts in the dataset, while keeping the number of cells in the dataset fixed to two groups of 16 cells. As shown in Figure 1C, BANDITS, DEXSeq and DRIMSeq scale poorly to datasets with many transcripts. From the right panel, satuRn scales linearly with increasing numbers of transcripts, albeit with a steeper slope than edgeR diffsplice, DoubleExpSeq and limma diffsplice. Note, that the scalability of DTU analyses can be improved through parallelization, if this is allowed by the underlying algorithm. Parallel execution is possible in satuRn and in all methods from the literature that were discussed in this manuscript, except for DoubleExpSeq and NBSplice.

In Extended data Figure S20,²⁵ we additionally evaluate the scalability of the DTU tools on different subsets of the bulk GTEx⁴¹ dataset. These results are in line with those of the single-cell scalability benchmark, with DEXSeq and DRIMSeq being markedly slower than the other methods. When directly comparing the scalability between the bulk and single-cell data on equally sized subsets (Extended data Figure S21²⁵), we found limited differences in scalability between the data types, except for DEXSeq. DEXSeq scales considerably better on bulk data than on single-cell data, suggesting that the estimation of the GLM parameters is slower when the data are sparse. However, as the scalability profile of DEXSeq is quadratic with respect to the number of cells/samples in the data, it is still infeasible to adopt DEXSeq in datasets with many cells/samples.

Altogether, we find that while several methods for DTU analysis exist, none are optimally suited for analyzing single-cell datasets. DRIMSeq, NBSplice, edgeR diffsplice and limma diffsplice have a much lower overall performance in our benchmarks. DEXSeq does not scale to large datasets. Finally, DoubleExpSeq does not support experimental designs that require an analysis with multiple additive effects, e.g. randomized complete block designs and designs where batch-effect correction is required, which are essential for many practical scRNA-Seq analysis settings.⁴⁶ In addition, it fails to control the FDR at the desired level, especially with increasing sample sizes. With satuRn, we overcome these issues and effectively unlock DTU analysis for single cell applications.

Case study

We use satuRn to perform a DTU analysis on a subset of the single-cell (SMART-seq2¹¹) RNA-seq dataset from Tasic et al.⁴² In addition, we analyze the same dataset with DoubleExpSeq and limma diffsplice, which are the only other DTU methods that scale to large scRNA-seq datasets and have a reasonable performance in our benchmarks. In the original publication, the authors studied differential gene expression between cell types originating from two areas at distant poles of the mouse neocortex; the primary visual cortical area (VISp), which processes sensory information with millisecond timescale dynamics^47–49 and the anterior lateral motor cortex (ALM), which displays slower dynamics related to short-term memory, deliberation, decision-making and planning.^50,51 Based on marker genes, Tasic et al.⁴² assigned all of the 23,822 cells from the scRNA-seq dataset to one of three cell classes; glutamatergic (excitatory) neurons, GABAergic (inhibitory) neurons or non-neuronal cells. The authors then further classified the neuronal cells into several subclasses based on their dominant layer of dissection and projection patterns (through a retrograde labelling experiment). Finally, these subclasses are further classified into cell types based on the expression of specific marker genes.

DGE analysis with edgeR

In their original DGE analysis, Tasic et al.⁴² obtained the largest number of differentially expressed genes between the cell types originating from the ALM and VISp regions of the glutamatergic L5 IT subclass (2,739 cells in total), where L5 refers to layer-of-dissection 5 and IT refers to the intratelencephalic projection type. Here, we first perform a DGE analysis with an edgeR-based workflow (see Methods) on the same comparisons between L5 IT cell types that were assessed by Tasic et al. Table 2 shows the number of differentially expressed genes between the groups of interest in column 4.

DTU analysis with satuRn

Next, we perform a DTU analysis for the same cell types using satuRn. In column 5 of Table 2, we display the number of differentially used transcripts for each comparison. We also show the number of unique genes in which we find evidence for changes in usage of at least one transcript (column 6). While the number of differentially used transcripts is lower than the number of differentially expressed genes in each of the contrasts, we did identify differentially used transcripts in all contrasts of interest. Most interestingly, we observe that the overlap between the differentially expressed genes and the genes in which we found evidence for DTU is very limited (Table 2, column 7). This shows that the information obtained from our DTU analyses are orthogonal to the results from the canonical DGE analyses, which has been reported previously for simulated bulk data.¹⁸

Gene set enrichment analysis

We perform a gene set enrichment analysis (GSEA⁵²) on the three comparisons with most DE genes and most genes with evidence for DTU (comparisons 5, 6 and 7). Similar gene ontology categories are returned for the set of DGE genes and the set of DTU genes, with many of the enriched processes being biologically relevant in the context of this case study. Enriched gene sets include the gene ontology classes, synapse, neuron projection, synaptic signaling and cell projection organization. This shows that the complementary information brought by the DTU analysis is indeed biologically relevant. For an extensive overview of the GSEA of the set of DGE genes and genes with evidence of DTU in comparisons 5, 6 and 7, we refer to the Extended data.²⁵

satuRn identifies biologically relevant DTU transcripts

To display the utility of satuRn for identifying and visualizing DTU transcripts in scRNA-seq datasets, we focus on comparison #6 of the DTU analysis and discuss the gene P2X Purinoceptor 4 or P2rx4 (Ensembl ID ENSMUSG00000029470), alsb gene which is part of a family of purinergic receptors that have been implicated in functions such as learning, memory and sleep. In the DGE analysis, no evidence for differential expression of P2rx4 was found at the gene level (FDR-adjusted p-value = 1). By contrast, in our DTU analysis the transcripts of P2xr4 displayed the highest statistical evidence for differential usage within the set of transcripts that could be assigned to the ontology class ‘neuron projection’.⁵² The mean usage of transcript ENSMUST00000081554 is estimated to be 28.9% in Tnc cells and 75.9% in Hsd11b1 Endou cells (FDR-adjusted p-value = 1.22E-13). For transcript ENSMUST00000195963, the transcript usage changes from 58.2% in Tnc cells and 16.6% in Hsd11b1 Endou cells (FDR-adjusted p-value = 1.79E-10). For the third transcript of P2rx4 that was assessed in our DTU analysis, ENSMUST00000132062, we found no statistical evidence for DTU (FDR-adjusted p-value = 0.534). In Figure 6, we show the output for the visualization of the transcript usages for P2rx4 as obtained from satuRn. Interestingly, the majority transcript in the Tnc cell type, ENSMUST00000195963, is not protein coding.⁵⁴ By contrast, the majority transcript in the Hsd11b1 Endou cell type, ENSMUST00000081554, is coding for the P2X purinoceptor protein (UniProt ID Q9Z256). As such, the changes in transcript usage between both cell types represent actual biological differences in the functionality of the gene products, which may be relevant to the process of neuron projection. This functional difference would have remained obscured when only performing a canonical DGE analysis.

Figure 6. Differential transcript usage in the P2rx4 gene.

Panels A-C show the usage of the different transcripts of the P2rx4 gene across cells of the Tnc and Hsd11b1 cell types. The size of each datapoint is weighted according to the total expression of the gene in that cell, i.e., the gene counts per cell. The yellow diamonds indicate the estimated mean usage of a transcript for each cell type, as estimated by satuRn. The cyan and dark green diamonds indicate mean and median gene expression levels per cell type, respectively. Panels A and B display the transcript usage across cells of the Tnc and Hsd11b1 Endou cell types of transcripts ENSMUST00000081554 and ENSMUST00000195963, respectively. The proportion of usage of the former transcript is clearly higher in Hsd11b1 Endou cells, while the latter transcripts is most abundant in Tnc cells. For the third transcript, ENSMUST00000132062 (panel C) there is no evidence for differential usage between both cell types. In addition, there is no evidence for differential expression of P2rx4 on the gene level (panel D). DTU and DGE significance levels are indicated in the figure headers.

Comparison to limma diffsplice

We also analyzed the case study dataset with limma diffsplice.²³ When running limma diffsplice with default settings, a large number of DTU transcripts was returned (Extended data²⁵ figure S22) and we observe that the p-values were shifted towards smaller values (Extended data²⁵ figures S23 and S24). Therefore, we adopted the same empirical null strategy as implemented in satuRn to post-process the results. While this dramatically decreased the number of significant DTU transcripts, limma diffsplice still identified more transcripts (i.e., true or false positives) than our method. However, when we inspected the transcripts that were highly ranked in the top DTU list of limma diffsplice but lowly ranked in our top list, we found that most of these transcripts either originate from genes that are lowly expressed, or they are transcripts with a large fraction of zero counts (i.e., zero expression in a large percentage of cells). Limma diffsplice thus claims differential usage more often for transcripts that only contain little information for assessing DTU. This is depicted in Figure 7.

Figure 7. Global comparison between differential transcript usage (DTU) transcripts uniquely identified by satuRn, uniquely identified by limma diffsplice or by both methods.

Left panel: Boxplots on the average gene-level count for the DTU genes identified by the respective methods. Transcripts uniquely identified by satuRn originate from genes that have a much higher gene-level count (averaged over cells) as compared to transcripts uniquely identified by limma diffsplice. Note that the y-axis is displayed on a log10 scale. Right panel: Violin plots indicating the fraction of cells in which the transcripts are expressed. Transcripts uniquely identified by satuRn are expressed, on average, in a much larger fraction of the cells. Conversely, transcripts identified as DTU uniquely by limma diffsplice often have no expression in a large fraction of the cells. The dark green diamond indicates the median fraction of cells in which the DTU transcripts are expressed.

This behavior can be expected. Limma diffsplice tests for DTU by comparing the log-fold change in expression of transcript t with the average log-fold change in the expression of all transcripts belonging to the same gene as transcript t. As such, limma diffsplice does not incorporate any information on the absolute gene expression levels. In contrast, our quasi-binomial GLM framework models the log-odds of drawing a particular transcript t from the pool of transcripts in the corresponding gene. As a consequence, transcripts belonging to lowly expressed genes are correctly considered less informative in satuRn and are thus less likely to be picked up. For example, in Figure 8A, we show that while our method estimates a mean usage of 7% in Tnc cells and 26% in Hsd11b1 Endou cells (indicated by the gold diamond), the transcript is not identified as differentially used, given the low abundance of the corresponding gene and the highly variable single-cell level observations.

Figure 8. Three examples displaying differential transcript usage (DTU) transcripts that are uniquely identified by satuRn or limma diffsplice.

Each panel shows transcript usage across cells of the Tnc and Hsd11b1 cell types. The size of each datapoint is weighted according to the total expression of the corresponding gene in that cell, i.e., the total gene count per cell. The yellow diamonds indicate the estimated mean usage of a transcript for each cell type, as estimated by satuRn. The cyan diamonds indicate the mean transcript expression levels per cell type. The header of each panel indicates the FDR-adjusted p-value and the rank of the DTU finding in the top lists by limma diffsplice and satuRn analyses. Panel A: Transcript uniquely identified as differentially used by limma diffsplice. The DTU claim by limma is driven by the difference in mean transcript usage between cell types. Given the low abundance of the corresponding gene and the highly dispersed single-cell level observations, satuRn does not identify the transcript as differentially used. Panel B: Transcript uniquely identified as differentially used by satuRn. Even though the mean difference in transcript usage between cell types is estimated to be 3%, satuRn claims significance given that the difference is stably supported by many cells with high gene-level expression levels. Panel C: Transcript uniquely identified as differentially used by limma diffsplice. The DTU claim by limma is driven by the raw mean difference in transcript usage between cell types. In contrast, satuRn takes into account that the Hsd11b1 Endou cells expressing the transcript at 0% usage have low gene-level count. The size of the dots (which represent individual cells) is weighted according to the total expression of the gene in that cell, i.e., the total gene count per cell. The yellow diamonds indicate the estimated mean usage of a transcript for each cell type, as estimated by satuRn. The cyan diamonds indicate the raw mean transcript usage levels per cell type.

Conversely, by looking at the transcripts that were highly ranked in our DTU list but lowly ranked in the top list of limma, we observe that our model is more likely to capture small changes in transcript usage that are stable across all cells and belong to genes that are highly expressed. An example of such a transcript is shown in Figure 8B. satuRn estimates a mean usage of 3% in Tnc cells and 6% in Hsd11b1 Endou cells. While this is only a minor change in transcript usage, satuRn still identifies this transcript as differentially used because the gene is highly expressed and the small change in usage is supported by a large number of cells. In case such small differences in usage are not considered biologically meaningful, it is possible to set a threshold on the minimal desired difference. Finally, by not taking into account gene abundances, limma is more influenced by outlying observations that have a low gene-level abundance (Figure 8C). Indeed, DTU claims by limma are driven by differences in raw mean usages of transcripts. In Figure 8C, the raw mean usage of the transcript is 77% in Tnc cells and 45% in Hsd11b1 Endou cells, as indicated by the cyan diamonds. By contrast, the mean usage estimate by satuRn, which takes into account that the Hsd11b1 Endou cells expressing the transcript at 0% usage have low gene-level count, is 83% for Tnc cells and 75% for Hsd11b1 Endou cells, as indicated by the gold diamonds.

We therefore argue that, given the above observations, the transcripts identified by satuRn should be considered more reliable, as they generally originate from genes containing more information for assessing DTU.

Comparison to DoubleExpSeq

We additionally analyzed the dataset by Tasic et al. with DoubleExpSeq.²⁰ DoubleExpSeq identified a large number of DTU transcripts in all eight comparisons between cell types, ranging from 335 to 4580 DTU transcripts (Extended data²⁵ figure S22). This is consistent with our performance benchmarks, which already suggested that DoubleExpSeq becomes overly liberal in single-cell datasets with a large number of cells (Figure 4, Extended data²⁵ Figures S7, S8 and S9). We therefore expect many of these transcripts to correspond to false positives. Furthermore, this is reflected in the pathological distribution of p-values obtained by DoubleExpSeq, where p-values tend to be small and therefore the analysis too liberal (Extended data²⁵ figure S25). Furthermore, as discussed in the benchmark studies, we could not adopt the empirical null strategy to improve the FDR control of DoubleExpSeq. Again, a large number of p-values equal 1 poses a problem for estimating the empirical null distribution (Extended data²⁵ figure S26).

While the results of DoubleExpSeq are likely to be overly liberal, the ranking of the transcripts (based on the p-values of the DTU analysis) might still be reasonable. In Figure 9 we observe a large overlap between the top 200 transcripts identified by satuRn in comparison #6 of the case study and the top 200 transcripts of DoubleExpSeq in that comparison. This overlap is considerably smaller with a limma diffsplice analysis.

Figure 9. Venn diagram displaying the degree of overlap of the top 200 transcripts in comparison #6 of the case study in three DTU analysis tools.

We observe that in the set of the top 200 transcripts identified by satuRn, 149 transcripts overlap with the top 200 list from DoubleExpSeq. In the top 200 list of limma diffsplice, 108 transcripts are present that were not in the top lists of satuRn or DoubleExpSeq.

Finally, we note that while DoubleExpSeq could still be used in this case study given the simple factorial design (using a single factor to assign each cell to a cell type), DoubleExpSeq cannot be used in multifactorial designs, for instance to compare expression levels across multiple cell types between multiple samples or treatment groups.

Analyses with complex experimental designs

The analysis of scRNA-seq experiments often requires incorporating covariates in the statistical model, for instance to correct for batch effects or to account for subject-level baseline characteristics like age or body mass index. satuRn and most other DTU tools allow for the analysis of complex experimental designs with multiple categorical and continuous covariates (see Table 1). Exceptions are DoubleExpSeq and BANDITS that do not support modeling multiple covariates, which limits their utility for scRNA-seq data analysis and bulk RNA-seq experiments with complex designs. This may result in suboptimal data analysis. First, covariates may be able to explain a larger part of the variability in the data, which will increase the statistical power to pick up the effect of interest, as the explained variability would otherwise be incorporated in the error term. Second, the omission of covariates may result in biased parameter estimates. Third, the included covariates themselves could be key to further unravel the studied biological process.

To demonstrate satuRn’s ability to analyze datasets with complex experimental designs, we reanalyzed the single-cell study from Tasic et al.⁴² but now we include a main effect for cell type, gender, and the continuous covariate age, as well as an interaction effect between cell type and gender.

We first assess DTU between the Hsd11b1 Endou and Tnc cell types and compare the results with the original analysis ignoring gender and age effects. As expected, the additional covariates affect the ranking of transcripts in the resulting top list. For instance, transcript ENSMUST00000198199 dropped from the 8^th to the 68^th position after correcting for gender and age (empFDR = 1.89e-11 without covariates, empFDR = 2.26e-3 with covariates). Second, we assessed the interaction effect to prioritize transcripts for which the DTU between cell types differs according to gender. The differential transcript usage of ENSMUST00000145174 between Hsd11b1 Endou and Tnc cell types, for instance, is significantly different between female and male mice (empirical FDR = 0.0036). Indeed, for female mice there is evidence for DTU (empirical FDR = 0.01258), while no difference in usage is observed for male mice (empirical FDR = 0.9830, Figure 10).

Figure 10. DTU between Hsd11b1_Endou cells and Tnc cells of transcript ENSMUST00000145174 is different between female and male mice.

The difference in the usage of transcript ENSMUST00000145174 between Hsd11b1 cells and Tnc cells is significantly different between female and male mice (empirical FDR = 0.0036). Left panel: In female mice, the usage the transcript is significantly higher in Hsd11b1_Endou cells than in Tnc cells (empirical FDR = 0.0126). Right panel: In male mice, no difference in usage is observed (empirical FDR = 0.9830).

Alternative input data types

Above, we have used satuRn for transcript-level differential usage analyses. However, there are two other differential usage analyses that we can perform at the sub gene level: differential usage of transcript compatibility counts,⁵⁵ sometimes also referred to as equivalence class counts,⁵⁶ and differential exon usage analyses. Below, we first repeat the case study on the single-cell dataset from Tasic et al.⁴² using a transcript compatibility count matrix as input to satuRn. Second, we perform a differential exon usage analysis on a subset of the “pasilla” bulk RNA-Seq dataset by Brooks et al.,⁵⁷ which is the dataset that was used in the exon-level differential usage analysis of the DEXSeq publication.¹⁹

Case study on transcript compatibility counts

To obtain transcript-level abundance estimates, pseudoalignment tools like kallisto¹ or Salmon² rely on an expectation-maximization (EM) algorithm⁵⁸ or similar for assigning reads to transcripts. This is a challenging task, given the largely overlapping RNA sequence of different gene isoforms. Therefore, it has been suggested to perform differential usage analyses on transcript compatibility counts (TCCs),^55,59,60 sometimes referred to as equivalence class (EC) counts,⁵⁶ instead. An equivalence class is the set of transcripts that a sequencing read is compatible with. As such, reads that are pseudo-aligned to the same set of transcripts are part of the same equivalence class. A TCC is the number of reads assigned to each of the different equivalence classes. Note, that TCCs can be obtained much faster than transcript-level counts, as the former does not require running the EM algorithm. Another important advantage of TCCs is that they can also be obtained from UMI-based scRNA-seq experiments, like the 10× Chromium protocol.

Here we perform a differential equivalence class usage (DECU) analysis on the single-cell dataset from Tasic et al.⁴² The DECU analysis is performed without correction for gender and age to allow for a comparison with our elaborate transcript-level analysis from the first case study analysis. We adopt the same filtering criteria as for the transcript-level analysis. However, ECs that correspond to more than one gene are additionally removed, as the relationship between ECs and genes needs to be determined unambiguously to infer within-gene EC usages.^56,59 In our analysis, approximately 15% of the reads were discarded for this reason, which is in line with previous reports.⁶¹ After these steps, 26.141 ECs are retained, as compared to 14.420 transcripts in the transcript-level analysis. Next, the same comparisons of interest are tested. The number of differentially used ECs for each comparison are displayed in Table 3.

For several comparisons, there is a clear discrepancy between the number of differentially used transcripts and equivalence classes (Table 3, columns 4 and 6) and in the number of unique genes in which evidence for differential usage was observed (Table 3, columns 5 and 7). The overlap between the genes found in the transcript-level and the EC-level analysis ranges from 53% in comparisons 6 and 8 to 100% in comparison 3.

For comparison 6 (53% overlap), we investigate the genes identified by both analyses. For some genes that are identified in both analyses, there is a strong agreement between the transcript-level and EC-level results. One example is gene Phosphodiesterase 4D or Pde4d (Ensembl ID ENSMUSG00000021699). As depicted in Figure 11, there is a one-to-one concordance in the usage profiles of the equivalence classes and the transcripts. Interestingly, the first equivalence class 52911|52918 has reads that are compatible to both transcripts ENSMUST00000134973 and ENSMUST00000138938. In the transcript-level analysis, transcript ENSMUST00000138938 was not retained after filtering, as most of these reads were assigned to transcript ENSMUST00000134973 by Salmon. Note that while the interpretation at the transcript level is more sensible from a biological perspective, it also assumes an appropriate attribution of reads to transcripts.

Figure 11. Concordance between a differential usage analysis at the equivalence class level and the transcript level for gene Pde4d.

Top row: Equivalence class-level analysis. Three equivalence classes passed feature-level filtering. For equivalence class 52911|52918, strong evidence for differential usage between cell type Hsd11b1_Endou and Tnc was observed (empirical FDR = 3.970e-05). Bottom row: Transcript-level analysis. Three transcripts passed feature-level filtering. Each transcript corresponds to the respective equivalence classes in the top row. Note that while equivalence class 52911|52918 is compatible with both transcripts ENSMUST00000134973 and ENSMUST00000138938, the latter did not pass filtering in the transcript-level analysis. The visualizations and empirical FDR values are strongly concordant between the EC-level and the transcript-level analyses. Note that the gene-level counts (size black dots) in the EC-level analysis are lower because equivalence classes that correspond to multiple genes are being removed.

For most other genes in which differential usage was identified in both analyses, the comparison between the transcript analysis and the equivalence class analysis is more ambiguous, as equivalence classes often do not allow to pinpoint the responsible transcript(s) driving the observed differential signal. For example, the gene P2X Purinoceptor 4 or P2rx4 (Ensembl ID ENSMUSG00000029470), which was discussed in Figure 6 of the transcript-level case study, was also identified in the EC-level analysis. Four P2rx4 equivalence classes passed feature-level filtering, which we will denote with EC1-EC4. EC1 uniquely corresponds to transcript ENSMUST00000195963, and both the transcript and EC are found to be differentially used in their respective analyses (Extended data Figure S27, panel A²⁵). The other ECs correspond to two, three and four different transcripts, respectively, but not to transcript ENSMUST00000195963 (Extended data Figure S27, panel A²⁵). For EC2 and EC3, strong evidence for differential usage was found. For EC4, there was no evidence for differential usage. However, based on these ECs alone, it would not be possible to unambiguously attribute these changes to single transcripts, complicating the biological interpretation (Extended data Figure S27²⁵). In contrast, while the transcript-level analyses do allow for clear functional interpretation, whether these results reflect true underlying biology depends on the ability of the quantification tools to assign reads to transcripts (see Discussion). Finally, none of the four ECs are compatible with transcript ENSMUST00000132062, the third transcript that passed feature-level filtering in the transcript-level analysis. In the EC-level analysis, there were many equivalence classes corresponding to this transcript but with a relatively low count, such that these individual equivalence classes did not pass filtering. The fact that the transcript was retained in the transcript-level analysis can be explained by Salmon assigning many of the reads of these low count equivalence classes that correspond to multiple transcripts to transcript ENSMUST00000132062.

Case study on exon-level counts

Finally, we performed a demonstrational case study on a subset of the Drosophila melanogaster “pasilla” bulk RNA-Seq dataset by Brooks et al.,⁵⁷ which is the dataset that was used in the exon-level differential usage analysis of the DEXSeq publication and its corresponding R package vignette.¹⁹ In brief, this dataset consists of two experimental groups of three and four samples each. In addition, these samples were sequenced partly in single-end runs and partly in paired-end runs, which is corrected for by including a covariate to the DEXSeq and satuRn models. The dataset was obtained from the pasilla R data package,⁶² which contains a count matrix for the seven samples and 498 exons originating from 46 genes.

Compared to the original DEXSeq analysis, satuRn identified fewer differentially used exons between the two experimental groups, with three statistically significant findings at 5% FDR for satuRn compared to 20 significant findings with DEXSeq. While we did not benchmark satuRn and DEXSeq on exon-level data, this could correspond to DEXSeq being overly liberal in transcript-level bulk analyses, especially with small sample sizes (Figure 2 and Extended data Figure S2²⁵). More importantly, however, the ranking of the exons based on their respective p-values is almost identical between the two methods (Extended data Figure S28²⁵). Again, this is supportive of a similar performance of DEXSeq and satuRn on small bulk RNA-Seq datasets, albeit with a different type 1 error rate. In Extended data Figure S29,²⁵ we additionally demonstrate how satuRn can be used to visualize changes in the usage of exons.

Underlying data

Zenodo: Datasets associated with this publication https://doi.org/10.5281/zenodo.6826603.⁶⁵

This project contains the following underlying data:

In addition, all folders except Scalability_analysis.zip contain intermediate DTU analysis results that are available as.Rdata files or, in the case of Case_study.zip, .Rds files. Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).

Extended data

Zenodo: Extended data for satuRn publication https://doi.org/10.5281/zenodo.6810615.²⁵

This project contains the following extended data:

License: Data are available under the terms of the CC-BY 4.0 license.