DeepSplice: A Transformer-Based Framework for Predicting Alternative Splicing Events from RNA-seq Data

Abstract

Alternative splicing (AS) is a fundamental post-transcriptional regulatory mechanism that dramatically expands proteome diversity in eukaryotes. Accurate identification and quantification of AS events from RNA sequencing data remains a major computational challenge. Here we present DeepSplice, a transformer-based deep learning framework that integrates raw RNA-seq read signals, splice-site sequence context, and evolutionary conservation scores to predict five canonical types of alternative splicing events. Benchmarked on three independent human cell-line datasets, DeepSplice achieves an average AUROC of 0.947 and outperforms state-of-the-art tools by 4-11% on F1 score.

1. Introduction

Alternative splicing enables a single gene to produce multiple mRNA isoforms by varying the selection of exons and introns during pre-mRNA processing. More than 95% of human multi-exon genes undergo alternative splicing, and dysregulation of this process is implicated in a wide spectrum of diseases, including cancer, neurodegeneration, and cardiovascular disorders [1, 2].

Current computational approaches for detecting AS events from RNA-seq data can be broadly divided into three categories:

1. Alignment-based methods (e.g., rMATS [3], DEXSeq [4]) that rely on read counts at annotated splice junctions.
2. Assembly-based methods (e.g., StringTie [5], Trinity [6]) that reconstruct full-length transcripts before quantification.
3. Machine learning methods (e.g., SplAdder [7], VAST-TOOLS [8]) that model splicing as a classification or regression problem.

Despite significant progress, existing tools still suffer from limited sensitivity for low-coverage events, high false-positive rates in repetitive genomic regions, and poor generalization across tissue types and species. Deep learning, particularly transformer architectures [9], has recently demonstrated superior capacity for capturing long-range dependencies in biological sequences [10, 11], motivating us to develop DeepSplice.

In this work we make the following contributions:

• A novel multi-modal transformer architecture that jointly encodes RNA-seq coverage profiles, primary splice-site sequences, and PhyloP conservation scores.
• A hierarchical attention mechanism that identifies the most informative read-level and nucleotide-level features for each AS event type.
• Comprehensive benchmarks on six public datasets spanning three human cell lines and two mouse tissues.
• A downstream application demonstrating clinically relevant splicing disruptions across 23 TCGA cancer cohorts.

2. Methods

2.1 Data Collection and Preprocessing

We obtained paired-end RNA-seq data (2x150 bp, 50M+ read pairs) for three human cell lines from ENCODE:

Reads were aligned to GRCh38 (GENCODE v43) using STAR 2.7.10a with default two-pass mode parameters. rMATS 4.1.2 was used to generate a gold-standard set of AS events with inclusion level difference |DELTA-PSI| > 0.1 and FDR < 0.05. This produced 187,432 high-confidence AS events distributed as follows:

• SE (exon skipping): 98,741 (52.7%)
• RI (intron retention): 34,218 (18.3%)
• A5SS: 22,651 (12.1%)
• A3SS: 21,934 (11.7%)
• MXE: 9,888 (5.3%)

Negative examples (constitutively spliced junctions) were sampled at a 2:1 ratio to positive events, stratified by gene expression level to avoid confounding.

2.2 Feature Engineering

For each candidate AS event, we extracted three complementary feature modalities:

2.2.1 Coverage Profile Tensor

RNA-seq read coverage was computed over a 400-nt window centered on each splice site using samtools 1.17. Coverage values were normalized per million mapped reads (RPM) and log-transformed: c_hat = log2(c + 1). The resulting 1D signal was discretized into 20-nt bins, producing a 20-dimensional vector per splice site, and both the upstream and downstream splice sites were concatenated to form a 40-dimensional coverage profile.

2.2.2 Sequence Context Embedding

The +/-200 nt genomic sequence flanking each splice site was one-hot encoded (4 channels x 400 positions). Additionally, six canonical splice-site features were extracted: GT-AG, GC-AG, AT-AC dinucleotides, branch point score (computed with SVM-BPfinder [12]), polypyrimidine tract length, and MaxEntScan [13] 5'SS / 3'SS scores.

2.2.3 Evolutionary Conservation

