Selecting differential splicing methods: Practical considerations

Introduction

Alternative splicing (AS) can be best described as fine-tuning gene expression by rearranging exons and introns in pre-mRNA. With 90-95% of human multi-exon genes estimated to possess some form of alternative splicing, it is a widespread regulatory process in cellular biology.¹ The cell utilises a large ribonucleoprotein (RBP) complex known as the spliceosome which is guided to target sites through the interaction of sequence elements (splice sites, enhancers & silencers and the polypyrimidine tract) and/or splicing factors. Pre-mRNA splicing can also occur without the splicesome as in the case of self-splicing group I & II introns, tRNA splicing and trans-splicing.² This ultimately results in genome-wide transcript diversity and subsequently, measurable changes to protein functionality.

Previous research has uncovered the phenotypic consequences of alternative splicing in disease. In humans, clinical research has shown alternative splicing (AS) as a key instigator in several forms of cancer and neurodegenerative disorders.^3–5 One notable discovery in Microtubule-associated protein tau’s (MAPT) possession of mis-spliced isoforms causing abnormal TAU accumulation progressing to Alzheimer’s disease.⁶ In cancer, numerous mis-spliced variants of tumour suppressors, apoptotic and angiogenic proteins have been discovered to contribute to tumour progression.^7,8 Beyond clinical research, the utility of alternative transcripts for bioengineering purposes has been explored. For example, an alternatively spliced version of the transcription factor X-box binding protein 1 (XBP1) coexpressed in production cell lines has been shown to increase productivity in the biomanufacturing of recombinant proteins.^9–11 In bio-agriculture, the CRISPR-mediated directed evolution of SF3B1 mutants (a spliceosomal component) in rice has improved crop traits through better resistance to splicing inhibitors.¹² Increasingly, the value of AS in both clinical and biotechnology applications has been recognised; highlighting the need for robust bioinformatics pipelines to identify variants.

For prospective researchers to investigate AS, the transcriptomic data is usually generated using next-generation sequencing. Short-read RNAseq is the most commonly used experimental technique to interrogate a transcriptome owing to its versatility and cost-effectiveness.^13,14 It involves sequencing short fragments of RNA molecules, providing insights into the respective expression levels of genomic features assembled from reference genomes. These features may be coding sequences, genes, transcripts, exons, introns, codons or even untranslated regions. A typical RNAseq pre-processing pipeline will consist of quality control (QC), read alignment & quantification before statistical analysis begins. QC assesses the quality of the raw fragmented reads using a standardised tool such as FastQC and trims low-quality reads or adaptor sequences.¹⁵ Then for alignment, a reference genome/transcriptome arranges the subsequent sequences into feature bins such as genes, transcripts, exons and coding sequences using software such as STAR or HISAT.^16,17 Alignment files (usually in the form of Sequence Alignment Maps: SAMs) can then be quantified to these features using a quantification tool such as HTSeq, Salmon or featureCounts usually normalising for library size and sequencing depth.^18–20 Depending on the purpose of analysis, normalisation may be scaled by total number of reads (CPM: Counts per Million), per length of transcript (TPM: Transcripts per Million), by paired-end fragments (RPKM: Fragments Per Kilobase of Transcript) or by using a median of ratios (DESeq2’s method).²¹ Commonly, a differential expression analysis will be performed at the gene or transcript level between groups of samples to identify statistically significant changes in expression. The pre-processing steps for RNA-seq have been extensively researched over many years, and there is a consensus within the community regarding the gold-standard set of tools. Projects like nf-core enable the execution of RNA-seq pre-processing pipelines with minimal intervention and limited bioinformatics expertise.²² However, these tend to be focused on the use-case of conventional differential expression rather than the more bespoke AS pipelines as discussed here.

