MultCPM: a multi-omics cancer recurrence prediction model utilizing a multi-head attention mechanism

MultCPM model framework

MultCPM, a deep learning model for cancer recurrence prediction, integrates clinical and multi-omics data, as depicted in Fig. 1. The model comprises four components: a biological pathway module, a pathway multi-head attention module, a hierarchical fusion module and a recurrence prediction module. First, the preprocessed multi-omics data, comprising mRNA, miRNA, and SNV, are associated with their respective biological pathways. Subsequently, three relationship matrices (mRNA-pathway, miRNA-pathway, and SNV-pathway) are constructed to extract key features. Next, three independent pathway modules are constructed using a multi-head attention mechanism to capture pathway features within each sample. Finally, the learned pathway features are hierarchically fused, and the fused features are input into a multilayer perceptron model to perform cancer recurrence prediction and model interpretability analysis. The model enhances prediction accuracy by capturing correlations between multi-omics using the described method.

Biological pathways module

As shown in Fig. 1A, the biological pathway module of the MultCPM model consists of an omics feature layer and a pathway layer. Using a feedforward neural network, this module maps multi-omics data to a functional pathway space, establishing structured associations between omics features and known biological pathway nodes. The design of this module builds upon existing research (Hao et al., 2018; Zhou et al., 2024a), aiming to fully leverage biological knowledge graphs to guide feature modeling and enhance the biological plausibility and interpretability of feature learning. Let $X^i\in {\mathbb{R}}^{N\times F}$ represent the raw feature matrix of the $i$-th type of omics data (mRNA, miRNA, or SNV) after preprocessing, where $N$ denotes the number of samples and $F$ denotes the number of features. We then construct a sparse binary mapping matrix $A^i\in \{0,1{\}}^{N\times L}$ to characterize the association between all features of this omics type and biological pathways, where $L$ denotes the total number of pathways (i.e., the number of pathway nodes in KEGG). The mapping matrix $A^i$ is defined as follows:

(1) $$A_{k,p}^i=\{\begin{array}{c}1,if{x}_k^i\in p\\ 0,otherwise\end{array}.$$

Here, $x_k^i$ denotes the $k$-th feature of the $i$-th type of omics data, and $p$ denotes a node in a pathway. When feature $x_k^i$ is annotated in the biological knowledge graph as a component of pathway $p$ (e.g., labeled as a member gene of the pathway in KEGG), the corresponding relationship matrix element $A_{k,p}^i$ is assigned a value of 1, indicating a structural association between the feature and pathway node $p$. Otherwise, it is assigned a value of 0, indicating no association between the two.

To further illustrate the construction process, a diagram of the mapping logic is provided in the Supplemental Materials, using the relationship between mRNA features and pathways as an example to demonstrate how matrix $A^i$ is generated.

Next, the feature propagation process from the omics feature layer to the pathway layer is described by the following formula:

(2) $${\mathrm{P}}^{\mathrm{i}}={\mathrm{X}}^{\mathrm{i}}{\mathrm{A}}^{\mathrm{T}}+\mathrm{b}.$$

In this formula, $p^i\in {\mathbb{R}}^{N\times L}$ represents the pathway feature representation corresponding to the $i$-th type of omics data, $X^i\in {\mathbb{R}}^{N\times F}$ is the input omics feature matrix, $A^i\in \{0,1{\}}^{N\times L}$ is the association matrix between omics features and pathway nodes, and $b\in {\mathbb{R}}^L$ is the bias term. This formula performs a linear mapping from the original omics feature space to the pathway space, aiding in the identification of omics information related to specific biological pathways.

Using the biological pathway module, relationship matrices were constructed between mRNA, SNV, and miRNA data and their corresponding pathway nodes, denoted as $p^{mRNA}\in {\mathbb{R}}^{N\times F_{mRNA}}$, $p^{SNV}\in {\mathbb{R}}^{N\times F_{SNV}}$, and $p^{miRNA}\in {\mathbb{R}}^{N\times F_{miRNA}}$, where $N$ denotes the number of samples, and $F_i$ denotes the dimension of the $i$-th type of omics feature after biological pathway mapping. In this study, the feature dimensions for each group are provided in Table 1. Each row of the relationship matrix corresponds to a sample, while the columns represent the feature representations of that sample across different pathway nodes, including gene expression levels, mutation status, and regulatory features, among other biological information.

Pathway multi-head attention module

To explore the potential correlations between samples, this study introduces the multi-head attention mechanism (Vaswani et al., 2017) to model the complex interactions between multi-omics data (mRNA, SNV, and miRNA) and biological pathway nodes. This mechanism captures biological dependencies at multiple levels by mapping input features to several subspaces, thereby extracting potential knowledge comprehensively. Multiple attention heads can simultaneously focus on the biological regulatory relationships between different omics features and their corresponding KEGG pathway nodes. Key genes in these pathways play a critical role in tumor progression and recurrence. Effective modeling of these dependencies aids in revealing potential recurrence drivers and enhances the model’s ability to comprehend the complexity of cancer biology.

As shown in Fig. 1B, the multi-head attention module consists of three parallel multi-head attention mechanisms, with inputs being pathway-level omics features mapped by the biological pathway module. These features incorporate KEGG pathway knowledge during construction, establishing a structured association between omics features and pathway nodes. Based on this, the attention mechanism effectively models feature dependencies between samples within the pathway space, deeply explores potential functional module synergies, and further enhance the model’s ability to represent key pathway features and provide biological interpretability.

For each attention mechanism, the feature matrix is linearly transformed to generate the query (Q), key (K), and value (V) matrices. The query identifies relevant pathway features, the key provides matching criteria, and the value contains the corresponding feature information. Each head computes the attention weights (W), which represent the correlations between different pathways, by calculating the dot product of the query matrix (Q) and the key matrix (K). The attention weights indicate the degree to which one pathway feature influences another. Higher attention weight implies that the features of one pathway are more dependent on those of another pathway.

Each attention head computes queries, keys, and values in distinct subspaces using independent linear transformations, allowing different heads to focus on various feature relationships, as expressed by the following formula:

(3) $$Attention\left(Q,K,V\right)=softmax\left({\displaystyle \frac{Q{K}^T}{\sqrt{d_k}}}\right)V$$where $d_k$ is the dimension of the key matrix, used for normalization.

The features of each attention head are derived through multi-head processing, with the $j$-th head computed as follows:

(4) $$hea{d}_j=Attention\left(Q{W_j}^Q,K{W_j}^K,V{W_j}^V\right).$$

Here, ${W_j}^Q$, ${W_j}^K$, ${W_j}^V$ represent the query, key, and value transformation matrices of the $j$-th head.

Finally, the outputs from each head are concatenated and subjected to a linear transformation to produce the final multi-head attention output. This output represents the combined result of all attention heads and encapsulates information about the complex interactions among multiple pathways, as calculated by the following formula:

(5) $$MultiHead\left(Q,K,V\right)=Concat\left(hea{d}_1,\dots ,hea{d}_h\right)W^o.$$

$W^o$ represents the training parameters, while $h$ denotes the number of heads. Here, we set $h$ to 16. For the impact of different numbers of heads on model performance, refer to the experimental analysis of MultCPM performance with varying attention heads in the experimental section.

Learning the correlations between samples using a multi-attention mechanism, extracting meaningful features from multiple perspectives, and generating high-dimensional feature representations that capture pathway feature correlations between samples can enhance the model’s predictive performance.

Hierarchical fusion module

In multi-omics cancer prediction tasks, capturing potential biological interactions between different omics is crucial for improving model performance and interpretability. Previous studies have shown that SNV mutations can regulate mRNA expression and stability by affecting promoter regions, splicing sites, and other mechanisms (Zhou et al., 2018). Additionally, SNVs may be located at miRNA binding sites or within the miRNA itself, thereby affecting miRNA production and function, and ultimately regulating target gene expression (Schneiderova et al., 2017; Urbanek-Trzeciak et al., 2020).

As shown in Fig. 1C, based on this biological motivation, we propose a “layered fusion strategy” to construct two typical regulatory pathways, “SNV-mRNA” and “SNV-miRNA”, and ultimately integrate them into a unified feature representation. The entire fusion process is as follows:

1. Encoding regulation fusion pathway

Concatenate the pathway features of SNV and mRNA to obtain the intermediate fusion representation $X_{SNV-mRNA}$ of the encoding layer:

(6) $$X_{SNV-mRNA}=concatenation\left(X_{SNV},X_{mRNA}\right).$$

This step simulates the regulatory mechanism of SNV on gene expression.

2. Non-coding regulatory fusion pathway

Concatenate the pathway features of SNV and miRNA to form a fusion representation of non-coding regulation $X_{SNV-miRNA}$:

(7) $${\mathrm{X}}_{\mathrm{S}\mathrm{N}\mathrm{V}-\mathrm{m}\mathrm{i}\mathrm{R}\mathrm{N}\mathrm{A}}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\left({\mathrm{X}}_{\mathrm{S}\mathrm{N}\mathrm{V}},{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{R}\mathrm{N}\mathrm{A}}\right).$$

This reflects how SNV indirectly regulates gene expression by affecting the miRNA regulatory axis.

3. Global fusion

Fuse the two intermediate representations again to obtain the final feature representation $X_{final}$:

(8) $$X_{final}=concatenation\left(X_{SNV-mRNA},X_{SNV-\mathrm{m}\mathrm{i}RNA}\right).$$

The final representation integrates both coding and non-coding regulatory pathways, enabling the model to simultaneously model the effects of mutations on both layers of regulatory mechanisms. The specific implementation process of this hierarchical fusion is further illustrated by the pseudocode as follows (Box 1):

In the above hierarchical fusion, SNV features were reused twice during the fusion process due to their role as a core upstream factor in cancer development, influencing multiple regulatory pathways simultaneously. This “hierarchical splicing” strategy explicitly models hierarchical regulatory pathways (e.g., SNV-mRNA, SNV-miRNA) structurally, helping to preserve local regulatory structures and enhancing the model’s ability to model interactions at high-level feature abstraction. Compared to direct splicing of the entire histology and bilinear fusion, this approach improves performance and biological interpretability. The performance comparison of each fusion method is presented in the “Comparison of different fusion methods” section of the experimental results.

Recurrence prediction module

As illustrated in Fig. 1D, we construct a four-layer multilayer perceptron (MLP) network for cancer recurrence prediction, utilizing the tanh activation function in the hidden layers and the Sigmoid activation function in the output layer. The objective function for MLP classification is defined as follows:

(9) $$\ell =-{\displaystyle \frac{1}{N}{\sum }_{i=1}^N\left[y_i\log ({\hat{y}}_i)+\left(1-y_i\right)\log (1-{\hat{y}}_i)\right].}$$

Here, $N$ represents the number of samples, $y_i$ is the true label of the $i$-th sample, where $y_i$ = 1 indicates recurrence and $y_i$ = 0 indicates no recurrence. ${\hat{y}}_i$ denotes the probability predicted by the model for recurrence of the $i$-th sample, which is finally output through the Sigmoid function. The features and clinical data obtained after hierarchical fusion are input into the network. After forward propagation through each hidden layer, the network outputs a recurrence probability for each patient. A higher probability value indicates a greater likelihood of recurrence, and vice versa.

Model interpretation

To evaluate the role of multi-omics features (mRNA, miRNA, and SNV) in predicting cancer recurrence, the DeepSHAP method (Chen, Lundberg & Lee, 2022; Eikså, Vatne & Lekkas, 2024) was used to analyze the interpretability of the trained MultCPM model. DeepSHAP combines Shapley value theory with the gradient decomposition mechanism of DeepLIFT (Bhattarai et al., 2024) to effectively measure the marginal contribution of each input feature in the deep neural network to cancer recurrence prediction. This method reveals the key features involved in the model’s decision-making process by quantifying their contributions to the prediction.

Figure 2 illustrates the complete process of DeepSHAP parsing the MultCPM model. First, mRNA expression data, miRNA expression data, and SNV mutation data are used as inputs, with the dataset divided into training and test sets in an 80%:20% ratio to ensure the independence of the training and evaluation processes. Subsequently, the three types of omics data are input into the MultCPM model for training, and prediction is performed on the test set to obtain the model output, which provides the basis for subsequent SHapley Additive exPlanations (SHAP) value calculation. In the DeepSHAP interpretation stage, the trained MultCPM model is combined with the test set inputs, and the SHAP value for each input feature quantifying its positive or negative contribution to the model output (i.e., cancer recurrence prediction) is calculated. Through this process, feature contribution matrices for mRNA, miRNA, and SNV data are obtained separately, revealing the decision-making basis and key drivers within the model.

In DeepSHAP’s model parsing, the Shapley value is calculated as follows:

(10) $${\varphi }_i(f,x)={\sum }_{S\subseteq N\{i\}}\frac{|S|!(|N|-|S|-1)!}{|N|!}[f(S\cup i)-f(S)]$$

$\varphi \left(f,x\right)$ represents the SHAP value of feature $i$ for model $f$ under input $x$. $S$ denotes the subset of features, and $N$ represents the full set of features. $f\left(S\right)$ is the predicted output of the model using only the features in subset $S$. This formula calculates the marginal contribution of each feature $i$ to the output when it is introduced into the feature subset $S$. The contribution value $C_i$ of the feature is calculated as follows:

(11) $$C_i={\displaystyle \frac{\mathrm{\partial }f\left(x\right)}{\mathrm{\partial }{x}_i}.\left(x_i-r_i\right).}$$

Here, $f\left(x\right)$ represents the output of the model, ${\displaystyle \frac{\mathrm{\partial }f\left(x\right)}{\mathrm{\partial }{x}_i}}$ denotes the gradient of the model at feature $i$, $x_i$ is an input feature, and $r_i$ is the reference value for feature $i$.

Through DeepSHAP’s model parsing, we not only mapped the model output back to the original feature layer but also quantified the importance of each feature in predicting cancer recurrence. This provided data support for the subsequent analysis of the biological mechanism and screening of potential markers.

Data processing

In this study, we used cancer datasets from the The Cancer Genome Atlas (TCGA) and Therapeutically Applicable Research to Generate Effective Treatments (TARGET) databases, along with the KEGG pathway dataset, to evaluate the predictive performance of our models.

First, for the KEGG pathway dataset, pathways and their associated gene sets were primarily obtained in this study using the getGenesets function from the R package EnrichmentBrowser (Geistlinger, Csaba & Zimmer, 2016), which extracted 359 KEGG pathways and their corresponding 8,842 gene entries. Meanwhile, the mirPath v.3 web server (Vlachos et al., 2015) was used to obtain the set of miRNAs associated with KEGG pathways, involving 238 pathways and 952 miRNA entries. To improve the model’s generalization ability and avoid the overfitting effect of specific cancer pathways on the training results, we preprocessed the extracted pathway data as follows:

1. Elimination of cancer-specific pathways: 62 pathways clearly labeled as related to specific cancers, such as “hsa05224-Breast cancer,” were excluded based on pathway annotations from the KEGG website, leaving 297 pathways for subsequent modeling.

2. Standardization of pathway names: Some pathways were inconsistently named or lacked detailed descriptions in different sources. For instance, the “Mitogen-Activated Protein Kinase (MAPK) signaling pathway” was sometimes listed as just “MAPK.” We standardized and supplemented the names of these pathways using the KEGG database (Kanehisa et al., 2024), resulting in the correction and improvement of 32 pathway names.

3. Gene and miRNA de-duplication: Due to duplicate annotations of genes and miRNAs across different pathways, we removed 768 duplicate genes and 140 duplicate miRNAs, resulting in 8,074 unique genes and 812 unique miRNAs.

These steps ensured the completeness, consistency, and non-redundancy of the pathway data, making the mapping of multi-omics data to pathway features more accurate and reliable.

For the TCGA datasets, we obtained data for breast cancer (BRCA), liver cancer (LIHC), and bladder cancer (BLCA) from the Xena TCGA Pan-Cancer platform (Goldman et al., 2020). These data include multi-omics data for mRNA (FPKM), miRNA, and SNV (MuTect2), as well as clinical data on patient survival and mortality. For sample selection, a rigorous data screening process was employed to ensure that each sample had complete multi-omics characterization and clinical recurrence information, thus ensuring the scientific rigor and reproducibility of model training and evaluation. The following screening process was applied to retain samples with complete mRNA, miRNA, and SNV data, clinical data with clear recurrence information, and samples with less than 10% missing data for each histology, in accordance with the methodology of a previous study (Dhillon & Singh, 2020). After this screening, the final number of retained samples is shown in the sample column of Table 2.

In this experiment, due to the different data formats of the original multi-omics data, it was necessary to normalize each histology dataset. For mRNA data, ENSGID identifiers were replaced with corresponding gene names, and only genes matching the KEGG gene set were retained, with gene expression values normalized. For miRNA data, only miRNAs associated with the KEGG pathway were retained, and their expression values were normalized. For SNV data, nucleotide variants in genes were marked as 1, and non-variants as 0. The feature dimensions before and after processing are shown in Table 2, which includes the number of raw and trained feature values. To fill missing values in mRNA and miRNA, we compared four commonly used methods: random forest, K-nearest neighbors (KNN), soft-impute, and Non-negative Matrix Factorization (NMF). The experimental results show that random forest is the optimal method for filling missing values, and the impact of this strategy on model performance is presented in Appendix B. Detailed information on the histology data of the TCGA dataset is provided in Table 2.

For the TARGET dataset, clinical data, mRNA, and miRNA expression data for acute myeloid leukemia (TARGET-AML) and nephroblastoma (TARGET-WT) were obtained from the GDC data portal (https://portal.gdc.cancer.gov/v1/). miRNA data were first processed using a log2 transformation to nonlinearly smooth their distribution, followed by normalization. mRNA data were pre-processed in the same manner as the TCGA dataset. For samples or genes/miRNAs with duplicate records, the final expression levels were obtained by averaging the values. To reduce the dimensionality of the multi-omics data and eliminate redundant features, the chi-square test was applied to screen features for each group. The missing value imputation strategy is consistent with TCGA and implemented using the random forest algorithm. The final dimensionality of the omics features in the TARGET dataset is detailed in Table 3.

Evaluation metrics

To evaluate the model’s performance in cancer recurrence prediction, we utilized standard evaluation metrics, including accuracy (Schneiderova et al., 2017), precision, recall, F1-score, Area Under the Receiver Operating Characteristic Curve (AUC-ROC), and area Under the Precision-Recall Curve (AUPR). Detailed definitions of these metrics are outlined below:

(12) $$Accuracy={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(13) $$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(14) $$\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}={\displaystyle \frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}}$$

(15) $$F1\text{-}score=2\ast {\displaystyle \frac{Precision\ast Recall}{Precision+Recall}.}$$

The predictive performance of the MultCPM model was assessed using these metrics.

Implementation details

The MultCPM model is built on the TensorFlow 2.6.0 framework and Python 3.9.13, with the primary development and runtime environment being Ubuntu 22.04. Model training is conducted on an NVIDIA RTX 4090 GPU (24 GB video memory). The optimizer used is Adam, with a learning rate of 1e−4, a batch size of 128, and 200 training epochs. The loss function employed is binary_crossentropy.

To ensure the robustness and generalization of model performance, we use 5-fold cross-validation repeated five times, with the average result serving as the final performance metric. Additionally, the number of attention heads, a key hyperparameter, is tuned. All code, data processing scripts, and running instructions are open sourced at https://github.com/dowell2016/MultCPM. The complete code conda run configuration file (README.md) is also provided.

Performance comparison with other advanced classification methods

To further assess the effectiveness of the proposed MultCPM model, its performance was compared against several widely used machine learning classifiers, such as Logistic Regression (LR) and Support Vector Machine (SVM), along with deep learning approaches like MOGONET (Wang et al., 2021a) and DeepKEGG (Lan et al., 2024). For training the LR and SVM models, the omics data connections were directly used as input features. The preprocessed data was used directly for training in the MOGONET method.

We comprehensively evaluated the performance of various machine learning and deep learning models on three cancer datasets BRCA, BLCA, and LIHC. Table 4 summarizes the comparison of LR, SVM, MOGONET, DeepKEGG, and our proposed MultCPM model across five key metrics: ACC, AUC, F1-score, AUPR, and REC.

The performance of the MultCPM model on the BRCA dataset is particularly noteworthy. The model achieves an ACC of 0.820 ± 0.004, significantly outperforming other methods, and an AUC of 0.900 ± 0.005, demonstrating clear advantages. The F1-score of 0.712 ± 0.013 is notably higher than that of DeepKEGG (0.672 ± 0.007) and MOGONET (0.596 ± 0.033). Additionally, the model achieves a REC of 0.632 ± 0.006 and an AUPR of 0.893 ± 0.006, which are significantly superior to those of DeepKEGG, MOGONET, and traditional machine learning models.

On the BLCA dataset, the MultCPM model demonstrated exceptional performance, achieving an ACC of 0.903 ± 0.004, slightly surpassing DeepKEGG 0.896 ± 0.003, and significantly outperforming traditional methods such as SVM and LR. The AUC reached 0.964 ± 0.011, markedly exceeding the performance of other models. The F1-score was 0.863 ± 0.013, outperforming DeepKEGG 0.857 ± 0.006, further emphasizing its advantage in handling imbalanced data. Additionally, the REC and AUPR were slightly higher than those of DeepKEGG and far exceeded other models, fully demonstrating its superior predictive capability.

The performance of MultCPM on the LIHC dataset is significantly improved, achieving an ACC of 0.937 ± 0.008, which is 6 percentage points higher than DeepKEGG 0.877 ± 0.006. The AUC reaches 0.979 ± 0.006, the highest among all models. The F1-score is 0.935 ± 0.005, substantially outperforming 0.870 ± 0.007 of DeepKEGG and 0.816 ± 0.011 of MOGONET. Additionally, its REC and AUPR are 0.877 ± 0.025 and 0.949 ± 0.005, respectively, fully demonstrating the superior capability of MultCPM in multimodal data integration.

In summary, MultCPM demonstrated outstanding classification performance across all three cancer datasets. Its superior results in ACC, F1-score, REC, AUC, and AUPR indicate that MultCPM effectively learns high-dimensional features in mRNA, miRNA, and SNV from multimodal data, significantly enhancing the model’s learning efficiency. By employing a hierarchical data integration approach, MultCPM fully captures the complex feature interactions among multi-omics data, substantially improving classification task performance. Notably, its exceptional performance on the LIHC dataset further validates the model’s strong adaptability and generalizability in handling heterogeneous biomedical data.

To further validate the generalization ability of the MultCPM model, this study conducted extended experiments on acute myeloid leukemia (AML) and nephroblastoma (WT) from the TARGET dataset. These experiments systematically evaluated the performance differences between MultCPM and other comparative methods in predicting recurrence across different cancer types. The results of the experiments are shown in Table 5.

The MultCPM model achieved the best performance on the AML dataset. Its ACC reached 0.778 ± 0.016, approximately 1 percentage point higher than DeepKEGG’s 0.768 ± 0.005. The AUC of MultCPM reached 0.851 ± 0.003, slightly higher than DeepKEGG’s 0.847 ± 0.001, marking the highest value among all models. Its F1-score was 0.796 ± 0.011, significantly better than DeepKEGG’s 0.790 ± 0.005 and MOGONET’s 0.718 ± 0.004, reflecting its stronger category balancing ability. In terms of recall (REC) and AUPR, MultCPM achieved 0.808 ± 0.006 and 0.883 ± 0.002, respectively, further confirming its robustness and superiority in handling high-dimensional heterogeneous data fusion.

MultCPM also delivered the best performance on the WT dataset. Its ACC was 0.849 ± 0.003, and its AUC reached 0.865 ± 0.004, both the highest values among all models. The F1-score of 0.906 ± 0.007 was significantly higher than DeepKEGG’s 0.892 ± 0.003, outperforming all baseline models. Additionally, MultCPM’s REC reached 0.920 ± 0.003, and AUPR reached 0.963 ± 0.008, demonstrating strong recurrence patient identification and stable prediction performance. These results further highlight the generalization advantage of this method in multi-omics integration.

Performance comparison of prediction for single-omics data and multi-omics integration

To further validate the effectiveness of integrating multi-omics data in cancer recurrence prediction, this study applied the MultCPM model to both multi-omics integration and single-omics input tasks. Using the BRCA, BLCA, and LIHC datasets, we systematically compared the classification performance of multi-omics integration with that of the three single-omics (SNV, mRNA, and miRNA). The results are shown on Fig. 3.

The experimental results show that in the BRCA and BLCA datasets, the multi-omics fusion model outperforms the single-omics model in terms of ACC, AUC, F1-score, REC, AUPR, and other metrics. This demonstrates that the MultCPM model effectively explores the potential complementarities between data from different groups, improving the prediction of complex diseases. However, in the LIHC dataset, the SNV single-omics model outperforms the multi-omics fusion model in some metrics, which raises our concern. To further investigate this, we evaluated the existing multi-omics model DeepKEGG on the same dataset and found that it also performs better with the SNV single-omics model than the multi-omics fusion model. We hypothesize that this may be due to the lower quality or limited information of mRNA and miRNA features in the LIHC dataset, preventing the multi-omics fusion from fully leveraging its advantages. In contrast, SNV, as the core mutation-level information, holds greater predictive value in the recurrence mechanism of LIHC cancers. This phenomenon suggests that, during multi-omics fusion modeling, the quality differences and information redundancy of data from different groups must be carefully considered. Future work will further analyze the biological characteristics of the various omics features in the LIHC dataset and their impact on optimizing fusion strategies.

Comparison of different fusion methods

To further assess the effectiveness of the proposed hierarchical fusion method, we conducted ablation experiments on three cancer datasets to systematically evaluate its performance. Specifically, the hierarchical fusion method was compared with two alternative data fusion methods: direct splicing (concatenate) and bilinear fusion (Wang et al., 2021b). To achieve bilinear fusion, this study designed a custom BilinearFusion layer, which takes as input the feature tensors of SNV and mRNA, and SNV and miRNA, processed by the multi-head attention module. The layer effectively models the bilinear interaction between the two feature sets through two tensor dot product operations. Implementation details are provided in Appendix C. To ensure a fair comparison with the hierarchical fusion approach, bilinear fusion is computed separately for SNV-mRNA and SNV-miRNA, with the fusion results of both sets concatenated to form the final input feature representation. The experimental comparison of different fusion strategies is shown in Table 6, further validating the performance differences of the respective methods in the multi-omics integration task.

The experimental results, presented in Table 6, show that the hierarchical fusion method significantly outperforms direct splicing and bilinear fusion across multiple evaluation metrics, including F1-score, ACC, AUC, REC, and AUPR. These findings indicate that the hierarchical fusion method effectively captures the potential complementarity and complex interaction between different omics data, thereby significantly enhancing prediction performance.

Performance of MultCPM with different attention heads

The number of attention heads $h$ should be balanced against model complexity and dataset size. When a smaller value of $h$ is chosen, the model has limited ability to learn complex interactions between the omics data, making it difficult to adequately capture key features. In contrast, a larger value of $h$ enhances the model’s expressive power but also increases computational complexity, potentially leading to overfitting, especially on datasets with limited sample sizes.

We conducted five 5-fold cross-validation experiments on three datasets BLCA, LIHC, and BRCA using several evaluation metrics, including F1-score, ACC, AUC, REC, and AUPR. Additionally, in conjunction with the cancer datasets (sample sizes ranging from approximately 200 to 400), we set $h$ to a small range of discrete values {8, 16, 32, 64} to evaluate model performance variation under different settings. Here, $h=8$ serves as the minimum attention-head configuration for assessing the lower-bound performance of models with limited expression capabilities. The final performance evaluation results (averaged over five 5-fold cross-validation experiments) are shown in Table 7.

The results above indicate that the F1-score, ACC, AUC, REC, and AUPR of the model follow a pattern of “initial increase, followed by stabilization, and then decrease” as the number of attention heads $h$ increases across the three datasets. When $h=8$, the model is underfitted, whereas when $h=32$ and $h=64$, the model exhibits overfitting to some extent. The best prediction performance on all three datasets was achieved when $h=16$, significantly outperforming other settings.

This finding suggests that an optimal number of attention heads can better capture dependencies between multi-omics data, thereby improving the model’s generalization and predictive power. Considering both model performance and computational cost, $h=16$ is selected as the optimal hyperparameter configuration for MultCPM.

DeepSHAP interpretability analysis

To better illustrate the contributions of individual features to the model’s prediction outcomes on the BRCA dataset, we visualized the model using SHAP values. By generating the SHAP Heatmap, we highlighted the distribution of SHAP values for the top mRNA, SNV, and miRNA features and examined how variations in these feature values affect the prediction results, as depicted in Fig. 4.

The analysis results revealed that MAP3K1, MYLK, and MAPK10 in the mRNA profile; ADCY5, ATP1A4, and IL6 in the SNV profile; and hsa-miR-363 and hsa-miR-20b in the miRNA profile had high SHAP values, indicating their significant impact on cancer recurrence prediction. The high SHAP values of these genes suggest that they play crucial positive or negative roles in prediction, thereby supporting the biological interpretability of the model.

The biological relevance of the high SHAP-value mRNAs, miRNAs, and SNV features identified in relation to cancer recurrence was validated through DeepSHAP’s interpretive analysis. For instance, the MAP3K1 and MYLK genes are well-documented as key regulators of cancer development and recurrence, and DeepSHAP analysis further highlights their significance in our predictive models (Liu et al., 2018). Similarly, hsa-miR-363 and hsa-miR-20b have been strongly associated with metastasis and recurrence across various cancer types, suggesting their critical regulatory roles in cancer recurrence mechanisms (Hui et al., 2020; Wang, Yang & Xiao, 2016). Additionally, the analysis of SNV loci provides valuable references for mutation-related studies, further underscoring their potential significance in cancer recurrence research (Kandoth et al., 2013).

In summary, DeepSHAP analysis not only provided a reliable basis for interpreting the model’s prediction results but also revealed the key molecules and their pathway mechanisms that the model relies on at the biological level, enhancing the model’s interpretability and potential for application in clinical cancer recurrence risk assessment.