1. Introduction

Isoforms represent distinct transcript variants generated from a single gene through alternative splicing (AS), serving as the primary executors of cellular functional diversity^[1,2]. A single gene can produce multiple protein-coding isoforms that carry distinct or even antagonistic biological functions, which are essential for understanding complex disease mechanisms and identifying therapeutic targets^[3]. For example, the TP53 gene produces alternatively spliced isoforms, such as p53 and p47, that regulate distinct stages of the cell cycle; specifically, p53 controls the G1 transition, whereas p47 induces G2 arrest^[4]. These cases highlight that functional analysis at the isoform level provides higher resolution and clinical value for understanding disease mechanisms compared to gene-level studies^[5]. Despite their significance, determining the specific functional roles of isoforms remains a formidable challenge due to both experimental and computational constraints. While advanced techniques like single-cell long-read sequencing can identify full-length transcripts at single-cell resolution^[6-8], the widespread application of these technologies is frequently hindered by low capture efficiency, high sequencing error rates, and the necessity for complex computational deconvolution^[9-12]. Furthermore, large-scale wet-lab validation of tens of thousands of isoforms is labor-intensive and costly^[13], making high-precision computational approaches imperative for functional characterization.

Beyond transcript-level quantification, the granular analysis of the specific splicing mechanisms that generate diverse isoforms has garnered significant research interest^[14]. Recent studies have demonstrated that splicing event patterns^[15], characterized by the Percent Spliced In metric, show significant specificity across cell types and developmental stages. These features have been successfully applied to downstream tasks such as cell-type clustering and cell-state transition analysis, often providing superior discriminative power compared to traditional gene expression profiles^[16,17]. Since the exon alterations caused by these events directly determine the domain composition and functional attributes of the resulting protein products^[18], the association between specific splicing events and Gene Ontology (GO) terms serves as a crucial entry point for understanding functional differentiation.

Existing computational methods for isoform function prediction can be broadly categorized into three groups, which include early multi-instance learning (MIL) methods, network and matrix factorization, and deep learning-based approaches. To address the severe scarcity of isoform-level labels, early approaches rely on the traditional MIL paradigm. The core idea is to formulate genes as bags and their alternatively spliced isoforms as instances, allowing gene-level GO annotations to serve as bag-level supervisory signals. Representative methods such as multiple-instance (MI) support vector machine (SVM)^[19], instance-oriented MI label propagation^[20], and weighted logistic regression-based MIL method^[21] employ SVMs, label propagation, and sparse regularization within this framework to iteratively assign functional labels to individual isoforms. However, these early MIL strategies typically assume that each specific function of a gene is independently driven by a single isoform, failing to account for the collaborative actions of multiple isoforms and the complex hierarchical structure of GO terms.

To incorporate biological topologies and alleviate the reliance on complete annotations, subsequent studies shifted toward heterogeneous networks and matrix factorization. These approaches project various biological entities into a unified graph or latent space to propagate functional signals across different omics layers. Notable algorithms include IsoFun^[22], DisoFun^[23], and IsoDA^[24], which utilize bi-random walks, collaborative matrix factorization, and domain adaptation to infer missing annotations based on co-expression and sequence similarities. While successfully capturing the global biological context, these methods rely on shallow and linear architectures, which restricts their capacity to model the differential importance of various biological associations or to uncover the deep, non-linear interactions between protein products and complex GO hierarchies.

Recently, deep learning and graph neural networks (GNNs) have been widely deployed in this field, offering powerful non-linear expressiveness to automatically extract features from multi-modal data. Methods like deep learning-based prediction of IsoForm FUnctions from sequences and expression^[25], deep multi-instance learning (DMIL)-IsoFun^[26], and FINER^[27] employ deep convolutional architectures to integrate sequence and expression profiles to guide more precise functional assignments. Moreover, graph-aware frameworks have been developed to model hierarchical structures and collaborative relationships among isoforms; for instance, GraphIsoFun^[28] and IsofunGO^[29] utilize topological embeddings, while CrossIsoFun^[30] further introduces cross-dataset transfer mechanisms to enhance generalization. Despite these advancements, existing deep models still face two critical bottlenecks. First, by ignoring the functional roles of local splicing events, holistic encoding of highly homologous isoforms yields nearly identical global embeddings, rendering models unable to distinguish their specific functions. Second, these methods often struggle with the extreme sparsity and severe class imbalance of isoform-level GO annotations, frequently leading to overfitting or prediction failure on highly specific, rare functional terms.