A growing repertoire of tools now annotate and quantify changes to splicing events. Quantification of features such as splice sites, and exon/intron junctions found in alignment files are commonly used to annotate splicing events. Although the true repertoire of splicing events is difficult to capture, conventional processes can be categorised into distinct groups. The most common events are exon skipping, retained introns, mutually exclusive exons, alternative 5′ and 3′ splice sites. More complex regulatory events involve genomic features beyond exons and introns, such as alternative promoter and polyadenylation sites, which result in varying mRNA 5′ and 3′ UTR ends. However, these events are seldom included in most bioinformatics analyses, tools such as CAGER (Cap Analysis of Gene Expression) and DaPars (Dynamic Analysis of Alternative PolyAdenylation from RNA-Seq) are available for niche research.^23,24 Visualisation of AS is predicated upon the level of detail required in the analysis. If a highly detailed analysis of individual gene structure is needed, splice graphs, sashimi plots and junction maps are commonly used.^25,26 To visualize changes to groups of transcripts, typically MA and Volcano plots are used much the same way as in differential expression level analysis.²¹

Current statistical methods for differential splicing

Commonly, researchers are interested in comparisons of two or more groups of samples known as differential analyses. Differential gene/transcript expression (DGE/DTE) of genes or transcripts involves taking raw read count data, normalizing or scaling it, and calculating whether the changes in expression levels between different biological groups are statistically significant. Differential transcript/exon usage (DTU/DEU), however, uses gene-level group modelling to assess whether the proportional use of the feature (exon or transcript) is statistically significant. Differential splicing events (DSE) on the other hand use a diverse array of statistical methods to quantify and infer splicing events. A comprehensive summary of differential splicing tools is described in the supplementary table (Supplementary Table 1) and in the following sections.

Parametric & mixed methods

Differential expression analysis tools began in the early 2000s coinciding with the development of high throughput technologies such as microarrays. An early example was LIMMA (Linear Models for Microarray Data), developed by Gordon Smyth and colleagues in 2003, which utilises a linear regression framework and empirical Bayes techniques to identify differentially expressed genes.²⁷ Whilst initially only utilised for microarrays, the functionality thus extended to RNASeq data and has been one of the most cited RNASeq methods. As the field shifted from microarray technology to RNASeq, methods were developed such as DESeq (Differential Expression Analysis for Sequence Count Data) and edgeR to capture the nature of count data better and improve modelling.^21,27,28 A major change incorporated in DESeq2 was empirical Bayes-based shrinkage to improve gene-wise variance estimation enhancing accuracy ( Figure 1). Secondly, GLMs (Generalized Linear Models) replaced the simple linear models as these were shown to adapt well to non-normally distributed count-based data.²¹ The flexibility of GLMs allowed algorithms to effectively deal with issues such as overdispersion, shrinkage, heteroscedasticity and covariates. To date, GLMs are usually fitted to the NB (Negative Binomial) distribution which confers some strong advantages. The NB distribution effectively captures overdispersion (the empirical variability of counts) and can handle a large excess of zero values commonly seen in transcript or exon-level count data. However, limma, DESeq2 and edgeR were not developed to specifically address the challenges of identifying AS.

Figure 1. Timeline of statistical methods in differential splicing tool development.

Methods are categorized into parametric and non-parametric approaches, grouped by methodological families. The classification is based on the underlying statistical procedures used for modelling or hypothesis testing, as detailed in Supplementary Table 1. Note that some methods incorporate elements of both parametric and non-parametric frameworks, resulting in overlapping features.

In 2014, DEXSeq was introduced by Michael Love and colleagues, a framework based on DESeq2’s GLM NB model becoming the de-facto tool for parametric splicing-based analysis. Instead of analysing gene-level differential expression, DEXSeq identifies exons within genes that exhibit significant changes in their usage across conditions. This is particularly useful for studying the exonic composition of alternatively spliced transcripts. The development of tools such as DSGseq, rDiff-parametric, JunctionSeq and SeqGSEA has expanded the functionality of the GLM NB family of differential splicing tools.^29–32 DSGseq utilises a holistic approach considering splicing events not as individual elements but as comprehensive gene-wise splice graphs that more accurately reflect complex splicing dependencies.²⁹ The tool rDiff-parametric on the other hand utilises isoform-specific loci such as restricted exonic regions to identify significant differences in isoform composition.³⁰ The proposed advantage of this approach is in the smaller exonic regions rather than full isoform deconvolution. Assigning reads to isoforms is challenging because these transcripts are practically identical, making it difficult to definitively attribute a read from an overlapping region to a particular region without supplementary data. Therefore, full isoform deconvolution is significantly biased against genes with many isoform variants.³³