Per-nucleotide PhyloP 100-way vertebrate conservation scores were downloaded from UCSC and averaged across the same 400-nt windows to generate a 20-dimensional conservation vector.

2.3 DeepSplice Architecture

DeepSplice employs a three-branch encoder followed by a cross-modal transformer fusion module.

Input Modalities
      |
+-----+---------------------+
|     |                     |
v     v                     v
Coverage  Sequence Context  Conservation
1D-CNN    BERT-style        MLP
(3 layers) Transformer     (2 layers)
|          (6 heads,d=256)  |
+-----------------------------+
               |
      Cross-Modal Attention
        (4 heads, d=512)
               |
       Classification Head
    (5 binary output neurons)

Coverage encoder: Three 1D convolutional layers (kernel sizes 3, 5, 7; 64 filters each) with batch normalization and ReLU activations, followed by global average pooling.

Sequence encoder: A 6-layer BERT-style transformer (hidden size 256, 8 attention heads, feed-forward size 1024) pre-trained on 50M human intron/exon sequences using masked nucleotide prediction.

Conservation encoder: A two-layer MLP (256 -> 128 -> 64 units) with dropout (p=0.3).

Fusion: The three branch representations are projected to a common 512-dimensional space and fused via cross-modal multi-head attention followed by a two-layer classification MLP with sigmoid output.

The loss function is a weighted binary cross-entropy to account for class imbalance:

$$\mathcal{L} = -\frac{1}{N}\sum_{i=1}^{N} \left[ w_+ y_i \log \hat{y}$$i + w- (1-y_i) \log (1-\hat{y}_i) \right]

where $w_+ = N / (2 N_+)$ and $w_- = N / (2 N_-)$ are class weights inversely proportional to class frequencies.

2.4 Training Details

Models were trained using AdamW (lr=3e-4, weight decay=0.01) with a cosine annealing schedule over 50 epochs. Batch size was 256. Early stopping with patience=10 was applied on the validation AUROC. All experiments used 5-fold cross-validation with chromosome-level splits to prevent data leakage. Training was performed on 4x NVIDIA A100 (80 GB) GPUs using PyTorch 2.1 with mixed-precision (FP16) training.

2.5 Baseline Methods

We compared DeepSplice against five published tools:

• rMATS 4.1.2: statistical model based on read counts at annotated junctions
• SUPPA2 2.3.3: likelihood-ratio framework leveraging transcript quantification
• SplAdder 3.0.0: graph-based augmented splice graph approach
• Whippet 1.6: lightweight probabilistic model
• DARTS 0.1: deep-learning model using only sequence features

3. Results

3.1 Overall Performance

DeepSplice achieves state-of-the-art performance across all five AS event types and all three cell lines. Table 1 summarizes the average metrics over 5-fold cross-validation.

Table 1. Performance comparison (mean +/- SD across 5 folds, GM12878+HepG2+K562 combined)

DeepSplice improves F1 score by 10.6% over rMATS, 11.3% over Whippet, and 5.3% over DARTS, demonstrating substantial gains across all comparison methods.

3.2 Per-Event-Type Analysis

Performance varies by event type, with exon skipping (SE) being the easiest to predict (AUROC=0.962) and mutually exclusive exons (MXE) the most challenging (AUROC=0.921):

3.3 Ablation Study

To quantify the contribution of each input modality, we trained DeepSplice variants with individual modalities removed.

Table 2. Ablation study -- AUROC on combined test set

Sequence context provides the largest individual contribution, followed by coverage profiles, consistent with the known importance of splice-site consensus sequences. Cross-modal attention fusion outperforms simple concatenation by 1.8%, validating the design choice.

3.4 Generalization to Mouse Tissues

To evaluate cross-species generalizability, we applied DeepSplice (trained on human data without retraining) to mouse hippocampus and liver RNA-seq data from ENCODE. The model achieves AUROC=0.921 (hippocampus) and AUROC=0.934 (liver), indicating strong zero-shot generalization to mouse splicing patterns.

3.5 Cancer Splicing Landscape