To address these challenges, we propose SpliceEM, a novel framework for high-precision isoform function prediction. SpliceEM addresses these limitations through three main components. First, it introduces an event-aware heterogeneous graph transformer (HGT) that explicitly integrates AS events into a unified latent space. Utilizing cross-modal attention, it disentangles fine-grained functional signals from highly homologous global sequences. Second, to tackle extreme class imbalance and complex hierarchical logic, SpliceEM constructs a joint optimization module featuring a gene-masked asymmetric loss (ASL)^[31] and hierarchical topological constraints. This module effectively filters gradient interference from the vast majority of unrelated background annotations while enforcing strict adherence to the biological ontology. Furthermore, to overcome the severe scarcity of ground-truth annotations, SpliceEM develops an Exponential Moving Average (EMA) teacher-guided self-knowledge distillation framework. By integrating MIL logic, this architecture synergistically leverages abundant macro-biological gene priors and limited precise isoform labels to safely explore the unannotated space.

Extensive experiments on a curated human dataset demonstrate that SpliceEM significantly outperforms state-of-the-art methods, including IsofunGO and CrossIsoFun. In particular, SpliceEM successfully characterizes cell-type-specific functional profiles and identifies key isoforms responsible for specific biological processes (BPs) within distinct cellular environments. By shifting the focus from global sequences to localized splicing events and utilizing hierarchical biological priors, SpliceEM provides a robust and interpretable algorithmic foundation for decoding transcriptomic functional diversity. These findings confirm that SpliceEM serves as a valuable tool for understanding how cellular contexts shape isoform functional activities, providing a baseline for deeper investigations into regulatory mechanisms and disease development.

2. Methods

As illustrated in Figure 1, the SpliceEM framework predicts isoform-level functions by integrating three core modules: 1) a heterogeneous graph construction module, which integrates multi-modal biological features to establish a unified topological network; 2) an HGT-based message passing module, which captures deep cross-modal interactions to extract highly discriminative isoform embeddings; and 3) an EMA-MIL joint optimization module, which utilizes teacher-student distillation to fuse broad macro-level gene priors with sparse but highly precise micro-level isoform annotations. Notably, within this optimization module, a gene-masked ASL explicitly leverages these precise isoform-level ground truths to combat extreme class imbalance, while hierarchical triplet constraints are applied to preserve the GO topology. Through these three interconnected stages, SpliceEM effectively translates heterogeneous biological features into accurate functional probability matrices.

Display Full Size

Download

Figure 1. Overall architecture of the SpliceEM framework. (A) Heterogeneous Graph Construction: Integrates multi-modal features (sequences, GO terms, and AS events) into a unified topological network; (B) HGT-based Message Passing: Captures cross-modal interactions via a Heterogeneous Graph Transformer to extract discriminative isoform embeddings; (C) EMA-MIL Joint Optimization: Employs a teacher-student framework with gene-masked ASL and hierarchical constraints to achieve robust isoform function prediction. GO: gene ontology; AS: alternative splicing; HGT: heterogeneous graph transformer; EMA: exponential moving average; MIL: multi-instance learning; ASL: asymmetric loss; BERT: bidirectional encoder representations from transformers; ESM: evolutionary scale modeling; BCE: binary cross-entropy.

2.1 Heterogeneous graph construction

To comprehensively capture the multidimensional interactions within the transcriptome, we model the biological system as a global heterogeneous graph G = (V,E) comprising four basic node types: isoforms, genes, GO terms, and AS events.

Graph Topology Construction: The graph edges (E) are rigorously established based on multi-source biological databases. Specifically, the gene-isoform subordinate relationships are constructed based on the Ensembl genome annotations; the isoform-AS event connections are identified and mapped using SUPPA2^[32]; and the directed hierarchical dependencies (GO-GO), such as is_a and part_of relationships, are constructed according to the directed acyclic graph (DAG) provided by the Gene Ontology Consortium^[33].

Feature Initialization and Mapping: Having established the graph topology, we extract specific raw biological features for each node type to assign them meaningful representations. The specific definitions and mapping processes for these features are detailed as follows.

For isoform and gene nodes, protein sequences for valid isoforms are retrieved from the UniProt database via Ensembl ID mapping. Global sequence representations are then extracted using the pre-trained protein language model ESM-2^[34]. Host gene features are aggregated from the representations of all known isoforms they encode. Both are independently projected into the hidden space through specific non-linear transformation modules:

(1)$h_{\text{iso}}^{\left(0\right)}=\mathrm{LayerNorm}\left(\mathrm{ReLU}\left(W_{\text{iso}}{x}_{\text{esm}}\right)\right),h_{\text{gene}}^{\left(0\right)}=\mathrm{LayerNorm}\left(\mathrm{ReLU}\left(W_{\text{gene}}{x}_{\text{gene}}\right)\right)$