A few newer methods such as DRIMSeq and DTUrtle have progressed onto non-parametric or mixed Dirichlet Multinomial Models (DMM) which have been argued to capture better the complex variability of count data and better estimate isoform abundance^34,35 ( Figure 1 & Supplementary Table 1). Other methods such as IsoformSwitchAnalyzeR and some custom DEXSeq workflows now incorporate modularity allowing users a selection of bioinformatics tools for filtering, hypothesis testing and posterior calculations.^36,37 An example of the usage of parametric analysis was in the discovery of a chimeric fusion transcript of PRKACA and DNAJB1 in a rare liver tumour FL-HCC (fibrolamellar hepatocellular carcinoma) using DEXSeq’s differential exon usage framework.³⁸ The discovery of differential exon usage of PRKACA’s exons 2-10 and subsequent decreased usage of DNAJB1’s exons 2-3 led the researchers to identify a chimeric transcript in FL-HCC patients. This demonstrated the utility of smaller exon-based analysis in identifying differences in transcript structure which would not be detected in larger gene or transcript-based analysis alone.

Popularity & developer maintenance of methods

To assess the academic popularity of tools, a citation and developer engagement analysis of original research articles within the Web of Science (WoS) domain and the respective GitHub website domains (if applicable). The assessment spanned from 2010 to 2024 and encompassed 19 original papers on various differential splicing tools. Notably, the citation counts for these splicing tools were considerably lower compared to conventional RNA-Seq differential expression analysis tools. For instance, while the general purpose DGE/DTE tool DESeq2 amassed a total of 35,887 citations during the same period, citations for differential splicing tools ranged from 7 to 1300 (Figure 2). This discrepancy may pose challenges for researchers seeking resources and workflows specific to differential splicing analysis. Additionally, the importance of developer support cannot be understated, as it directly influences the usability and longevity of software tools. Notably, differential splicing tools such as DEXSeq, EBSeq, rMATS, SUPPA2, and MAJIQ^{26,39,40,47,48} have shown increasing usage and ongoing developer engagement, as evidenced by their growing citation counts and sustained support ( Figure 3; Figure 4). One possible explanation for the lower citation rates observed in exon/transcript-based methodologies could be the broader adoption of general-purpose differential expression workflows, like DESeq2 that can employ DTE.²¹ Researchers may prefer more explicit splicing event-based tools for targeted splicing analyses and defer to DTE for transcript-based analyses. While the nuances between DTU and DTE may not be a primary focus for many researchers, it is a distinction worth noting in the context of differential splicing analysis.

Figure 2. Citation counts of differential splicing tools (2010–2024) from Web of Science (WoS) Data.

Total citation counts for surveyed differential splicing tools (2010–2024) from the Web of Science Data Portal (WoS). Tools are categorized by analysis level: event, exon, or transcript. DRIMSeq’s original paper was excluded from the citation frequency analysis as it was not indexed in WoS. Certain data included herein are derived from Clarivate Web of Science. © Copyright Clarivate 2023. All rights reserved.Total citation counts for surveyed differential splicing tools (2010–2024) from the Web of Science Data Portal (WoS). Tools are categorized by analysis level: event, exon, or transcript. DRIMSeq’s original paper was excluded from the citation frequency analysis as it was not indexed in WoS. Certain data included herein are derived from Clarivate Web of Science. © Copyright Clarivate 2023. All rights reserved.

