1 Introduction↩︎

Large-scale analyses of molecular profile data have the potential to transform cancer genomics and are pivotal for the advancement of precision medicine and personalised therapies [1], [2]. However, the integration of high-dimensional, heterogeneous molecular data, such as genomics, transcriptomics, and proteomics, pose significant challenges, particularly in balancing predictive accuracy with interpretability in machine learning models [3], [4]. While deep learning has demonstrated remarkable success in drug response prediction [5], [6], phenotypic classification [7], [8], disease subtype classification [9], [10], and survival prognosis [11], its translational impact is often limited by the inherent opacity of its black-box representations [3], [4].

Interpretability is critical in medical applications, where decision-making must be transparent and biologically plausible [12]. Post-hoc explainability methods (e.g. [13], [14]) explain individual features retrospectively, but do not constrain the internal factors in the model and often require re-analysis for each prediction or dataset [15]–[17], which limits their utility in clinical applications. A recent shift has been the incorporation of biological prior knowledge into neural network architectures, giving rise to so-called Biologically Informed Neural Networks (BINNs), also called visible neural networks (VNNs). These models embed known structures, such as biochemical pathways or gene-sets, into the model architecture, thereby enhancing interpretability and potentially improving generalisability. Recent studies ([3], [4]) highlight the ability of BINNs to bridge the gap between predictive performance and mechanistic insight, with models that are explainable by design and do not lag in performance to traditional models.

The models developed recently [18] fall inside the class of BINNs that use Biochemical pathways as the main prior-knowledge source for our architectural constraints, which, in this paper, we use to provide an avenue to explore the multi-omic interactions from a pathway-centric perspective. The main contributions of this paper are:

2 Materials and Methods↩︎

In this section, we briefly outline the PAAE framework in Subsection 2.1, and highlight recent methodological advances over previous work [18], as illustrated in Figure [fig:diagram]. We then detail the repeatability and robustness analyses (Subsection 2.3), describe the multi-omics integration strategies employed (Subsection 2.4), and present the methodology used for identifying individual-layer contributions to the downstream tasks (Subsection 2.5). Finally, we summarise the settings used for the downstream tasks in Subsection 2.2.

2.1 The PAAE Framework↩︎

The PAAE framework [19], shown in Figure [fig:diagram], is a Biologically-Informed Neural Network [4] framework based on the Autoencoder [20] and Variational Autoencoder (VAE) deep learning architectures. It works by incorporating pathway information in its architecture, building one pathway encoder block \(E_{p}\) for each pathway \(p\) in a pathway set \(P\). Each of these pathway encoders receive as input \(x_{p}\) only the molecular features which map to a gene that belongs inside the pathway \(p\), and then output a single pathway activity score \(a_{p} = E_{p}(X_{p})\). Thus, we generate a pathway activity space \(a\) by concatenating \(\|\) all individual pathway activities, as in Equation 1

\[\label{eq:paae-a} a = \left\|_{p\in P}a_{p}\right. = \left\|_{p\in P}E_{p}(x_{p}) \right.\tag{1}\]

While in the original formulation, this pathway activity space is further encoded into a latent space \(z\), given that we are interested in only in the interpretable pathway activity space, we can consider that \(a\) gets further decoded by a decoder \(D\), which reproduces the full molecular input \(x\). That is, we have a reconstruction \(\hat{x} = D(a)\), which is then used to optimise the network parameters for all \(E_{p}\), \(\Omega\), and for \(D\), \(\Gamma\) using the a distance function \(d\). In practice, we use the mean squared error as \(d\), since we normalise our input, giving us the optimisation objective in Equation 2 .

\[\label{eq:paae-opt} \min_{\Omega,\Gamma}\left(\hat{x}-x\right)^{2}\tag{2}\]

2.2 Downstream Tasks↩︎

For classification, we used the best-performing classifier found during internal validation on the RNAseq study (i.e. Logistic Regression LR), using as input either the latent space \(z\) or pathway activity space \(a\), to predict the cancer subtype. Model quality is assessed with usual classification metrics (Accuracy, Precision, Recall, and the Area Under the Receiver Operating Characteristic Curve, ROC AUC), macro-averaged and, in the case of ROC AUC, one-vs-rest (OvR) was used. Cancer subtype identification is also performed in an unsupervised fashion through clustering and clustering performance is evaluated using Mutual Information (MI) between the generated clusters and the target classes. Additionally, MI is also used to analyse latent and input features without training a classifier, either by measuring the MI between the feature and the discrete class label (as was done in the RNAseq-only study), or between a feature and a binary indicator variable for OvR analyses.

For survival analyses, we build univariate Cox Proportional Hazard (Cox-PH) models to assess whether individual features are significantly associated with survival. When a feature is significant, values are stratified into tertiles (low, medium and high), and logrank tests are performed between the low and high groups. In cases where two datasets are compared, the consistency of the direction of survival association (i.e. sign matching between low and high groups) is also taken into consideration. Additionally, we may use the Concordance-Index (CI) of a Cox-PH model [21] as a ranking metric for survival relevance, analogous to how MI is used for classification tasks. For further details, see Supplementary Material Sec. 11.

2.3 Repeatability and Robustness Analyses↩︎