where $x_{\mathrm{esm}}\in R^{d_{\mathrm{esm}}}(d_{\mathrm{esm}}=1,280)$ is the raw isoform sequence feature vector extracted by ESM-2; $x_{\text{gene}}\in R^{d_g}$ is the aggregated raw gene feature vector; $W_{\text{iso}}\in R^{d_h\times d_{\mathrm{esm}}}$ and $W_{\text{gene}}\in R^{d_h\times d_g}$ are learnable projection matrices specific to isoform and gene nodes, respectively; ${\text{h}}_{\text{iso}}^{\left(0\right)},{\text{h}}_{\text{gene}}^{\left(0\right)}\in {\text{R}}^{{\text{d}}_{\text{h}}}$ are the unified representations at the initial layer (Layer 0). ReLU(·) and LayerNorm(·) denote the Rectified Linear Unit activation and Layer Normalization, enhancing non-linear expressiveness and numerical stability.

To represent AS events, we utilize the seven basic types of AS events extracted by SUPPA2: skipping exon (SE), retained intron (RI), alternative 3’ splice site, alternative 5’ splice site (A5), alternative first exon (AF), alternative last exon (AL), and mutually exclusive exons. These events are discretely encoded as one-hot feature vectors and mapped into the same representation space:

(2)$h_{\text{splice}}^{\left(0\right)}=\mathrm{LayerNorm}\left(\mathrm{ReLU}\left(W_{\text{splice}}{V}_{\text{type}}\right)\right)$

where $v_{\text{type}}\in \{0,1{\}}^{d_{\mathrm{as}}}(d_{\mathrm{as}}=7)$ is a one-hot vector representing the specific AS category, and $W_{\text{splice}}\in R^{d_h\times d_{\mathrm{as}}}$ is the learnable projection matrix for splicing event nodes.

Finally, the raw semantic features of GO terms are initialized using pre-trained language models and projected through a corresponding transformation matrix to obtain $h_{\mathrm{go}}^{\left(0\right)}\in R^{d_h}$.

By applying these specific non-linear transformations, all heterogeneous nodes are successfully aligned within the same d_h-dimensional latent space. This unified alignment ensures that subsequent network layers can process multimodal inputs with varying raw dimensions.

2.2 HGT-based message passing

Due to the inherent heterogeneity of nodes and edges in the graph, standard GNNs are intrinsically limited in modeling connections with diverse semantics. To address this, SpliceEM employs a multi-layer HGT architecture. Since HGT explicitly utilizes distinct parameter matrices for different node and edge types, it can effectively capture deep cross-modal interactions and complex topological dependencies. For the message-passing process at layer 1(l ≥ 1):

(3)$\mathrm{Attention}(s,t)=\mathrm{Softmax}(\frac{\left(W_Q^{\left(t\right)}{h}_t^{(l-1)}\right)W_{\mathrm{att}}^{\left(e\right)}{\left(W_K^{\left(s\right)}{h}_s^{(l-1)}\right)}^{\top }}{\sqrt{d_k}}\cdot {\mu }^{\left(e\right)})$

where s and t denote the source and target nodes; e represents the specific heterogeneous relationship edge type; $h_s^{(l-1)},h_t^{(l-1)}\in R^{d_h}$ are the node feature vectors from the previous layer. Specifically, when l = 1, they correspond to the initialized features (e.g., ${\text{h}}_{\text{iso}}^{\left(0\right)},{\text{h}}_{\text{splice}}^{\left(0\right)}$). $W_Q^{\left(t\right)},W_K^{\left(s\right)}\in R^{d_k\times d_h}$ are the target-specific query and source-specific key projection matrices; d_k is the channel dimension per attention head; $W_{\mathrm{att}}^{\left(e\right)}\in R^{d_k\times d_k}$ is the interaction weight matrix for edge type e; and μ^(e) ∈ R is the learnable attention scaling factor.

Source node information is then mapped via a type-specific value matrix and aggregated to the target node based on the computed attention weights:

(4)$\mathrm{Msg}(s,t)=\mathrm{Attention}(s,t)\cdot W_V^{\left(s\right)}{h}_s^{(l-1)}$

where $W_V^{\left(s\right)}\in R^{d_k\times d_h}$ is the source-specific value projection matrix. Upon completing aggregation from all semantic neighbors, the feature of target node t is updated to $h_t^{\left(l\right)}$. After L layers of message passing, the network maps the initial feature set to the final representations H^(L), which comprehensively fuse global topological and multimodal functional signals.

2.3 EMA-MIL joint optimization

In isoform function prediction, the optimization process must address two critical bottlenecks: the extreme class imbalance of GO terms and the severe sparsity of isoform-level annotations. To resolve this, SpliceEM constructs a joint optimization module consisting of four core components: a self-knowledge distillation framework, a gene-masked ASL, GO hierarchical constraints, and a unified inference strategy.