Figure 3. Citation trends of differential splicing tools (2010–2024) from Web of Science (WoS) Data.

Annual citation frequency for current differential splicing tools (2010–2024) from Web of Science (WoS). Tools are categorized by analysis level: event, exon, or transcript. DRIMSeq’s original paper is excluded as it is not indexed in WoS. Certain data included herein are derived from Clarivate Web of Science. © Copyright Clarivate 2023. All rights reserved.

Figure 4. Developer maintenance of differential splicing tools.

Annual GitHub repository commits (2010–2024) by category, highlighting community-led maintenance of differential splicing tools. Tools without GitHub pages (MAJIQ, MISO, DSGseq, and dSpliceType) were excluded from the analysis.

The decision between exon/transcript-level (typically parametric) and event-level (typically non-parametric) analyses hinges on several factors, including the particular scientific inquiry, data accessibility, and the level of granularity required to address the research goal. In certain scenarios, integrating both methodologies could offer a more holistic understanding of splicing control mechanisms and their biological significance.

Benchmarking of methods is difficult

To evaluate the quality of differential splicing bioinformatics tools, several benchmarks have been conducted to date. Benchmarking either the scientific accuracy or the computational power of methods can be challenging due to several factors. The main issue is the lack of ground truth to set as a reference to compare measurements to. Commonly, a small subset of experimentally validated splicing events is used as a gold standard to compare against. This was demonstrated in a recent systematic evaluation of 10 differential splicing tools in 2019, where a total of 62 qPCR-validated differentially spliced genes were tested.⁴⁹ The results from this benchmark revealed weak consensus over tool quality as the performance was markedly different across the 4 human and mouse cancer datasets. This demonstrates another issue with these evaluations: inherent heterogeneity in RNASeq data. Often, the performance of methods will depend on the upstream RNASeq pre-processing steps such as in library size, sequence depth, positional bias and annotation quality. To mitigate these issues, some papers use simulated data to i) increase the number of differentially spliced genes to reference and ii) achieve finer control over ground truth and variability within the data.^50–52 One such benchmark used RSEM-based simulated data based on a human prostate cancer dataset (GSE22260⁵³).⁵¹ Another comparison utilised a combination of experimental and simulated Arabidopsis heat shock RNASeq datasets using the Flux Simulator tool.⁵⁴ However, it is important to note that simulated data lacks the complexity of typical biological data. Confounding factors such as outliers, and technical/procedural biases cannot be modelled in current simulations.

The consensus drawn from these three benchmarks is that the performance of differential splicing tools exhibits considerable variability depending on the outlined factors. The ongoing evolution and upkeep of tools by developers introduce a time-dependent aspect to benchmarking. Community-led maintenance efforts consistently enhance the functionality and reliability of tools over time. Rather than aiming for a singular optimal tool for differential splicing analysis, researchers should contemplate employing a suite of tools tailored to address specific inquiries.

Method recommendations

A diagram outlining optimal tool selection is provided to guide prospective alternative splicing (AS) researchers ( Figure 5). Initially, researchers should evaluate the scope and objectives of their analysis. For instance, if the aim is to identify known transcripts, it is advisable to opt for a parametric transcript-based tool like DEXSeq or DRIMSeq and execute a DTU study following Michael Love’s protocol.³⁷ Nonetheless, variations in experimental parameters such as sample size or covariate inclusion may necessitate alternative approaches.

Figure 5. Guideline for differential splicing tool selection based on experimental parameters.

Decision tree for differential splicing analysis, categorized by three branches based on the level of analysis. Transcript-based methods are represented in blue, exon-based methods in pink, and event-based methods in yellow.