For robustness and repeatability analyses, we measure the Pearson Correlation between neurons across different repeated random initialisation and training runs, following work reported in [22] and [23]. We expand on this work by incorporating other correlation metrics, such as Kendall’s \(\tau\)-b coefficient [24] and Centered Kernel Alignment (CKA, [25], [26]) similarity as a measure of overall repeatability, by analysing both the aggregate ranking of downstream-task metrics and the ranking agreement between models through Average Precision K (\(AP@K\) [27], and by performing experiments building aligned consensus models from multiple repeats.

2.3.1 Aligned Consensus Models↩︎

When performing repeatability and robustness analyses, we calculated Aligned Consensus Models (ACM) in order to study whether they were more repeatable than our non-consensus models. To build an ACM, one takes the pathway activity vector \(a_{j}\) of a repeat \(j \in [1,\ldots,k]\) out of \(k\) repeats. We then take the first model \(a_{1}\) as a reference, and calculate the aligned average pathway activity vector \(\overline{a}_{i}\) for pathway \(i\) as in Equation 3 .

\[\label{eq:aligned-consensus-model} \overline{a}_{i} = \frac{1}{k} \sum_{j\leq k} \text{sgn}(\text{corr}({a}_{j_{i}},{a}_{1_{i}})){a}_{j_{i}}\tag{3}\]

That is, we take the average of the pathway activity vector \(a_{j_{i}}\) for pathway \(i\) and repeat \(j\), multiplied by the sign \(\text{sgn}(c_{{1,j}_{i}})\) of a correlation metric \(\text{corr}\) between it and the reference model’s pathway activity vector \(c_{{1,j}_{i}} = \text{corr}({a}_{j_{i}},{a}_{1_{i}})\). In practice, we take the Pearson Correlation Coefficient for \(\text{corr}\).

2.3.2 Ranking Agreement↩︎

To compute the ranking agreement between two individual repeats \(j,m \in [1,k]\), we take a metric \(f(a,y)\) that takes a pathway activity vector \(a\) and a target \(y\), and computes a score \(s_i\) for each pathway \(i\) w.r.t. to that metric. We then take the two rankings \(R_j\) and \(R_m\), where \(R_j[r] = i\) means that pathway \(i\) in repeat \(j\) was the \(r\)-highest score among the pathways. We calculate the Average Precision K (\(AP@K\)) metric [27], modified to only calculate relevant items up to \(K\) as shown in Equation 4 .

\[\label{eq:acc-at-k} AP@K_{j,m} = \frac{1}{K} \sum_{r<K} P(j) \cdot \mathbb{1}(R_{j}[r] \in R_{m}[:K])\tag{4}\]

Where \(P(j)\) is the precision at rank \(j\), \(\mathbb{1}\) is an indicator function for whether the item is in the relevant items, and \(R_{m}[:K]\) are all items from the beginning up to \(K\) in ranking \(R_{m}\). Note that, we assume that \(K \leq |P|\).

2.4 Early and Late Multi-Omics Pathway Activity Inference and Analysis↩︎

Multi-omics integration is categorised into sequential, late and early (joint) integration methods [10]. Late integration methods analyse each omics layer independently before combining results, capturing shared patterns across layers. Early integration methods analyse all layers simultaneously, enabling detection of complex inter-omic interactions. In this section, we extend a single-omics framework to multi-omics using both early and late pathway-based integration. We will consider \(|O|\) views (i.e. Omics layers), where each omics profile \(o\) is represented as \(x_{o}\), assuming pre-processing has been completed.

In our early integration approach, we combine all omics profiles into one comprehensive profile, \(x\) (Equation 5 ). Each pathway encoder, \(E_{p}\), then specifically focuses on the molecular features linked to the genes (\(G_{p}\)) that are relevant to its particular pathway (\(p \in P\)).

\[\label{early-input} x = \|_{o \in O} x_{o}\tag{5}\]

PAAE/PAVAE is optimised to reconstruct the unified profile \(x\), enabling each pathway activity score \(a_{p}\) to capture complex cross-omics interactions and jointly represent correlated features while reconstructing the full input.

For late integration, we train \(|O|\) separate PAAE/PAVAE models with distinct parameters, each producing pathway activity scores \(a_{o}\) (Equation 6 ). If an omics layer lacks genes for pathway \(p\), we exclude it from that model, denoted as \(p \not\in A_{o}\).

\[\label{eq:omic-pathway-activity} a_{o_{p}} \approx E_{o_{p}}(x), E_{o_{p}}(x) = ({E_{o_{p}}}_1 \cdot \ldots {E_{o_{p}}}_k)(x_{o})\tag{6}\]

After generating pathway activities from each model, we integrate using two methods: late-mean integration, and late-concat integration. Late-mean integration takes the arithmetic average of the pathway activities of each model to yield an integrated pathway activity vector \(\overline{a}\) (Equation 7 ):

\[\label{eq:late-mean} \overline{a}_{p} = \frac{1}{\sum_{o \in O, p\in A_{o}}1} \sum_{o \in O, p\in A_{o}} a_{o}\tag{7}\]

For late-concat integration, individual pathway activity vectors from each omics layer are concatenated, producing an integrated pathway activity vector \(\overline{a}\) (Equation 8 , where \(\|\) denotes concatenation along features).

\[\label{eq:late-concat} \overline{a} = \|_{o \in O} a_{o}\tag{8}\]

2.5 Individual Omic Layer Contributions↩︎

To evaluate the contribution of each omics layer on downstream performance, we adapt the Shapley value \(\varphi_{i}(v)\) [28], which approximates a player’s average marginal contribution in a cooperative game.

The definition of the Shapley value \(\varphi_{i}(v)\) for a player \(i \in N\) in a game with \(|N|\) players, is calculated by a weighted average of the value function \(v\), measuring the the outcome of the game when a subset players that excludes \(i\), \(S \subseteq N \setminus \{i\}\) participate, and subtracting it from the value if player \(i\) had joined the game, as is shown in Equation 9 below.

\[\label{eq:shapley} \varphi_{i}(v) = \frac{1}{n} \sum_{S \subseteq N \setminus \{i\}} \binom{n-1}{|S|}^{-1} (v(S \cup\{i\}) - v(S))\tag{9}\]

We then adapt this idea to the concept of a multi-view (i.e., multi-omics) machine learning problem, similarly to that in Shapley regression values [29]. Given a set of views (omics layers) \(O\) and a performance metric \(v\), we calculate the average marginal contribution a view \(o\) as a function \(f\) of the difference between the contribution of every subset of features \(S\) such that \(S \subseteq O \setminus \{o\}\), as in Equation 10 .

\[\label{eq:marginal-contrib} \varphi_{i}(v) = f(\{(v(S \cup\{o\}) - v(S)) \forall S \subseteq O \setminus \{i\}, S \not= \emptyset\})\tag{10}\]

In our specific case, we use \(f\) as the median value of all marginal contributions, which we define using a cross-validated metric for \(v\), such as ROC AUC or Concordance-Index, meaning we train a model based on view combination \(S\) and \(S \cup \{o\}\). We ignore the dummy case (i.e. where there is no input, \(S=\emptyset\)). We also display the distribution of the marginal contributions in plots.

3 Experimental Setup↩︎

3.1 Repeatability Studies↩︎

To assess model repeatability, which in this context refers to the ability of our models to produce consistent results when run multiple times under identical conditions on the same dataset, we use the MO-TCGA-BRCA dataset instead of the single-omics TCGA-BRCA dataset, as this work is intended to support downstream multi-omics analyses. All models were trained using the same hyperparameters established during internal validation of the single-omics analysis. To ensure computational efficiency, we restricted our experiments to PAAE models trained on the Hallmark Genes pathway, which exhibited the lowest memory and runtime requirements, thereby enabling more extensive evaluation. We build models for pathway activity dropout values of \(Pr(\text{dropout}(a)) \in \{0,0.1,0.2,\ldots0.9\}\), with the highest dropout level resulting in an average of approximately five active pathways per training epoch. For each dropout level, 16 repeats per aligned consensus model were conducted and 16 consensus models were built, yielding \(k=256\) models per dropout value, and thus 2560 models in total.

We employ scikit-learn v1.0.2 for mutual information (MI), using feature_selection.mutual_info_classif to estimate the mutual information between pathway activities and PAM50 subtypes, and metrics.mutual_info_score for one-vs-rest (OvR) MI. For survival analysis, we use lifelines v0.27.0, specifically CoxPHFitter.concordance_index_ to compute the Concordance Index (CI) of individual pathways or models. Correlation metrics are computed using scipy v1.8.1 stats.pearsonr for Pearson’s correlation coefficient, stats.kendalltau for Kendall’s \(\tau\)-b Coefficient, and stats.rankdata to derive pathway rankings \(R_{j}\). Ranking agreement is assessed using an internally developed implementation of average precision at rank \(AP@K\) (Equation 4 ).

3.2 Interpretable Pathway Activity-Based Subtype Classification↩︎

For our intrepretability study, we use a PAAE model based on the KEGG pathway set previously trained on the single-omics TCGA dataset. Given that, our goal is to propose a highly interpretable model to aid clinical application for subtype classification, we first identify the top pathway for each PAM50 subtype based on the highest OvR Mutual Information, and restrict our feature set to these selected pathways. We then generate UMAP representations of the resulting feature space, and produce featuremaps as visual representations of pathway activity. Subsequently, we train a decision tree classifier on these features, constraining the maximum number of leaf nodes to the number of subtypes. To address class imbalance and ensure that each subtype is represented in the terminal nodes, we apply Synthetic Minority Oversampling Technique (SMOTE) [30] prior to training the decision tree.

The decision rules derived by the decision tree nodes are interpreted as logical expressions. These logical conditions are used to subset samples within the UMAP space and identify approximate regions associated with each subtype. For each leaf node, we compute kernel density estimates (KDEs) over the UMAP embeddings by taking the union of all decision path conditions leading to that node. This approach enables us to visualize the spatial concentration of samples corresponding to specific subtypes in the reduced-dimensional space. We provide these visualisations as well as confusion matrices and ROC curves for our models on both the training (TCGA) and the external validation/test (Metabric) datasets, to demonstrate that this simplified, fully-interpretable model, is capable of achieving reasonable classification performance.

3.3 Multi-Omics Analyses↩︎

The MO-TCGA-BRCA dataset is employed for multi-omics analysis, and a PAAE model is trained with the same hyperparameters selected in the internal validation of single-omics analysis, with the Hallmark Genes pathway. We train models for every possible combination of the 6 omics layers available (63 combinations in total). External validation is not performed in this setting, as the Metabric dataset lacks several omics layers available in the TCGA dataset and is therefore incompatible with multi-omics integration. To ensure stability in survival analysis, features with variance lower than 0.01 over survival statuses are excluded¹, and apply an \(l^2\) norm of 0.001.

Mapping from gene-based omics layers (such as gene expression, copy number variation, mutation and RPPA) to biological pathways is straightforward. For methylation, we used the CpG site-to-gene mapping provided by the manufacturer for the Illumina HumanMethylation450 BeadChip. For the \(\mu\)RNA layer, we used the miRTarBase2025 miRTarBase_SE_WR database of \(\mu\)RNA target genes which includes target gene interactions supported by strong experimental evidence [31]. Alternative strategies were considered and described in Supplementary Material 6.

4 Results and Discussion↩︎

4.1 Feature Space Visualization for Intepretability Analysis↩︎

With respect to interpretability, pathway activity scores enable comparisons across subtypes, supporting clinically relevant insights (see discussion in Supplementary Material Section 8). As shown in Fig. 1, clustering of pathway activity scores with cosine distance yields clear separation aligned with clinical labels, demonstrating that augmenting a model with pathway prior knowledge enhances both interpretability and predictive power, even in an unsupervised setting. Furthermore, mutual information analysis shows that most of the 32 pathways identified as most informative in the training set remain informative in the test dataset, highlighting the robustness of the learned representations.

Using the KEGG pathway set, the top 5 pathways with highest one-vs-rest mutual information for each PAM50 subtype (Basal, Her2, Luminal A, Luminal B, and Normal) in the TCGA dataset were Sphingolipid Metabolism (SM), ABC Transporters (ABC-T), Valine Leucine and Isoleucine Degradation (VLID/VLI DEG), Glycosaminoglycan Biosynthesis Chondroitin Sulfate (GBCS) and Ether Lipid Metabolism (ELM). In the Metabric dataset, Tryptophan Metabolism, P53 Signaling Pathway, Glutathione Metabolism, and Mismatch Repair appeared among the top pathways instead of ABC-T, VLID, GBCS, and ELM. However, the TCGA pathways remained highly ranked at 8th, 3rd, 19th and 2nd place out of 186 total pathways, respectively, in the Metabric dataset indicating strong cross-cohort consistency. We show the feature map of this representation as in the 2D UMAP projection of our pathway activity space, as shown in Figure 2.

Key pathways prioritised by our methodology, i.e. Sphingolipid Metabolism [32], [33], ABC-Transporters [34], [35], GBCS [36], [37], and ELM pathways [38], [39], have been documented in tumorigenesis and cancer progression, either specifically in breast cancer [32]–[34], [37], [39], or more broadly across various cancer types [32], [35], [36], [38]. Importantly, we note that the identified Valine Leucine and Isoleucine Degradation (VLID) pathway is relatively under-explored in the context of cancer biology. Recent studies [40] report a significant association between elevated circulating branched-chain amino acid (BCAA) levels and lower breast cancer risk in the NHSII cohort. This observation suggests that reduced BCAA catabolic activity may be linked to disease progression. The identification of the VLID pathway as a discriminative feature in our model could reflect underlying metabolic shifts, highlighting its potential as a novel biomarker for cancer prognosis.

Similarly, recent literature underscores the key role of tryptophan metabolism in breast cancer progression and immune evasion through the immunosuppressive kynurenine pathway [41], [42]. Single-cell transcriptomic analyses show macrophages in breast cancer tissues exhibiting elevated tryptophan metabolic activity, which correlates with M1 type polarisation and cytolytic CD8⁺ T cell infiltration, thereby serving as both an immunotherapy response predictor and a marker of subtype specific immune microenvironments [43]. Taken together, apart from well-studied pathways in cancer (e.g. p53 or mismatch repair), our study highlights pathways of emerging importance (i.e. metabolic and transport pathways above) that act as key drivers of tumour biology. These pathways play crucial roles in resistance mechanisms and immunosuppressive networks in breast cancer, making them promising targets for future research and therapeutic intervention.

4.2 Survival Analysis↩︎

For each of above-mentioned pathways (SM, ABC-T, VLID, GBCS, and ELM), 10 genes were identified by highest Absolute Neural Path Weight (ANPW) ranking [10] and statistical tests for survival analyses on all 5 pathways and these genes were performed, selecting those with both a significant logrank separation and concordant survival directionality between low and high expression groups at the cutoff date (See Supplementary Material Section 11.2 for details). Fig. 3 shows the Kaplan-Meier curves for the upper and lower third percentiles of pathway activity scores or TPM/IPM expression values, along with logrank test p-values. It is noted that the VLID pathway shows clear survival stratification, particularly in the Metabric dataset.

4.3 Ablation study on the effect of dropout on pathway activity models↩︎

Previous studies in linear BINNs showed improved robustness when applying dropout, specifically in the latent layer [3], [22], [23]. Here, we extend this investigation to the case of a non-linear Visible Neural Network in the form of our PAAE model. As shown in Figure 4, the Centered Kernel Alignment (CKA) representation is stable across a range of dropout values. Additionally, the absolute Pearson correlation between individual pathway activity vectors across independently trained models increases as the dropout rate also increases, showing that applying dropout on the latent and, in our case, in the pathway activity layer has a positive impact on robustness. Interestingly, a phase transition was observed between 70% and 80% dropout, where the Pearson correlation increases sharply, suggesting a shift in model behaviour.

Further investigation revealed that, although higher dropout resulted in models becoming more repeateable, representational power declined, especially beyond the previously mentioned phase transition, as seen in Figure 12 for ROC AUC and Concordance Index. Heatmaps to visualise models with different dropout rates (Figure 13), illustrate that, although increasing dropout makes features more reproducible, it unsurprisingly also “bleaches” the signal from pathways, generating pathways that are not only more repeatable across different training runs, but also more similar among themselves. This likely reflects a higher noise-to-signal ratio during the training process, which drives the model toward more generic representations.

Figures 14 and 15 suggest that representations generated beyond the phase transition are significantly different from those at lower dropout. Models with 0-70% dropout exhibit comparable similarity to those at 50% dropout, while being markedly different from 90% dropout. Additionally, the average signal each individual pathway has with respect to survival and classification changes depending on the amount of dropout used especially on the 80% and 90% dropout models (see Figures 17 and 16). There is an especially marked difference with respect to the Mutual Information with the PAM50 subtypes. When analysing the ranking agreement across models, higher dropout seems to reduce consistency in mutual information (Figure 19), while agreement in survival ranking remains stable, with a marked increase for the highest dropout value (Figure 18).

4.4 Interpretable Subtype Classification using the Pathway Activity Space: A Case Study↩︎

Model interpretability tools can assist medical practitioners in understanding and contextualising predictions. However, inherently interpretable models, such as BINNs, enable the generation of direct visualisations owing to the internal representations generated by the model itself. In this subsection we report a case study demonstrating a fully-interpretable subtype prediction model. Importantly, the objective here is not to optimise predictive accuracy, but to illustrate the interpretability of our model. For this analysis, we use a KEGG-based PAAE model, but similar interpretability techniques could also be applied to other models. The top most informative pathways with respect to Mutual Information for each subtype in the TCGA dataset are selected, and our model is subset to include only these pathways. A decision tree classifier is trained with the maximum number of nodes equal to the number of subtypes. To ensure that all classes are represented in the leaf nodes, SMOTE oversampling was used [30]. The resulting decision tree is shown in Figure 5.

This approach yields a clear, human-interpretable model that separates samples into each of the five subtypes. Furthermore, by adapting the UMAP space to use only the 5 top pathways rather than the full pathway set, a smoother UMAP featuremap in terms of pathway activities that drive these subtypes, is generated. In Figure 21 (top four plots) illustrates UMAP embeddings generated for all KEGG pathways, as well as for the top 5 pathways. In Figure 6 two sets of six plots are shown, analogous to those in Figure 2, showing featuremaps and KDE contours showing the area over the UMAP space. These indicate regions where samples that satisfy the decision tree criteria for classification as Basal in both the TCGA (middle) and Metabric (bottom) datasets. As expected, a model that is constrained to this level of interpretability performs worse than full models, as can be seen from confusion matrices and receiver operating characteristic curves in Figure 20. While this case study demonstrates the potential for direct interpretability in our approach, it also highlights the complexity of cancer and indicates that accurate and robust prediction requires the integration of a larger set of interacting pathways.

4.5 Multi-Omics for PAM50 Classification and Survival Prediction↩︎

This subsection details our multi-omics integration analysis using the PAAE model. Our primary objectives are to demonstrate how multi-omics data can be integrated with PAAE, to assess the general benefits of using multi-omics approaches over single-omics, and to quantify the individual contribution of each omics layer to the integrated model’s performance. We illustrate the effectiveness of early versus late integration strategies, as defined in Subsection 2.4, for multimodal data in the context of PAM50 subtype classification and survival prediction, and determine the effect of specific modalities in data integration using the PAAE method. As shown in Figure 7, early integration rarely outperformed either of the late integration methods for the classification of PAM50, achieving superior performance in only 6 of 63 cases. In most cases, the inclusion of mutation data led to a marked decrease in performance, suggesting that this layer may interact antagonistically with others in the early integration framework. It should be noted that this adverse effect was not observed in the two late integration tests.

The highest median ROC AUC was achieved using late-mean integration of all omics layers except methylation. This was closely followed by late-concat integration of all layers except mutation and of all layers combined. Although other individual layers generally performed worse than gene expression alone, the multi-omics integration approach yielded synergistic improvements in performance, particularly through late-integration approaches.

We also evaluated additional metrics, including classification accuracy (Figure 22), and performed unsupervised clustering using the K-Means algorithm, assessing performance with mutual information (Figure 23) and Rand Index (in Figure 24). The overall conclusion, i.e. that multi-Omics integration was beneficial, was consistent across these metrics. Although the best-performing methods varied, the best performing models were observed in integrating 4-6 omics layers.

For survival prediction our results also indicate that multi-omics integration can enhance performance (Figure 8). However, the late-concat approach was less consistently effective in this setting, and the negative contribution of the mutation and CNV layers was more pronounced, particularly given that these layers performed worse than random chance.

4.6 Omics Layer Contributions↩︎

We then analysed the marginal contribution of each omic layer on a previously existing combination, as described in Subsection 2.5, and shown in Figures 9 and 25, to investigate which omics layer had the greatest individual impact on each task. The mutation layer consistently exhibited a strong negative impact on PAM50 classification (ROC AUC) when using early integration (Figure 9). This effect, however, was largely mitigated with late-concat integration, likely because this method allows the classifier to effectively ignore less informative layers. Late-mean integration also reduced the negative influence of the mutation layer. As expected, the gene expression layer was the only one to consistently improve PAM50 classification, likely because the PAM50 molecular subtypes were based on gene expression profiles.

The \(\mu\)RNA and RPPA layers showed the next highest median performance increase, with the RPPA layer exhibiting a positive median effect only in late integration. These layers consistently contributed positively only when late integration was used. We hypothesise that this is due to these integration methods, especially late-concat, allowing noisy pathways within these layers to be ignored or mitigates the effect of these pathways.

As shown in Figure 25, for survival prognostication, the late-mean integration method consistently outperformed all others, and was the only integration method that allowed for an consistently positive impact across all omics layers. The most beneficial layers, in order of median impact, were RPPA, gene expression and \(\mu\)RNA.

Interestingly, late-concat integration demonstrated worse performance than late-mean integration in this context. While this may appear to contradict our previous observation that late-concat can mitigate the impact of noisy features, further analysis suggests that one of the issues with survival information also comes from how noisy labels were, and how often features seemed to dictate entirely the survival status of a patient, which degraded, and sometimes impeded, convergence.

5 Discussion↩︎

In this work, we further evaluate our PAAE [18] deep learning framework that incorporates prior biological knowledge, with a view to improving patient stratification and prognostic prediction via analysis of multi-omics data, repeatability studies, and by providing changes to our validation. Our results reinforce that our pathway activity autoencoder frameworks, although focusing on unsupervised dimensionality reduction, translate into better latent representations, while also providing an intermediate latent representation that is directly translatable into a concept relevant to medical practitioners.

Our analysis of the effect of dropout on our BINN model (PAAE) confirms findings by [22] and [23], showing increased correlation between outputs of independently trained models with higher dropout rates. However, our results further demonstrate that excessive dropout can lead to a collapse in model performance, particularly in non-linear BINNs. In other words, although representations become more similar, they exhibit reduced predictive power for both downstream tasks. Furthermore, we observe that dropout not only affects overall performance, but also alters which pathways contribute in predictive tasks, thereby caution in the use of this technique is advisable. We also report an extreme example of interpretability with our model, where we select the pathway that has the most OvR mutual information per subtype to build an extremely simplified model to aid clinical interpretability. We show how to build a decision tree that classifies samples with minimal redundancy in its leaves, and analyse performance, while also providing visualisation tools to aid clinicians in deriving specific pathway combinations that contribute to distinct patient subtypes.

Importantly, our framework works in multi-omic data analysis, and provides key insights into contributions of specific modalities in subtype classification. Our study showed that, although Gene Expression, RPPA and \(\mu\)RNA layers had the most positive marginal contribution overall, for none of the tasks, a combination composed of exclusively these three layers was the best, and more often than not adding other layers acted synergistically. Furthermore, there seemed to be a disconnect between PAM50 classification and Survival Regression, with RPPA being a better indicator for the latter, while Gene Expression more closely aligned with PAM50. However, we believe that the results reflect the particular characteristics of our model and experimental setting, and should not be extrapolated without further validation.

Overall, our study has achieved

.

As is common in complex modelling studies, our approach comes with certain limitations. We note the non-convex nature of deep neural network optimisation, which is a well-recognised limitation in the field and an area of on-going research. Although we have thoroughly studied the repeatability of our model and identified both artifact-specific and globally-significant pathways, it is important to acknowledge that deep learning models trained on relatively small datasets cannot be expected to predict the full complexity of biological relationships underlying a disease. Our study on latent space dropout, even though echoing results from previous work in terms of model repeatability [3], [22], [23], also highlights limitations of such regularisation techniques. This observation opens up new avenues for investigating more principled approaches to improving model robustness and biological fidelity. Another promising direction for future work is the training of omics foundation models (e.g.[44]–[47]), which leverage larger datasets to enhance model stability and generalisability. However, these models are mostly trained on single-cell data, and their applicability to bulk modalities remains an open question. Similarly, tabular foundation models (e.g.TabPFN [48], [49]) mark recent progress, but their deployment in high-dimensional omics context may face scalability challenges due to increased memory requirements and the comparatively lower dropout rates characteristic of bulk data.

Acknowledgements↩︎

We thank Roman Laddach for help with data acquisition, and Roman Laddach and Wai Yee Wong for fruitful discussions. MW and ST acknowledge funding from King’s College London and the A*STAR Research Attachment Programme (ARAP) to PHdCA. MW acknowledges funding by the AI, Analytics and Informatics (AI3) Horizontal Technology Programme Office (HTPO) seed grant (grant no: C211118015) from A*STAR, Singapore. ST acknowledges funding from the British Skin Foundation (006/R/22) and the UK Royal Society (IES\R2\222084).

Supplementary Material for Incorporating Prior Knowledge in Deep Learning Models via Pathway Activity Autoencoders↩︎

6 Alternatives for Omics-Pathway Mappings↩︎

We briefly considered other alternative databases for non-gene-based omics layers. While the dataset provided by [50] uses a “merged” Methylation layer, we considered whether to use Illumina’s 27K or 450 CpG-to-gene mapping. Since all annotations in 27k were also present in the 450 annotation, we decided to proceed with the more up-to-date annotation file, especially since the 450 there are multiple genes for the same CpG site.

For \(\mu\)RNA, we briefly considered using predicted target genes instead of the experimentally validated genes of miRTarBase2025 miRTarBase_SE_WR [31]. That is, we considered using information from either miRDB [51], [52], TargetScanHuman [53], or miRNet [54]–[56]. However, we noticed that when using predicted targets the number of target genes for a significant portion of \(\mu\)RNA measurements in our dataset grew significantly, and compounded over each gene possibly being into more than one pathway. This resulted in every \(\mu\)RNA being associated, on average, with 2.3 pathways when following [52]’s recommendations. This translated in the memory requirement for training an early integration model almost doubling when adding the \(\mu\)RNA to a model with only Gene Expression, CNV, Mutation, and RPPA, despite the number of unique features in the \(\mu\)RNA layer being less than 5% of the total unique features.

7 Autoencoder Architecture and Training Environment For RNA-Seq Experiments↩︎

We train the AE, VAE, PAAE, and PAVAE models using batch training for 1024 epochs with a learning rate of \(10^{-4}\) using the Adam optimiser under the pytorch v1.12.0 and the skorch v0.11.0 frameworks, all other components of our pipeline are built using scikit-learn v1.0.2. These models are trained with a dropout rate of 50%.

These models used rectified linear units as nonlinearities after each hidden layer except the last layer in each module, i.e., the last layer in the pathway encoders, and the latent space encoder and decoder are all simple linear layers without any nonlinearity. We do this to allow our model to output any real value, and, due to our network initialization method, to produce roughly normal outputs.

For the Variational models we use the loss scheduling/warm-up procedure described in the main paper, using \(T_{s}=32\) and with annealing duration \(128\) for \(T_{e} = T_{s} + 128\) in our annealed scheduling function \(S_{T_{s},T_{e}}(t)\). We use a gaussian distribution defined by \(\mu\) and \(\sigma\), from which the values of \(z\) are sampled using the reparametrization trick to maintain differentiability [57].

For our baselines, we either used the default values when it was available (PASL) or performed grid search within a certain budget. For PACL we trained both with the original sizes, adapting the last layer to have a dimensionality of 64, as our models’ latent space also had 64, while also adding layer sizes that mirrored our choices for the AE model. Since we could not fit the entire PACL model in memory given our setup, we limited the batch size to 64. For PathME we used a hidden size that mirrored our choice for the PAAE model’s pathway activity sizes, while performing grid search on a smaller subset of the hyperparameter grid proposed by the authors. We had to perform a search on a smaller space since we performed a more exhaustive grid search instead of employing bayesian optimisation, which allowed us to explore more combinations (54) than in the original paper (50), while also removing redundant searches (such as searching with a dropout of 100%). Due to a software incompatibility with the CUDA libraries in PathME’s library, we had to run it using the CPU, and thus we limited the batch-size to the highest batch size tested in their model (32).

For all pipelines, we first perform an internal 8-fold cross validation for hyperparameter tuning, selecting the model with highest ROC AUC from a grid search in the parameters shown in Tab. 1, where the classifier was applied to the hidden representation vector (i.e., \(z\) or \(\mu\)), which resulted in the following hyperparameters being chosen:

• Autoencoder: Encoder with 2 layers with 128, and 64 units; using a Logistic Regression classifier.

• Pathway Activity Autoencoders: Pathway Activity Encoders with 2 layers with a 32-unit hidden layer; an Encoder with only a single 64-unit layer; and using a Logistic Regression classifier.

• Variational Autoencoder: Encoder with 2 layers with 128, and \(2\times64\) units (64 for \(\mu\) and 64 for \(\sigma\)); using a Logistic Regression classifier.

• Pathway Activity Variational Autoencoders: A linear Pathway Activity Encoder with only a single linear layer; Encoder with 2 layers with 128, and \(2\times64\) units (64 for \(\mu\) and 64 for \(\sigma\)); and using a Support Vector Machine classifier.

• PACL with \(k=1\), and 1024 epochs. For the Hallmark Genes pathway set the best model had two hidden layers with 128 and 64 units and a Logistic Regression classifier was used. For the KEGG pathway set the best model had three hidden layers with 500, 200 and 64 units and a Support Vector Machine classifier was used.

• PathME with 16 hidden units, \(\alpha=1\), \(\lambda=10^{-9}\) and \(0.001\) with a Logistic Regression Classifier.

We train the Logistic Regression models with L2 regularization with \(C=1\), for a maximum of \(100\) iterations using the lbfgs solver, and using the softmax function for the multi-class cases. For the Support Vector Machine classifiers we use \(C=1\) as the l2 regularization parameter and a Radial Basis Function kernel. For the Random Forest Classifier we build 100 Trees using Gini impurity, considering \(\sqrt{d_{input}}\) features for splitting, boostrapping samples when building trees.

After the hyperparameters were selected, we train the models in a completely unsupervised fashion on our training dataset (TCGA) and then use the models to compress the input to their latent representation, using \(z\) for the AE and PAAE models and the posterior means \(\mu\) for the VAE and PAVAE models, as well as a separate compression of our input to its pathway activity score representation \(a\) for Pathway Activity models, then using this compressed latent representation as the input for training the classifier.

After training, for the AE, VAE, PAAE, and PAVAE models, we then normalize the input of our test dataset using the same method as was used for the training dataset, since our test dataset is on a completely different scale and was produced by a different gene expresssion measuring technology. We then use this as input for the all models to compress the features, and use the compressed representation as input for the classifiers that produce predictions of each sample’s cancer subtypes.

We adapted the original code provided by PASL and PathME to run their models, and use a reimplementation of PACL to allow us to use the probability values instead of binary vectors as the internal representations, which allowed for better performance than the original code. It is of note that for PASL we used standard normalisation (z-scoring) since that is what was used internally in the model, and when it came to applying the model to an external dataset we z-scored the dataset and then multiplied it with the inverted model matrix, instead of using the models’ learned mean and average. For PACL, since RBMs need a probability input, we z-scored the values as described in the original paper, but then performed 0-1 scaling so that each value could be treated as a probability for the training sampler. We could not run the PathME model on our GPU since the library it was developed on was incompatible with out setup. Furthermore, PathME had stability issues and produced nan or inf values for some pathways, whenever this happens we drop these pathways during the internal and external cross validation.

For all tests we consider the macro-averaged Area Under the Receiver Operating Characteristic Curve (ROC AUC) as the main metric and, unless otherwise stated, we apply the Wilcoxon Rank-Sum test for testing statistical significance between different medians. We also provide other metrics such as accuracy, precision, recall and F1 score. Note that our Pathway Activity methods can’t access the full set of input genes when building their internal representation, while at the same time having to rebuild the entirety of the input.

7.1 Model sizes and alternatives↩︎

Although nominally some of our baselines seem to have smaller parameter counts than our models, that is mostly due to the fact that these models either (1) only reconstruct the input partially, which is what is done by PathME, or (1) don’t reconstruct the input, i.e. PACL and PASL. One can see in Tables 2, 3 and 4 the model parameter sizes for every hyperparameter configuration. In Tables 3 and 4, we show a comparison with other models when we consider

For brevity’s sake, we discuss here only the results for Table 3, since those in Table 4 are equivalent. One can see that if we were to reconstruct the full gene matrix using only pathway information [lst:sz-cmp-full], our PAAE model is vastly more efficient than using a PathME model, since the PathME model builds one separate autoencoder per pathway. If we were to use our PAAE model to only reconstruct genes that are contained within the pathway set [lst:sz-cmp-pway], then our PAAE model continues to be competitive w.r.t. to PathME, and we can still keep the performance from reconstructing the full gene matrix [lst:sz-cmp-full] as we’d only need to prune the decoder, only needing these extra parameters for training. If we only need the encoder [lst:sz-cmp-enc], while our more complex models requires more parameters than PACL, they still vastly outperform it. Furthermore, our model requires less parameters than PASL, even though PASL’s parameter count is pruned to only include nonzero values. Again, our model achieves this while still being able to keep the performance from reconstructing the full gene matrix [lst:sz-cmp-full] since we’d only need to prune our model and discard these parameters after training.

7.2 PACL and PASL external validation discussion↩︎

Both PACL and PASL experienced catastrophic failures when translating to another dataset. While PACL’s failure on its latent space was expected due to the information loss caused by using locally-trained RBMs, its failure on the pathway space was unexpected. The same is also true for the PASL results. We posit that this is due to the normalisation methods used by these methods. Since PACL used z-scoring and since we also needed to 0-1 normalise the input, we chose to keep z-scoring and min-max scaling as its default scaling method. PASL similarly uses z-scoring, which we attempted to use natively using its own codebase, but which caused numerical errors and made the model fail altogether. We then proceeded to apply the model in Python instead of Matlab, which avoided the numerical errors but still produced subpar results. We left a comparison of improved versions of these methods’ pipelines to future work due to computational budget constraints.

8 Interpreting the Pathway Activity Scores↩︎

As we can see in the original PAAE description [18], the pathway activity space that we visualise is built of the individual pathway activity scores of each pathway. Each of the pathway activity scores is built by stacking fully connected layers using the normalised gene expression of the genes belonging to a specific pathway as an input. The encoder in \(E_{p_{j}}\), if it has \(k\) layers, can be described as in Eq. 11 :

\[\label{eq:paae-activity-encoder} \begin{align} a_{j} = E_{p_{j}}(x) &= h_{j,k} \\ h_{j,k} &= W_{j,k} \times h_{j,k-1} + b_{j,k} \\ h_{j,i} &= f_{j,k}(W_{j,i} \times h_{j,i-1} + b_{j,i}) \\ h_{j,1} &= f_{j,k}(W_{j,1} \times x_{:,p_{j}} + b_{j,1}) \end{align}\tag{11}\]

Where \(W_{j,i}\) is a \(\Re^{d_{j,i-1} \times d_{j,i}}\) matrix, and \(f_{j,k}\) is a nonlinear function, not present in the last layer so that our output is roughly normal-shaped.

Since the weights \(W_{j,i}\) can be both positive and negative, one should not expect the final result \(a_{j}\) to have a positive/negative directionality w.r.t. the pathway’s input \(x_{:,p_{j}}\), and one should always consider the pathway activity scores in terms of contrasts between samples. If there is a biologically known directionality, then that can be used to either confirm or align the directionality of the pathway activity score by inverting its sign.

The pathway activity score \(a_{j}\) is the value that we use in both our clustermap (Figs. 11, 10) and featuremap (Fig. 2) visualisations, and thus both these visualisations should be interpreted with the caveat above. Furthermore, in cases where nonlinear layers are involved in the calculation of the pathway activity score \(a_{j}\), one should not expect the values to behave linearly.

9 Clustermaps↩︎

For clustermaps, in Figs. 11 and 10, we cluster along both rows and columns, using the cosine distance metric, as it is frequently used for deep learning embeddings, and we use seaborn’s clustermap function.

10 Mutual Information↩︎

We calculate non-OvR Mutual Information using scikit-learn’s mutual_info_classif function, inside the feature_selection module, to measure the mutual information between a continuous variable and the discrete cancer subtype labels. We do not filter for mutual information amongst variables, so multiple features might have a high MI w.r.t. to one class, instead of having multiple features with each having high MI w.r.t. to different classes. This differs from the OvR mutual information used in other parts of the paper.

11 Survival Analyses↩︎

For survival plots and statistical tests (Fig. 3) a time window limit of 5 years (\(365\times5=1825\) days) was set. We use the lifelines library for both the tests and the plots. All of the analyses are done in an almost-unsupervised fashion, being that the only supervision we use is that we choose the top 5 pathways w.r.t. mutual information with the PAM50 subtypes and, since some subtypes might have different survival outcomes, this might influence the survival curves’ relationship with survival.

For survival curves we select the most important genes from each pathway (as per ANPW, see 11.2). For each of these genes, we divide the samples into “high” and “low” expression, which are defined as the upper and lower third quantiles of the TPM/IPM values (See  11.1 for details). We then perform logrank tests and select those genes that have a significant logrank p-value on the training dataset. Then we further filter to those genes which also have both a significant logrank p-value on the test dataset, and matching low-high survival sign.

For the low-high survival sign, we measure the surival rate at the end of the time window limit. If the “high” expression group has higher survival rate than the “low” expression group in the training dataset, then it has a matching survival sign if it also has higher survival rate in the test dataset. The same logic follows that if the “high” expression group has lower survival rate than the “low” expression group in the training dataset, then it has a matching survival sign if it also has lower survival rate in the test dataset.

11.1 Preprocessing for Gene Expression Comparisons↩︎

For the analyses involving raw input features, such as those in Fig. 3, we re-normalise the genes from \(\log(FPKM+1)\) to \(\log(TPM+1)\) in the TCGA dataset, and from \(\log(Intensity+1)\) to \(\log(IPM+1)\), where FPKM stands for Fragments Per Kilobase Million, TPM for Transcripts Per Million, Intensity is the raw intensity value measured in the Metabric dataset, and IPM is “Intensity Per Million”, where the intensity is normalised per sample as if it was a raw count value with \(IPM_{g} = \frac{I_{g}\times10^{6}}{\sum_{g}I_{g}}\).

11.2 Feature Importance↩︎

Furthermore, we can get the most important features for each pathway using either neural path weights (NPW, Eq. 12 ) or absolute neural path weight (ANPW, Eq. 13 ) [10].

\[\label{eq:npw} \operatorname{NPW}_{j} = \prod_{i=1}^{k} W_{j,i} ~,~ \operatorname{NPW}_{j} \in \Re^{|p_{j}|}\tag{12}\]

\[\label{eq:anpw} \operatorname{ANPW}_{j}[g] = |\operatorname{NPW}_{j}[g]|\tag{13}\]

Note, however, that since our models might compute the pathway activity score value through a combination of nonlinear functions, these values should not be interpreted as a linear influence on the value, and combinations with other values may suppress or excite the score more than simply increasing/decreasing the genes with highest ANPW values.

Also note that one can do this interpretability technique for any point in any neural network, but we only define this for our pathway activity score for simplicity sake.

12 Supporting Figures↩︎

References↩︎

---

1. This helps avoid features which completely determine the survival status over the cross validation runs↩︎