2.3.1 Self-knowledge distillation framework

To fully leverage macro-biological priors at the gene level and mitigate the extreme sparsity of isoform annotations, we construct a self-knowledge distillation system based on the Mean Teacher paradigm. In this framework, the student and teacher are strictly initialized with the identical HGT architecture and weights, yet they undergo distinct update mechanisms to form a self-reinforcing learning loop.

Specifically, while the student model is actively optimized via gradient descent using the limited ground-truth isoform labels, it inherently struggles to generalize across the vast unannotated space. To safely bridge the gap between abundant gene-level priors and these extremely sparse isoform annotations, the teacher model is strictly detached from backpropagation and serves as a temporal consensus extractor. Its parameters θ_teacher are passively updated via EMA from the student’s trajectory: θ_teacher←αθ_teacher + (1 - α)θ_student. By ensembling the student’s historical states, the teacher avoids extreme fluctuations and generates highly stable predictions to guide the student in uncertain regions.

For candidate samples where the host gene possesses a specific function but the exact isoform attribution remains unknown, we define the candidate activation mask ${\text{M}}_{\text{cand}}=({\text{Y}}_{\text{gene}}^{\text{expand}}==1)\wedge ({\text{Y}}_{\text{iso}}==0)$. Within this space, the teacher model utilizes MIL logic to generate a normalized soft distribution, characterizing the relative potential of different isoforms under the same gene:

(5)${\hat{y}}_{\text{soft}}^{(i,j)}=\frac{\exp \left(z_{\text{teacher}}^{(i,j)}/T\right)}{\sum _{k\in \mathrm{Iso}\left(g\right)}\exp \left(z_{\text{teacher}}^{(k,j)}/T\right)}$

The student model then continuously approximates this stable soft distribution via binary cross-entropy, yielding the distillation loss L_distill. Through this self-distillation mechanism, the model essentially uses its own smoothed historical knowledge to reliably translate broad gene-level priors into specific isoform functions without introducing false-positive noise.

(6)$L_{\text{distill}}=\frac{1}{\left|M_{\text{cand}}\right|}\sum _{(i,j)\in M_{\text{cand}}}\mathrm{BCE}(z_{\text{student}}^{(i,j)},{\hat{y}}_{\text{soft}}^{(i,j)})$

2.3.2 Gene-masked ASL

Protein function prediction inherently suffers from extreme class imbalance, where negative samples outnumber positive ones by orders of magnitude. Standard loss functions treat all errors equally, causing the sparse positive gradients to be entirely overwhelmed by massive well-classified negative samples. To solve this, we introduce an ASL. By decoupling the decay rates for positive and negative classes and introducing a probability margin, the model strictly focuses on learning rare functional signals while gracefully ignoring those correctly classified redundant negatives.

To further mitigate false-negative interference caused by unlabeled isoforms, we combine ASL with gene-level biological masking. If a host gene strictly lacks a function, its subordinate isoforms definitively do not possess it, forming highly reliable deterministic negative samples. The mask M activates only known positive samples and these deterministic negatives: $M=(Y_{\text{iso}}==1)\vee (Y_{\text{gene}}^{\text{expand}}==0)$. Within this space, the ASL is defined as:

(7)$L_{\mathrm{ASL}}=\frac{1}{\left|M_{+}\right|}\sum _{(i,j)\in M_{+}}{L}_{+}^{(i,j)}+\frac{1}{\left|M_{-}\right|}\sum _{(i,j)\in M_{-}}{L}_{-}^{(i,j)}$

(8)$L_{+}=-(1-p{)}^{\gamma +}\log (p),L_{-}=-max(p-m,0{)}^{\gamma -}\log (1-p)$

where p ∈ [0,1] is the predicted probability.

2.3.3 GO hierarchy constraints

Unlike traditional independent flat labels, GO terms are intrinsically correlated through a DAG hierarchy, where specific functions inherit semantics from broader parent terms. To enable the feature space to internalize these topological priors and ensure logically consistent predictions, we introduce a hierarchical triplet margin loss L_hier:

(9)$L_{\text{hier}}=\frac{1}{N_{\mathrm{go}}}\sum _{i=1}^{N_{\mathrm{go}}}max(0,d(e_i,e_j^{+})-d(e_i,e_k^{-})+m_h)$

where $e_j^{+}$ and $e_k^{-}$ are embeddings of a sampled ancestor and non-ancestor node of the anchor term e_i.

2.3.4 Joint optimization and final inference

Integrating the aforementioned asymmetric classification loss, self-knowledge distillation, and hierarchical constraints, the final optimization objective of the SpliceEM framework is defined as:

(10)$L_{\text{total}}=L_{\mathrm{ASL}}+{\lambda }_{\text{dist}}\left(t\right)\cdot L_{\text{distill}}+{\lambda }_{\text{hier}}\cdot L_{\text{hier}}$

where λ_dist(t) is the distillation weight dynamically adjusted via a Gaussian ramp-up strategy (as implemented in our training pipeline) to smoothly introduce soft labels, and λ_hier balances the hierarchical constraints.

During the inference phase, the EMA teacher model, which acts as a robust temporal ensemble of the student’s training trajectory, serves as the final stable inference engine. Based on its generalized isoform embeddings $h_{\text{iso, teacher}}^{\left(L\right)}$, the final functional prediction matrix is calculated as $\hat{Z}=\sigma \left(\mathrm{FC}\left(h_{\text{iso, teacher}}^{\left(L\right)}\right)\right)$. Detailed implementation settings are provided in the Supplementary materials.

3. Results

3.1 Compared methods and experimental setups

To comprehensively evaluate the effectiveness of our proposed framework, we benchmark SpliceEM against five representative state-of-the-art methods for isoform function prediction. Specifically, IsoFun^[22] employs bi-random walks on a heterogeneous network; DisoFun^[23] utilizes collective matrix factorization on expression and gene-GO matrices; GraphIsoFun^[28] applies graph convolutional networks on integrated functional association graphs; IsofunGO^[29] combines hierarchical network embedding with attention-based MIL; and CrossIsoFun^[30] achieves deep cross-omics fusion leveraging generated isoform-level interactome data. For each baseline, we directly utilized the recommended parameter configurations that were provided in their original publications. We modified their publicly available implementations into unified Python code to ensure a fair and consistent evaluation on our dataset.

For our experiments, we curated a human dataset utilizing raw sequencing data obtained from Shiau et al.^[35] (Gene Expression Omnibus accession: GSE212945). We processed the raw reads through HISAT2 and StringTie to identify valid expressed transcript identifiers. These resulting Ensembl IDs were subsequently mapped to the UniProt database to retrieve corresponding protein sequences and ground-truth GO annotations. Additionally, the GO hierarchical structure was acquired from the Gene Ontology Consortium to model the topological dependencies between functional terms, while SUPPA2 was employed to extract AS events. Following this multi-source data integration, our primary dataset encompasses 8,472 genes and 16,100 valid isoforms. To ensure rigorous evaluation and prevent data leakage, all baseline methods were evaluated using the exact same data splits as SpliceEM. Specifically, we partitioned the dataset using a strict gene-centric strategy with an 8:1:1 ratio for training, validation, and testing sets, which ensures that all isoforms encoded by the same gene are partitioned into the same collection. Furthermore, to systematically assess prediction robustness across different levels of data availability, we stratify GO terms into three categories based on their annotation frequencies: sparsely annotated terms ([3, 30)), moderately annotated terms ([30, 100)), and richly annotated terms ([100, 300)). This stratified approach is particularly crucial for assessing robustness in data-sparse scenarios commonly encountered in functional genomics research. We adopt four representative metrics for evaluation: area under the receiver operating characteristic curve (AUROC), area under the precision-recall curve (AUPRC), F_max (maximum F-measure), and S_min (minimum semantic distance). Higher values of AUROC, AUPRC, and F_max indicate superior performance, whereas a lower value is desirable for S_min^[36].

3.2 Results of isoform function prediction

Table 1 evaluates SpliceEM against baselines across the BP, cellular component, and molecular function (MF) categories. Due to the inherent scarcity of highly annotated terms in the MF ontology, data for the [100, 300] range is insufficient (denoted as “—”). It is worth noting that the absolute AUPRC values remain generally low across all evaluated methods. This is an inherent challenge in single-cell isoform function prediction caused by the extreme sparsity of positive labels distributed across a massive GO term space. Nevertheless, under these overwhelmingly harsh conditions, SpliceEM consistently achieves state-of-the-art performance, securing the highest relative gains across all metrics and frequency intervals. This superiority manifests in three key aspects:

First, SpliceEM significantly improves AUPRC and F_max, demonstrating precise discrimination of homologous isoforms. Early baselines like IsoFun and DisoFun perform poorly due to shallow architectures. CrossIsoFun achieves competitive results by performing intermediate fusion of sequence, expression, and inferred isoform interaction graphs via a heterogeneous graph network. However, without explicitly modeling localized splicing variations, highly homologous isoforms remain difficult to fully distinguish in the latent space. SpliceEM overcomes this problem by integrating seven AS events as distinct topological nodes within a unified space. Using cross-modal attention, it effectively extracts splicing-specific functional signals from global sequences, accurately capturing micro-variations that drive functional differentiation.