If the objective is to uncover novel transcripts, an exon-based parametric approach might be better suited. This choice circumvents the challenges associated with isoform deconvolution and the breadth of transcript annotation, given the smaller exonic regions. For general-purpose differential exon usage (DEU) analysis, DEXSeq remains the preferred protocol due to its robust and flexible statistical methods, as well as its actively maintained software.²¹ However, again intricacies within the data may prompt the usage of more specialised alternatives. Transcript and exon-based methods offer top-down visualizations such as MA/Volcano plots, heatmaps and proportional transcript/exon graphs. If the analysis aims to visualise the movement of exons/introns and splice sites, then an event-based protocol would be more appropriate. Generally, tools such as rMATs, SUPPA2 and MISO offer comprehensive and detailed splicing event analysis.^39,40,55

Commonly, sashimi plots are the best method to visualise splice junctions from aligned data with events annotated, although this can also be plotted separately in IGV.⁵⁶ For user-friendly visualization, MAJIQ offers a summative HTML-based visualizer for complex events such as exitrons or orphan junctions.²⁶ Another factor to consider is the annotation quality of the organism/tissue being studied. If researchers are not confident in the quality of annotations and would like annotation-free analysis, methods such as LeafCutter are a good alternative to conventional methods.^57,58 Overall, event-based methods are more suited to advanced programmers owing to their use of command-line tools over interpreters that use IDEs (Integrated development environments). For most analyses, however, a DEU or DTU-based analysis is recommended for simple interpretability and robustness. Optional steps for AS-specific analyses can also be performed to enhance the data quality. For example, Portcullis enables the accurate filtering of false splice junctions that are often incorrectly characterized by common aligners.⁵⁹

Discussion

While the repertoire of tools to accommodate differential splicing analysis has grown in the past two decades, they are ultimately limited by the capabilities of the RNASeq technology available to date. Since 2010 however, the development of nanopore sequencing technology such as Oxford Nanopore Technologies (ONT) and PacBio’s single-molecule real-time (SMRT) has facilitated the development of long-read RNAseq.^60–62 Long read lengths typically fall within the range of 10kb to 100kb, with ultra-long read lengths now up to 1-2 Mb.⁶³ The main benefit this technology confers is the ability to bypass the aforementioned deconvolution issue stemming from multiple mapping and reconstruct full-length transcript isoforms in a single read. This can not only more accurately identify known transcripts but also novel or splice variants as well as fusion genes. Most current parametric DS tools can therefore be utilised in long-read-based analyses. A recent study utilized IsoformSwitchAnalyzeR’s DEXSeq-based DTU workflow on ONT long reads, demonstrating the capability of current long-standing methods on long-read data.^21,36,64 This was facilitated through long-read custom annotation of the transcriptome using TALON to identify novel transcripts.⁶⁵ Additionally, specific novel technologies such as LIQA have been developed to analyse long reads.⁶⁶

Long-read RNAseq still possess notable disadvantages, however. Early on, long-read RNASeq possessed error rates of 10-20%.^67,68 The development of HiFi sequencing by PacBio using circular consensus sequencing has since reduced the error rate to a reported 0.5%.⁶⁹ While the development of deep-learning algorithms such as DeepConsensus has sought to push HiFi accuracy further bringing it on par to short read.⁷⁰ However, this is still highly dependent on the depth of sequencing. The most efficient error correction method involves hybridising the analysis with short-read RNAseq methods.⁷¹ This ultimately means that while accuracy can now be brought to close to 99.5%, error correction drives the cost of long-read RNAseq methods up significantly. The field is progressing towards optimal error correction and is now focusing on lowering costs which is currently the largest hurdle for practical use for common research.

As interest in alternative splicing grows, researchers have access to an expanding array of tools. Advances in statistical methods and longer RNA sequencing read lengths are overcoming technical limitations. This leads to more precise transcript alignment and reduces the need for complex computational steps. With workflows becoming streamlined and modular, platforms like Nextflow enable researchers to create tailored pipelines for their specific goals and data types.⁷² These developments promise a brighter future for alternative splicing analysis, facilitating a deeper exploration of transcriptomic regulation and its functional significance.

Ethical approval and consent statement

Ethical approval and consent were not required.