We applied DeepSplice to 9,328 tumor samples across 23 TCGA cancer types. DeepSplice identified 12,847 recurrent tumor-specific AS events (present in >= 10% of samples in >= 2 cancer types) not detected by rMATS. Pathway analysis revealed significant enrichment (adjusted p < 0.001) of splicing alterations in apoptosis regulators (BCL2L1, CASP9), cell cycle checkpoints (BRCA1 exon 11, MDM2 exon 3), and chromatin remodeling (BRD4, KDM6A).

Notably, DeepSplice recovers known oncogenic splice variants:

• EGFRvIII exon 2-7 skipping in glioblastoma (predicted PSI=0.41 vs. RT-PCR measured PSI=0.39, Pearson r=0.97)
• MET exon 14 skipping in lung adenocarcinoma (F1=0.93)
• CD44 variable exon switching in breast cancer (AUROC=0.944)

3.6 ESE Mutation Impact Prediction

We evaluated DeepSplice on 3,219 clinically annotated exonic splicing enhancer (ESE) variants from the Human Splicing Finder database. DeepSplice correctly predicts splicing disruption for 83.4% of pathogenic ESE mutations vs. 71.2% for SpliceAI [14] and 67.8% for MaxEntScan, suggesting that the multi-modal architecture provides complementary information beyond sequence context alone.

4. Discussion

DeepSplice demonstrates that integrating multiple evidence streams through a cross-modal transformer architecture leads to substantial improvements in AS event detection. The pre-trained sequence encoder likely captures complex splice regulatory elements (SREs) beyond the canonical GT-AG dinucleotides, including distant branch points and exonic splicing silencers.

Several limitations should be noted. First, DeepSplice requires at least 20x coverage depth for reliable predictions; performance degrades for low-input or single-cell RNA-seq data. Second, the current model does not explicitly model tissue-specific splicing factor expression, which could be incorporated as an additional conditioning signal in future work. Third, while cross-species transfer to mouse works reasonably well, performance in more distant organisms warrants further investigation.

Future directions include: (1) extending to long-read RNA-seq (PacBio/Oxford Nanopore) to resolve complex isoform structures; (2) incorporating protein structural context for functional impact scoring; (3) developing a fine-tuning protocol for clinical samples with limited data.

5. Conclusion

We present DeepSplice, a multi-modal transformer framework for predicting alternative splicing events from RNA-seq data. DeepSplice achieves state-of-the-art performance on standard benchmarks, generalizes across species, and identifies clinically relevant splicing alterations in cancer. The model and training code are released under MIT License.

References

[1] Pan, Q. et al. Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing. Nat. Genet. 40, 1413-1415 (2008).

[2] Wang, E. T. et al. Alternative isoform regulation in human tissue transcriptomes. Nature 456, 470-476 (2008).

[3] Shen, S. et al. rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-seq data. Proc. Natl. Acad. Sci. USA 111, E5593-E5601 (2014).

[4] Anders, S. et al. Detecting differential usage of exons from RNA-seq data. Genome Res. 22, 2008-2017 (2012).

[5] Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290-295 (2015).

[6] Grabherr, M. G. et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat. Biotechnol. 29, 644-652 (2011).

[7] Kahles, A. et al. SplAdder: identification, quantification and testing for differential splicing using RNA-seq data. Bioinformatics 32, 1840-1847 (2016).

[8] Tapial, J. et al. An atlas of alternative splicing profiles and functional associations reveals new regulatory programs and genes that simultaneously express multiple major isoforms. Genome Res. 27, 1759-1768 (2017).

[9] Vaswani, A. et al. Attention is all you need. Adv. Neural Inf. Process. Syst. 30 (2017).

[10] Avsec, Z. et al. Effective gene expression prediction from sequence by integrating long-range interactions. Nat. Methods 18, 1196-1203 (2021).

[11] Chen, K. M. et al. A sequence-based global map of regulatory activity for deciphering human genetics. Nat. Genet. 54, 940-949 (2022).

[12] Corvelo, A. et al. Genome-wide association between branch point properties and alternative splicing. PLoS Comput. Biol. 6, e1001016 (2010).

[13] Yeo, G. & Burge, C. B. Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals. J. Comput. Biol. 11, 377-394 (2004).

[14] Jaganathan, K. et al. Predicting splicing from primary sequence with deep learning. Cell 176, 535-548 (2019).