Furthermore, SpliceEM consistently yields the lowest S_min, indicating stricter adherence to GO hierarchical logic. Traditional methods treat GO terms as flat labels, leading to severe misclassifications. IsofunGO performs competitively by constraining hierarchical relationships using static Lin similarity. Specifically, the triplet loss is utilized to explicitly distinguish highly similar functions in the latent space. Meanwhile, the ASL filters out training interference from the vast majority of clearly unrelated GO terms. Through this combined approach, SpliceEM robustly maintains hierarchical consistency while enhancing fine-grained functional boundaries.

Finally, SpliceEM exhibits remarkable robustness under extreme data sparsity. Under such conditions, existing baselines like CrossIsoFun and IsofunGO suffer from severe performance degradation because they fail to synergistically integrate the abundant macro-biological gene priors with the limited but highly precise isoform-level annotations. SpliceEM elegantly solves this via an EMA teacher-guided distillation framework. Instead of discarding either source of supervision, SpliceEM utilizes MIL to safely fuse broad gene-level priors with the sparse ground-truth isoform labels. Furthermore, establishing the temporally ensembled EMA teacher as the final inference engine effectively mitigates stochastic training fluctuations, ensuring robust exploration of the unannotated space and yielding highly stable isoform-level predictions even under extreme data scarcity.

3.3 Ablation study

Figure 2 illustrates the performance changes of SpliceEM on the human dataset after systematically removing four key architectural components. In this experiment, we utilize AUPRC as the primary evaluation metric. Unlike AUROC, which can present an overly optimistic assessment in highly imbalanced datasets due to the dominance of true negatives, AUPRC focuses exclusively on the predictive distribution of the positive class. Thus, it serves as a far more rigorous indicator of the model’s actual discriminative power in capturing rare, isoform-specific functional signals.

Display Full Size

Download

Figure 2. Ablation study of SpliceEM performance. ESM: evolutionary scale modeling; EMA: exponential moving average; ASL: asymmetric loss.

The ablation results reveal that removing the margin-based asymmetric loss (w/o ASL) triggers the most significant systemic drop in AUPRC. This underscores its role as a key component for addressing extreme class imbalance in isoform-level functional annotations. By applying a gene-masking strategy, the ASL effectively suppresses gradient perturbations from the vast majority of clearly unrelated GO terms, allowing the model to precisely establish robust decision boundaries for rare positive features.

Similarly, removing the teacher-student distillation framework (w/o EMA) leads to a substantial 13.55% decrease in performance. This drop validates the necessity of introducing self-supervised signals in sparse annotation scenarios. The EMA mechanism integrates consistent historical knowledge across the training cycle to build a stable teacher model, enabling safe knowledge transfer into the unannotated sample space while effectively filtering stochastic training noise.

Furthermore, replacing the pre-trained ESM-2 features with traditional sequence embeddings (w/o ESM-2) degrades performance by 11.74%. This confirms that the evolutionary semantic priors provided by advanced protein language models are crucial for establishing a baseline understanding of global structures and functional motifs, which the localized splicing event encoder subsequently refines and decouples.

Finally, removing the hierarchical triplet loss (w/o Hierarchical) causes an 8.81% performance loss. Although its direct contribution to the flat AUPRC metric is the smallest among the four components, its fundamental necessity was previously corroborated by our S_min evaluation in the main results. Rather than merely boosting raw classification precision, this constraint acts as a vital structural regularizer. It explicitly internalizes the topological priors of the DAG, ensuring that the model’s predictions strictly adhere to biological hierarchical logic and that they prevent structural conflicts in the semantic space.

3.4 Splicing-function analysis

Beyond isoform function prediction, this study explores the feature weights captured by the attention mechanism within the deep graph neural network to decode the latent patterns by which AS drives protein functional divergence. We extracted attention scores from the event-to-isoform relational edges in the final layer and analyzed the model’s predictive preferences across three scales: transcriptome-wide patterns, signaling pathway enrichment, and single-gene cases.

Transcriptome-Wide Splicing-Function Attention Mapping: First, we aggregated the average attention weights of different splicing types for various GO BPs to construct a global, data-driven attention heatmap (Figure 3). By filtering out semantically sparse nodes, we focused on the most active functional terms. The results demonstrate a striking functional division among distinct splicing types. While events like A5 and RI show highly localized attention in basic metabolic processes, SE and AF exhibit profound and broad-spectrum regulatory weights across the transcriptome. Remarkably, in critical cellular survival and signaling functions, such as the apoptotic process (GO:0006915) and monoatomic ion transport (GO:0006811), AF and SE emerge as the nearly exclusive drivers, with most other splicing types contributing completely silent signals. Since apoptosis evasion and transmembrane receptor trafficking are essential processes in cancer pathogenesis and cellular adaptation^[37], this striking dominance naturally prompts us to further investigate the specific cell signaling networks driven exclusively by SE and AF.

Display Full Size

Download

Figure 3. Attention weight heatmap mapping splice event types to GO function predictions at the transcriptome level. GO: gene ontology. A3: alternative 3’ splice site; A5: alternative 5’ splice site; AF: alternative first exon; AL: alternative last exon; MX: mutually exclusive exons; RI: retained intron; SE: skipping exon.

SE/AF-Specific Biological Pathway Enrichment Analysis: Having identified the prominent and exclusive role of SE and AF in dictating key functional differentiation, we analyzed the distribution of target genes heavily driven by these two events. To control for statistical false positives, we performed Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis exclusively on the high-attention target genes against a custom multi-isoform background (Figure 4). Specifically, high-attention SE/AF-driven genes were defined as host genes linked to SE or AF events with an Attention_Score > 0.1 and a Max_GO_Prob ≥ 0.5. These specific thresholds were both algorithmically and biologically motivated, as a Max_GO_Prob ≥ 0.5 serves as the standard binary classification boundary for a positive GO annotation, while an Attention_Score > 0.1 signifies a remarkably strong regulatory contribution, given the extensive number and diversity of semantic neighbors competing for attention weights within the heterogeneous graph. After mapping isoform-level records to host genes and removing duplicates, this procedure yielded 74 unique genes in the current dataset. These genes were tested for KEGG enrichment against a background universe of 7,433 mapped genes, with pathways ranked by false discovery rate adjusted $p$-values.

Display Full Size

Download

Figure 4. KEGG pathway enrichment analysis of target genes driven by high-attention SE/AF events. KEGG: Kyoto encyclopedia of genes and genomes; SE: skipping exon; AF: alternative first exon; MAPK: mitogen-activated protein kinase; JAK: Janus kinase; STAT: signal transducers and activators of transcription.

The analysis reveals that target genes significantly driven by SE and AF are profoundly enriched in classic oncogenic and intracellular signaling cascades, prominently including Pathways in cancer, the PI3K-Akt signaling pathway, the MAPK signaling pathway, the Ras signaling pathway, and the JAK-STAT signaling pathway. To rigorously verify that this observed enrichment could not be explained by background expectation, we performed 100 matched random gene-set controls by repeatedly sampling 74 genes from the same background universe. The top enriched cascades remained significantly more prominent than the random controls, yielding an empirical P = 0.0099 for the leading pathways. This enrichment perfectly matches known tumor adaptive evolution strategies. In cancer pathogenesis, cells frequently exploit SE and AF events to truncate negative regulatory domains or alter kinase interaction interfaces, thereby sustaining the hyperactivation of these exact RTK/MAPK and PI3K cascades^[38]. The model’s unsupervised attention mechanism independently highlights these classic pathological pathways, strongly suggesting its value in identifying actionable splicing targets in oncology.

Predictive Functional Divergence Validated by Context-Specific Expression: To investigate whether the functional divergence predicted by our model manifests in biological contexts, we performed cross-validation on target genes representing distinct regulatory patterns. The aspartate-β-hydroxylase gene (ENSG00000198363) serves as a primary example of an antagonistic isoform switch (Figure 5). Our model predicts a clear functional dichotomy between its two main transcripts. Isoform 1 (ENST00000379454) is predicted to govern nuclear and epigenetic events, including chromosome organization (GO:0051276) and mRNA transcription (GO:0009299). In contrast, the predicted functional spectrum for Isoform 2 (ENST00000517856) shifts toward membrane transport processes, notably calcium ion transmembrane transport (GO:0070588) and monoatomic cation transport (GO:0098655).

Display Full Size

Download

Figure 5. Functional divergence and context-specific expression profiles of the ASPH gene. The panels illustrate the predicted functional dichotomy and the antagonistic expression switch of the two primary transcripts across different biological environments. ASPH: aspartate-β-hydroxylase.

Examining their expression profiles across different states reveals an inverted regulatory pattern. Isoform 1 demonstrates higher expression in rapidly proliferating cell lines, whereas Isoform 2 becomes the dominant transcript in specific differentiated environments, such as kidney tissue. This functional divergence is driven by AS mechanisms that alter the protein structure, shifting its biological role from a nuclear regulator to a transmembrane transporter.

4. Discussion

The SpliceEM framework offers an innovative framework to address two primary challenges in isoform function prediction: multi-modal feature fusion and extreme annotation sparsity. Compared to previous approaches, SpliceEM offers a more detailed perspective on how transcriptomic functional diversity is modeled. Early MIL methods and shallow network algorithms (such as IsoFun and DisoFun) successfully leveraged broad gene-level information; however, they often struggled to capture complex biological hierarchies. Recent deep graph models like IsofunGO and CrossIsoFun have significantly improved representation capacity, but they typically encode highly similar isoforms as whole sequences. This whole-sequence approach sometimes fails to capture the specific functional variations caused by local exon changes. SpliceEM addresses this issue by integrating isoforms, genes, AS events, and GO terms into a multi-layer HGT architecture, using cross-modal attention to separate detailed functional signals from the overall sequences.

To navigate the extreme sparsity of isoform annotations and address the limitations of previous models in predicting rare functional terms, SpliceEM utilizes a teacher-student distillation framework guided by EMA and MIL. This mechanism synergistically fuses broad gene-level priors with the limited but highly precise isoform ground truths, safely guiding the exploration of the unannotated space while ensuring stable training. Additionally, combining gene-masked ASL with hierarchical triplet constraints effectively filters negative gradient interference and enforces GO topological logic. Our experimental results show that these strategies yield significant improvements in both AUPRC and S_min, establishing highly precise decision boundaries while ensuring ontologically coherent predictions.

Beyond improved computational metrics, the interpretability of SpliceEM facilitates practical downstream biological analysis. By extracting attention weights from event-to-isoform relational edges, we identified key splicing-function relationships, highlighting the prominent roles of SE and AF events in oncogenic signaling cascades. Furthermore, the single-gene validation of context-specific functional switching, which is supported by cross-sample expression profiles, demonstrates the model’s utility in exploring how AS drives adaptive functional reshaping.

5. Limitations

Despite its robust performance, the application of SpliceEM faces certain practical boundaries. First, its feature extraction relies heavily on pre-trained ESM-2 embeddings, which are primarily optimized on canonical proteins. This inherent bias may limit the model’s sensitivity to highly localized, non-canonical splicing micro-variations. Furthermore, predictions remain constrained by the historical annotation bias and incompleteness of the GO database. Second, while EMA distillation effectively mitigates label sparsity, severe dropout events in low-depth transcriptomic sequencing can still degrade the broad gene-level priors required for reliable soft supervision. Third, processing full-transcriptome heterogeneous graphs via multi-layer HGT incurs high memory overhead, presenting scalability challenges for massive multi-species or pan-cancer datasets. Finally, it is important to clarify the functional scope of SpliceEM. The framework is designed to predict the intrinsic biological functions of transcript variants at an isoform-level resolution, using sequence and localized splicing event features. It does not directly ingest single-cell expression profiles as inputs to predict dynamic, cell-type-specific functional shifts of a singular isoform. However, its high-precision isoform annotations serve as a critical foundation for downstream single-cell analyses, enabling researchers to decode how selective isoform expression drives broad functional transitions across different cell types and states, as demonstrated by our context-specific functional switching validations.

6. Conclusion

In this paper, we proposed SpliceEM, a deep learning framework designed to tackle the critical challenges of multi-modal feature fusion and extreme annotation sparsity in isoform function prediction. To capture the intrinsic dependencies among different biological entities, we constructed a unified heterogeneous graph processed by multi-layer HGT, explicitly incorporating AS events as distinct topological nodes. To solve multi-resolution label imbalance, we introduced a joint optimization mechanism that synergizes an EMA-guided distillation with MIL, to safely navigate unannotated spaces, while employing gene-masked ASL and hierarchical triplet constraints to establish precise decision boundaries. Experimental results demonstrate that SpliceEM consistently achieves state-of-the-art performance across multiple evaluation metrics, even under extreme data scarcity. Furthermore, downstream analyses proved its interpretability by successfully identifying key splicing events driving oncogenic pathways and decoding context-specific functional switching at the single-gene level. Future research will focus on exploring dynamic graph sampling for scalable pan-cancer or cross-species applications, ultimately advancing its utility in clinical biomarker discovery.

Supplementary materials

The supplementary material for this article is available at: Supplementary materials.

Acknowledgements

The authors declare that Gemini was used solely for language polishing during the manuscript preparation process. All research content, including study design, data analysis, interpretations, figures, and tables, is original and was not generated using AI tools.

Authors contribution

Gu T: Data curation, software, writing-original draft.

Wang J: Conceptualization, supervision, writing-review & editing.

Conflicts of interest

Jun Wang is an Editorial Board Member of Computational Biomedicine. The other author declares no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The human isoform expression dataset used in this study is available from the Gene Expression Omnibus (GEO) repository under accession number GSE212945 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE212945). The code is available at https://www.sdu-idea.cn/codes.php?name=spliceEM.

Funding

This work was supported by the Taishan Scholars Program (Grant No. tsqn202408317).

Copyright

© The Author(s) 2026.