arXiv CSopen access

Biologically Informed Deep Neural Networks for Multi-Omic Integration, Pathway Activity Inference and Risk Stratification in Cancer

arXiv:2607.05306v1 [cs.LG] 6 Jul 2026

Biologically Informed Deep Neural Networks for Multi-Omic Integration, Pathway Activity Inference and Risk Stratification in Cancer Pedro Henrique da Costa Avelar kcl,astar,uom , Le Ou-Yang Tsoka kcl

szu

, Min Wu astar , and Sophia

kcl

Department of Informatics, King’s College London, London WC2B 4BG, United Kingdom astar Institute for Infocomm Research, Agency for Science, Technology and Research (A*STAR), Singapore 138632, Singapore szu College of Electronics and Information Engineering, Shenzhen University, Shenzhen 518054, Guangdong, China uom Current Affiliation: Division of Informatics, Imaging and Data Science, Faculty of Biology, Medicine and Health, University of Manchester 2025

Abstract

model robustness and consistency, excessive regularisation can reduce predictive performance. Finally, viIntegrating complex, multi-omics data presents sig- sualizations of the learned feature space illustrate the nificant challenges. Existing approaches often face framework’s intrinsic transparency and clinical relea trade-off between model interpretability and repre- vance. The results underscore the value of multi-omic sentational capacity, with most either relying on post- integration and delineate the impact of individual hoc interpretation or use linear models that may over- omics layers, establishing practical guidelines for intelook complex interactions. We report Pathway Activ- gration within our framework. Overall, our pathway ity Autoencoders for the multi-omics setting, which activity autoencoder frameworks yield superior latent embed prior knowledge via pathway-informed archi- representations that are biologically meaningful and tectural constraints, fostering interpretability, while are directly translatable into clinically relevant inpreserving representational power. Our multi-omic sights. Code related to this work will be made freely framework is applied in the context of breast can- available at github.com/phcavelar/pathwayae Multi-Omic Integration; Autoencer and is evaluated in survival prediction and sub- Keywords: coder; Biologically-Informed Neural Networks; Intertype classification with results indicating a positive pretable Deep Learning; Breast Cancer effect of integration. We conduct analysis of individual omics layer impact on end-task performance, revealing that gene, protein, and microRNA expression layers provide the strongest contribution. Repeatability studies indicate that, while dropout improves 1

Unseen genes

reconstruction

Copy Number Variation DNA Mutation DNA Methylation mRNA Expression microRNA Expression RPPA

subset

PAAE

Unsupervised Representation Learning

concatenate

Multi-Omics Integration

external classifier

mutual information

absolute neural path weights

2

0

2

4

6

8

10

12

PAM50 Normal LumA LumB Basal Her2

0.2

0.1

0.0

0.1

0.2

0.3

gex

meth

mutation rppa omics_layer

cnv

mirna

0

2

4

10

LumB

Her2

12

integration early late-concat late-mean

Omics Layer SHAP Impact Analysis

LumA

8

SM<=-0.60 True

VLID≤0.7

6

P53≤-0.32

Normal

PC≤0.07

GBKS≤-0.6

2

Prediction Conditions for Basal on Metabric

Interpretable Pathway Space Basal

2

0

2

4

6

8

Cancer Subtype Prediction, Survival Analysis, Feature Importance

Figure 1: Diagram of the PAAE model. The representation learning module (left) receives a subset of the input gene expression matrix and constrains the architecture to reflect specified pathway gene sets, each producing a single pathway activity value a. These activities are further transformed into an internal latent representation z, which is then used to reconstruct the entire gene set. Results from downstream tasks are shown on the right. Heatmaps correspond to the best-performing configuration identified for the “Hallmark Genes” pathway set. Matrices are drawn to scale and subset for readability. The PAAE framework yields directly interpretable pathway activity scores and higher-level latent representations that can be leveraged for downstream tasks such as classification, clustering, and survival analysis. Each pathway activity is independently interpretable, and its features can be inspected for fine-grained biological insights.

subset

subset

subset

subset

accuracy_effect

2

1

Introduction

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: (i) An update of the methodology for identifying pathways relevant to subtype classification by using OvR Mutual Information, ensuring that even rare subtypes are represented in the selected pathways; (ii) Thorough repeatability and robustness studies, with 2560 models having been evaluated in total showing that while dropout, a common form of regularisation used to promote repeatability, does improve correlation between models, it might also have a negative impact on performance; (iii) A fully automated pipeline to extract a completely interpretable model is provided, resulting in a decision tree with the same number of leaf nodes as there are subtypes; (iv) Prior-knowledge-enhanced multi-omics integration is demonstrated to be beneficial to using single-layers; and (v) A thorough analysis of the impact of individual omics layers in end-tasks for the TCGA Breast Cancer dataset is performed through an adapted Shapley Value analysis.

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 [Bode and Dong, 2017, Chen et al., 2021]. However, the integration of highdimensional, heterogeneous molecular data, such as genomics, transcriptomics, and proteomics, pose significant challenges, particularly in balancing predictive accuracy with interpretability in machine learning models [Selby et al., 2024, 2025]. While deep learning has demonstrated remarkable success in drug response prediction [Xie et al., 2024, Abbasi et al., 2024], phenotypic classification [Van Hilten et al., 2021, Yap et al., 2021], disease subtype classification [Chaudhary et al., 2018, Uyar et al., 2021], and survival prognosis [Ching et al., 2018], its translational impact is often limited by the inherent opacity of its black-box representations [Selby et al., 2024, 2025]. Interpretability is critical in medical applications, where decision-making must be transparent and biologically plausible [Kaur et al., 2022]. Post-hoc explainability methods (e.g. [Lundberg and Lee, 2017], [Ribeiro et al., 2016]) explain individual features retrospectively, but do not constrain the internal factors in the model and often require re-analysis for each prediction or dataset [Li et al., 2022, Dwivedi et al., 2023, Abdrakhimov et al., 2024], 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 (Selby et al. [2024, 2025]) 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 [da Costa Avelar et al., 2024] fall inside the class of BINNs that use

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 [da Costa Avelar et al., 2024], as illustrated in Figure 1. 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 individuallayer 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 [da Costa Avelar et al., 2023], shown in Figure 1, is a Biologically-Informed Neural Network [BINN, Selby et al., 2025] framework based on the Autoencoder [AE, Hinton and Salakhut3

dinov, 2006] and Variational Autoencoder (VAE) deep learning architectures. It works by incorporating pathway information in its architecture, building one pathway encoder block Ep for each pathway p in a pathway set P . Each of these pathway encoders receive as input xp only the molecular features which map to a gene that belongs inside the pathway p, and then output a single pathway activity score ap = Ep (Xp ). Thus, we generate a pathway activity space a by concatenating ∥ all individual pathway activities, as in Equation 1 a = ∥p∈P ap = ∥p∈P Ep (xp )

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 [Cox, 1972] as a ranking metric for survival relevance, analogous to how MI is used for classification tasks. For further details, see Supplementary Material Sec. F.

(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 x̂ = D(a), which is then used to optimise the network parameters for all Ep , Ω, and for D, Γ using the a distance function d. In practice, we use 2.3 Repeatability and Robustness the mean squared error as d, since we normalise our Analyses input, giving us the optimisation objective in Equation 2. For robustness and repeatability analyses, we measure the Pearson Correlation between neurons across 2 min (x̂ − x) (2) different repeated random initialisation and trainΩ,Γ ing runs, following work reported in Fortelny and Bock [2020] and Seninge et al. [2021]. We expand on this work by incorporating other correlation metrics, 2.2 Downstream Tasks such as Kendall’s τ -b coefficient Kendall [1945] and For classification, we used the best-performing classi- Centered Kernel Alignment (CKA, Cristianini et al. fier found during internal validation on the RNAseq [2001], Kornblith et al. [2019]) similarity as a meastudy (i.e. Logistic Regression LR), using as input sure of overall repeatability, by analysing both the either the latent space z or pathway activity space aggregate ranking of downstream-task metrics and a, to predict the cancer subtype. Model quality is the ranking agreement between models through Avassessed with usual classification metrics (Accuracy, erage Precision K (AP @K [Manning et al., 2008], Precision, Recall, and the Area Under the Receiver and by performing experiments building aligned conOperating Characteristic Curve, ROC AUC), macro- sensus models from multiple repeats. averaged and, in the case of ROC AUC, one-vs-rest (OvR) was used. Cancer subtype identification is 2.3.1 Aligned Consensus Models also performed in an unsupervised fashion through clustering and clustering performance is evaluated us- When performing repeatability and robustness analying Mutual Information (MI) between the generated ses, we calculated Aligned Consensus Models (ACM) clusters and the target classes. Additionally, MI is in order to study whether they were more repeatable 4

than our non-consensus models. To build an ACM, one takes the pathway activity vector aj of a repeat j ∈ [1, . . . , k] out of k repeats. We then take the first model a1 as a reference, and calculate the aligned average pathway activity vector ai for pathway i as in Equation 3. ai =

1X sgn(corr(aji , a1i ))aji k

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 xo , 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, Ep , then specifically focuses on the molecular features linked to the genes (Gp ) that are relevant to its particular pathway (p ∈ P ).

(3)

j≤k

That is, we take the average of the pathway activity vector aji for pathway i and repeat j, multiplied by the sign sgn(c1,j i ) of a correlation metric corr between it and the reference model’s pathway activity vector c1,j i = corr(aji , a1i ). In practice, we take the Pearson Correlation Coefficient for corr. 2.3.2

x = ∥o∈O xo

Ranking Agreement

PAAE/PAVAE is optimised to reconstruct the unified profile x, enabling each pathway activity score ap 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 ao (Equation 6). If an omics layer lacks genes for pathway p, we exclude it from that model, denoted as p ̸∈ Ao .

To compute the ranking agreement between two individual repeats j, m ∈ [1, k], we take a metric f (a, y) that takes a pathway activity vector a and a target y, and computes a score si for each pathway i w.r.t. to that metric. We then take the two rankings Rj and Rm , where Rj [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 [Manning et al., 2008], modified to only calculate relevant items up to K as shown in Equation 4.

AP @Kj,m =

(5)

aop ≈ Eop (x), Eop (x) = (Eop 1 · . . . Eop k )(xo )

1 X P (j) · 1(Rj [r] ∈ Rm [: K]) (4) K

(6)

r<K

After generating pathway activities from each model, we integrate using two methods: late-mean Where P (j) is the precision at rank j, 1 is an indiintegration, and late-concat integration. Late-mean cator function for whether the item is in the relevant integration takes the arithmetic average of the pathitems, and Rm [: K] are all items from the beginning way activities of each model to yield an integrated up to K in ranking Rm . Note that, we assume that pathway activity vector a (Equation 7): K ≤ |P |. X 1 ap = P ao (7) 2.4 Early and Late Multi-Omics Patho∈O,p∈Ao 1 o∈O,p∈Ao

way Activity Inference and Analysis

For late-concat integration, individual pathway activity vectors from each omics layer are concatenated, Multi-omics integration is categorised into sequen- producing an integrated pathway activity vector a tial, late and early (joint) integration methods [Uyar (Equation 8, where ∥ denotes concatenation along et al., 2021]. Late integration methods analyse each features). 5

a = ∥o∈O ao

2.5

3

Experimental Setup

3.1

Repeatability Studies

(8)

Individual Omic Layer Contributions

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 singleomics 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 val −1 n−1 1 X ues of P r(dropout(a)) ∈ {0, 0.1, 0.2, . . . 0.9}, with the (v(S∪{i})−v(S)) (9) φi (v) = highest dropout level resulting in an average of ap|S| n S⊆N \{i} proximately five active pathways per training epoch. For each dropout level, 16 repeats per aligned conWe then adapt this idea to the concept of a sensus model were conducted and 16 consensus modmulti-view (i.e., multi-omics) machine learning probels were built, yielding k = 256 models per dropout lem, similarly to that in Shapley regression values value, and thus 2560 models in total. [Lipovetsky and Conklin, 2001]. Given a set of views (omics layers) O and a performance metric v, we calWe employ scikit-learn v1.0.2 culate the average marginal contribution a view o as a for mutual information (MI), using function f of the difference between the contribution feature selection.mutual info classif to of every subset of features S such that S ⊆ O \ {o}, estimate the mutual information between as in Equation 10. 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, φi (v) = f ({(v(S ∪ {o}) − v(S))∀S ⊆ O \ {i}, S ̸= ∅}) specifically CoxPHFitter.concordance index to (10) compute the Concordance Index (CI) of individIn our specific case, we use f as the median value ual pathways or models. Correlation metrics are of all marginal contributions, which we define using computed using scipy v1.8.1 stats.pearsonr for a cross-validated metric for v, such as ROC AUC or Pearson’s correlation coefficient, stats.kendalltau Concordance-Index, meaning we train a model based for Kendall’s τ -b Coefficient, and stats.rankdata on view combination S and S ∪ {o}. We ignore the to derive pathway rankings Rj . Ranking agreement dummy case (i.e. where there is no input, S = ∅). is assessed using an internally developed impleWe also display the distribution of the marginal con- mentation of average precision at rank AP @K tributions in plots. (Equation 4).

To evaluate the contribution of each omics layer on downstream performance, we adapt the Shapley value φi (v) Shapley [1953], which approximates a player’s average marginal contribution in a cooperative game. The definition of the Shapley value φi (v) for a player i ∈ 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 ⊆ N \ {i} participate, and subtracting it from the value if player i had joined the game, as is shown in Equation 9 below.

6

3.2

Interpretable Pathway Activity- combination of the 6 omics layers available (63 combinations in total). External validation is not perBased Subtype Classification formed 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 excluded1 , and apply an l2 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 siteto-gene mapping provided by the manufacturer for the Illumina HumanMethylation450 BeadChip. For the µRNA layer, we used the miRTarBase2025 miRTarBase SE WR database of µRNA target genes which includes target gene interactions supported by strong experimental evidence [Cui et al., 2025]. Alternative strategies were considered and described in Supplementary Material A.

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) [Chawla et al., 2002] 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.

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 C). As shown in Fig. 2, 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 3.3 Multi-Omics Analyses of the learned representations. The MO-TCGA-BRCA dataset is employed for multiUsing the KEGG pathway set, the top 5 pathomics analysis, and a PAAE model is trained with ways with highest one-vs-rest mutual information for the same hyperparameters selected in the internal 1 This helps avoid features which completely determine the validation of single-omics analysis, with the Hallmark Genes pathway. We train models for every possible survival status over the cross validation runs 7

PAM50 DORSO_VENTRAL_AX MISMATCH_REPAIR PYRIMIDINE_METAB STEROID_HORMONE_ VASOPRESSIN_REGU GLYCOSPHINGOLIPI VALINE_LEUCINE_A GLIOMA ARRHYTHMOGENIC_R NUCLEOTIDE_EXCIS ETHER_LIPID_META PANCREATIC_CANCE ENDOMETRIAL_CANC GLYCOSAMINOGLYCA FATTY_ACID_METAB HYPERTROPHIC_CAR OOCYTE_MEIOSIS RENIN_ANGIOTENSI ECM_RECEPTOR_INT PRIMARY_IMMUNODE P53_SIGNALING_PA ONE_CARBON_POOL_ PENTOSE_PHOSPHAT TRYPTOPHAN_METAB N_GLYCAN_BIOSYNT NICOTINATE_AND_N DNA_REPLICATION PROSTATE_CANCER RNA_DEGRADATION SPHINGOLIPID_MET CHRONIC_MYELOID_ PROXIMAL_TUBULE_

PAM50Call_RNAseq HISTIDINE_METABO FRUCTOSE_AND_MAN GLYCOSPHINGOLIPI VALINE_LEUCINE_A GLIOMA LONG_TERM_DEPRES GLYCOSAMINOGLYCA AMINOACYL_TRNA_B PROXIMAL_TUBULE_ SPHINGOLIPID_MET CHRONIC_MYELOID_ GLYCEROPHOSPHOLI LEUKOCYTE_TRANSE MISMATCH_REPAIR FATTY_ACID_METAB LIMONENE_AND_PIN PANCREATIC_CANCE ETHER_LIPID_META GLYCOSPHINGOLIPI GLYCOLYSIS_GLUCO P53_SIGNALING_PA ONE_CARBON_POOL_ TRYPTOPHAN_METAB RNA_DEGRADATION EPITHELIAL_CELL_ HYPERTROPHIC_CAR PEROXISOME RENIN_ANGIOTENSI DORSO_VENTRAL_AX GLYCOSAMINOGLYCA PHENYLALANINE_ME ABC_TRANSPORTERS SampleID

(a) TCGA (train)

(b) Metabric (test)

Figure 2: Clustermap based on cosine distances between sample-level inferred pathway activity vectors in the KEGG PAAE pathway activity space (see Sec. D). Colours denote BRCA clinical phenotypes (Normal in blue, Luminal A in orange, Luminal B in green, Basal in red, Her2 in purple). For clarity, only the 32 pathways with highest mutual information with respect to the class labels are shown (see Supplementary Material Sec. E for details). Fig. 11 shows similar results for the Hallmark Genes pathway set across all 50 pathways.

8

SM

VLID 0.8 0.4 0.0 0.4 0.8

0.8 0.4 0.0 0.4 0.8

ABC-T 0.8 0.4 0.0 0.4 0.8

GBCS 0.8 0.4 0.0 0.4 0.8

ELM 0.8 0.4 0.0 0.4 0.8

PAM50 Normal LumA LumB Basal Her2

Figure 3: Subtype distribution overlaid on 2D UMAP projection of the inferred pathway activity vectors (top-left), as well as featuremaps showing the intensity of each sample inferred pathway activity scores for KEGG-based PAAE pathway activity space in the Metabric dataset. For clarity, only the 5 pathways with the highest mutual information in TCGA subtypes are shown. 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 3. Key pathways prioritised by our methodology, i.e. Sphingolipid Metabolism [Ryland et al., 2011, Corsetto et al., 2023], ABC-Transporters [Xiang et al., 2011, Muriithi et al., 2020], GBCS [Huang et al., 2021, Yen et al., 2024], and ELM pathways [Benjamin et al., 2013, Yu et al., 2021], have been documented in tumorigenesis and cancer progression,

either specifically in breast cancer [Xiang et al., 2011, Ryland et al., 2011, Yu et al., 2021, Corsetto et al., 2023, Yen et al., 2024], or more broadly across various cancer types [Ryland et al., 2011, Huang et al., 2021, Benjamin et al., 2013, Muriithi et al., 2020]. 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 Zeleznik et al. [2021] 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 [Girithar et al., 2023, Ma et al., 2025]. Single-cell transcriptomic analyses show macrophages in breast cancer tissues exhibiting 9

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 [Xue et al., 2023]. 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

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 13 for ROC AUC and Concordance Index. Heatmaps to visualise models with different dropout rates (Figure 14), 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.

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 [Uyar et al., 2021] 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 F.2 for details). Fig. 4 shows the Kaplan-Meier curves for the upper and lower third percentiles of pathway activity scores or TPM/IPM expression valFigures 15 and 16 suggest that representations genues, along with logrank test p-values. It is noted that erated beyond the phase transition are significantly the VLID pathway shows clear survival stratification, different from those at lower dropout. Models with particularly in the Metabric dataset. 0-70% dropout exhibit comparable similarity to those at 50% dropout, while being markedly different from 4.3 Ablation study on the effect of 90% dropout. Additionally, the average signal each dropout on pathway activity mod- individual pathway has with respect to survival and classification changes depending on the amount of els dropout used especially on the 80% and 90% dropout Previous studies in linear BINNs showed improved models (see Figures 18 and 17). There is an especially robustness when applying dropout, specifically in the marked difference with respect to the Mutual Inforlatent layer [Fortelny and Bock, 2020, Seninge et al., mation with the PAM50 subtypes. When analysing 2021, Selby et al., 2024]. Here, we extend this investi- the ranking agreement across models, higher dropout gation to the case of a non-linear Visible Neural Net- seems to reduce consistency in mutual information work in the form of our PAAE model. As shown in (Figure 20), while agreement in survival ranking reFigure 5, the Centered Kernel Alignment (CKA) rep- mains stable, with a marked increase for the highest resentation is stable across a range of dropout values. dropout value (Figure 19). 10

KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION TCGA (log10p = 2.0)

1.00 0.95 0.90 0.85 0.80 0.75 0.70 0.65 0.60

KEGG_SPHINGOLIPID_METABOLISM

Metabric (log10p = 10.8)

TCGA (log10p = 1.1)

1.00

low high

Metabric (log10p = 9.7)

low high

0.95

low high

0.90 0.85 0.80 0.75 0.70

low high 0

500

1000 timeline

1500

0

1000 timeline

1500

0

1000 timeline

1500

0

500

1000 timeline

1500

(a)

(b)

KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION - MCCC1

Metabric (log10p = 8.7)

low high

0.9

0.8

0.8

0.7

0.7

0.6

0.6 500

1000 timeline

1500

0

500

1000 timeline

TCGA (log10p = 3.0)

1.0

low high

0.9

0

500

KEGG_SPHINGOLIPID_METABOLISM - DEGS2 TCGA (log10p = 1.4)

1.0

0.65

500

1500

Metabric (log10p = 3.9) low high

low high 0

500

1000 timeline

(c)

1500

0

500

1000 timeline

1500

(d)

KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION - OXCT1 TCGA (log10p = 3.0)

1.0

Metabric (log10p = 1.4)

low high

low high

0.9 0.8 0.7 0.6 0

500

1000 timeline

1500

0

500

1000 timeline

1500

(e)

Figure 4: Kaplan-Meier Curves for the pathways (4a,4b) and genes (4c-4e) that showed statistically significant differences (logrank test) between the upper and lower tertiles in expression in both the TCGA and Metabric datasets, with consistent survival directionality (see Supplementary Material Sec. F for details). These 4 genes represent 57.14% of the 7 genes that were significant in the TCGA logrank test, drawn from the 10 most informative genes (by OvR mutual information with the classification target) within each of the 5 most relevant pathways. 11

Figure 5: Plot illustrating the average Centered Kernel Alignment (CKA) similarity and the average Pearson Correlation between pathways obtained from models trained with different random seeds on the MO-TCGA-BRCA dataset. Models are labelled as either Aligned Consensus Model or not. Correlation increases as dropout rate increases, with sharp increase between 70 and 80% dropout, and an accompanying decrease in CKA similarity, indicating that while pathway-specific features become more consistent across models, global representational similarity slightly declines. Shaded areas indicate variance across model runs.

4.4

Interpretable Subtype Classifica- nodes, SMOTE oversampling was used [Chawla et al., tion using the Pathway Activity 2002]. The resulting decision tree is shown in Figure 6. 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

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 22 (top four plots) illustrates UMAP embeddings generated for all KEGG pathways, as well as for the top 5 pathways. In Figure 7 two sets of six plots are shown, analogous to those in Figure 3, 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 21. While this case study demonstrates the potential for direct in-

12

GBKS <= -0.6 gini = 0.8 samples = 2165 value = [433, 433, 433, 433, 433] class = Basal PC <= 0.07 gini = 0.22 gini = 0.77 samples = 425 samples = 1740 value = [375, 30, 7, 13, 0] value = [58, 403, 426, 420, 433] class = Basal class = Normal VLID <= 0.7 gini = 0.24 gini = 0.7 samples = 468 samples = 1272 value = [28, 0, 31, 4, 405] value = [30, 403, 395, 416, 28] class = Normal class = LumB P53 <= -0.32 gini = 0.46 gini = 0.63 samples = 377 samples = 895 value = [27, 269, 12, 59, 10] value = [3, 134, 383, 357, 18] class = Her2 class = LumA gini = 0.61 samples = 397 value = [1, 88, 219, 76, 13] class = LumA

gini = 0.56 samples = 498 value = [2, 46, 164, 281, 5] class = LumB

Figure 6: Decision tree generated after training it on pathway activity vectors of the PAAE KEGG model trained on the single-omics TCGA dataset. The tree is constrained to 5 leaf nodes, and SMOTE oversampling[Chawla et al., 2002] is used to ensure that each PAM50 subtype is represented in the leaf nodes. terpretability 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 8, early integration rarely outperformed either of the late in-

tegration 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 23), and performed unsupervised clustering using the K-Means algorithm, assessing performance with mutual information (Figure 24) and Rand Index (in Figure 25). The

13

Prediction Conditions for Basal on TCGA SM<=-0.60 True

12 10 8

PAM50Call_RNAseq Normal LumA LumB 4 Basal Her2 2

ABC-T 0.8 0.4 0.0 0.4 0.8

6

12 10 8 6 4 2

VLID 0.8 0.4 0.0 0.4 0.8 2

GBCS 0.8 0.4 0.0 0.4 0.8 0

2

4

6

8

2

ELM 0.8 0.4 0.0 0.4 0.8 0

2

4

6

8

2

0

2

4

6

8

Prediction Conditions for Basal on Metabric

8

SM<=-0.60 True

6 4 PAM50 Normal LumA LumB Basal Her2

2 0 2

ABC-T 0.8 0.4 0.0 0.4 0.8

8 6 4 VLID 0.8 0.4 0.0 0.4 0.8

2 0 2

2

0

2

4

6

8

10

12

GBCS 0.8 0.4 0.0 0.4 0.8 2

0

2

4

6

8

10

12

ELM 0.8 0.4 0.0 0.4 0.8 2

0

2

4

6

8

10

12

Figure 7: 2D UMAP featuremaps for TCGA (top six) and Metabric (bottom six) datasets for pathway activity vectors derived from the top 5 pathways with highest one-vs-rest mutual information. Red KDE contours are overlaid on the SM pathway featuremap in both datasets, highlighting regions that satisfy the decisions rules for the Basal subtype from decision tree depicted in Figure 6. For clarity, the same contours are displayed in black on the Class featuremaps.

14

gex mirna rppa cnv meth mutation rppa,mirna rppa,cnv gex,rppa gex,mirna gex,mutation meth,rppa gex,cnv meth,mirna gex,meth cnv,mirna mutation,mirna mutation,rppa meth,cnv meth,mutation mutation,cnv gex,rppa,mirna meth,rppa,mirna mutation,rppa,mirna rppa,cnv,mirna gex,rppa,cnv gex,mutation,rppa mutation,rppa,cnv gex,meth,rppa meth,rppa,cnv meth,mutation,rppa gex,meth,mirna gex,cnv,mirna gex,mutation,mirna meth,mutation,mirna gex,mutation,cnv mutation,cnv,mirna meth,cnv,mirna gex,meth,cnv gex,meth,mutation meth,mutation,cnv gex,mutation,rppa,mirna meth,mutation,rppa,mirna gex,meth,rppa,mirna meth,rppa,cnv,mirna gex,meth,rppa,cnv gex,rppa,cnv,mirna mutation,rppa,cnv,mirna gex,mutation,rppa,cnv gex,meth,mutation,rppa meth,mutation,rppa,cnv gex,meth,cnv,mirna gex,meth,mutation,mirna gex,mutation,cnv,mirna gex,meth,mutation,cnv meth,mutation,cnv,mirna gex,meth,mutation,rppa,mirna meth,mutation,rppa,cnv,mirna gex,meth,mutation,rppa,cnv gex,mutation,rppa,cnv,mirna gex,meth,rppa,cnv,mirna gex,meth,mutation,cnv,mirna gex,meth,mutation,rppa,cnv,mirna

roc_auc_ovr 1.0 n_omics_layers = 1 n_omics_layers = 3

15

n_omics_layers = 5

0.9

0.8

0.7 integration late-concat late-mean early

0.6

pway_def

Figure 8: ROC AUC of one-vs-rest logistic regression for all combinations of omics layer in early, lateconcat and late-mean integration. The gray line denotes the highest median ROC AUC observed across all configurations tested.

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 9). 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

layers, in order of median impact, were RPPA, gene expression and µ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

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 10 and 26, 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 10). 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 µ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 26, 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

Discussion

In this work, we further evaluate our PAAE [da Costa Avelar et al., 2024] 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 Fortelny and Bock [2020] and Seninge et al. [2021], 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

16

1.0

n_omics_layers = 1

n_omics_layers = 3

n_omics_layers = 5

concordance_index concordance index

0.8

0.6 integration late-mean early late-concat

0.4

0.0

rppa meth mirna gex mutation cnv gex,rppa gex,mutation meth,rppa cnv,mirna rppa,mirna meth,cnv mutation,cnv rppa,cnv gex,mirna mutation,rppa meth,mutation gex,meth mutation,mirna meth,mirna gex,cnv gex,meth,cnv mutation,rppa,cnv meth,rppa,cnv gex,mutation,rppa gex,rppa,mirna meth,rppa,mirna meth,mutation,mirna gex,mutation,mirna mutation,rppa,mirna gex,meth,mirna meth,cnv,mirna mutation,cnv,mirna gex,meth,mutation rppa,cnv,mirna gex,cnv,mirna gex,mutation,cnv gex,rppa,cnv gex,meth,rppa meth,mutation,rppa meth,mutation,cnv meth,mutation,rppa,mirna gex,meth,rppa,cnv meth,mutation,rppa,cnv gex,mutation,rppa,mirna gex,meth,mutation,mirna gex,meth,mutation,rppa gex,meth,rppa,mirna gex,meth,mutation,cnv gex,mutation,rppa,cnv gex,mutation,cnv,mirna meth,rppa,cnv,mirna gex,rppa,cnv,mirna mutation,rppa,cnv,mirna meth,mutation,cnv,mirna gex,meth,cnv,mirna gex,meth,mutation,rppa,cnv gex,meth,mutation,cnv,mirna gex,meth,mutation,rppa,mirna gex,meth,rppa,cnv,mirna gex,mutation,rppa,cnv,mirna meth,mutation,rppa,cnv,mirna gex,meth,mutation,rppa,cnv,mirna

0.2

Figure 9: Concordance Index (C-index) of a Coxpway_def proportional hazards model for all omics layer combinations across for three tested integration methods, i.e. early, late-concat and late-mean. The gray line indicates the highest median C-index among all tested configurations (dotted line marks concordance of 0.5).

0.4 0.3 roc_auc_ovr_effect

roc_auc_ovr_effect

0.02 0.00 0.02 0.04 0.06 0.08

integration early late-concat late-mean gex

meth mutation rppa omics_layer

0.4

integration early late-concat late-mean

0.2 0.1 0.0 0.1

mirna

0.2 0.1 0.0 0.1

0.2 cnv

integration early late-concat late-mean

0.3 roc_auc_ovr_effect

0.04

0.2 gex

meth mutation rppa omics_layer

cnv

mirna

gex

meth mutation rppa omics_layer

cnv

mirna

Figure 10: Marginal contribution of each omics layer to the classification performance, measured by ROC AUC of a one-vs-rest logistic regression model, across the three tested integration methods: early, late-concat and late-mean. The left panel shows the median marginal contribution, the right panel displays a swarm plot of values across all combinations excluding that layer, and in the middle the same information is shown as a boxplot.

17

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 µ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 (i) a robust evaluation of the effect of dropout in a non-linear BINN; (ii) integrated interpretability, with pathway activity scores enabling downstream analyses such as feature importance, clustering, and survival prediction; (iii) a fully interpretable classification case study; and (iv) a thorough multi-omics extension, revealing the positive impact of multi-omics integration, the impact of individual omics layers, and establishing practical guidelines for integration within our framework.. 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 [Fortelny and Bock, 2020, Seninge et al., 2021, Selby et al., 2024], 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.[Yang et al., 2022, Chen et al., 2023, Cui et al., 2023, Liu et al., 2024]), 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 [Hollmann et al., 2023, 2025]) 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).

References Maryam Abbasi, Filipa G. Carvalho, Bernardete Ribeiro, and Joel P. Arrais. Predicting drug activity against cancer through genomic profiles and SMILES. Artificial Intelligence in Medicine, 150:102820, April 2024. ISSN 09333657. doi: 10.1016/j.artmed.2024.102820.

18

URL https://linkinghub.elsevier.com/ retrieve/pii/S0933365724000629.

2002. ISSN 1076-9757. doi: 10.1613/jair. 953. URL https://www.jair.org/index.php/ jair/article/view/10302.

Bulat Abdrakhimov, Emmanuel Kayewa, and Zhiwei Wang. Prediction of Acute Cardiac Rejec- Feng Chen, Michael C. Wendl, Matthew A. Wyczation Based on Gene Expression Profiles. Journal lkowski, Matthew H. Bailey, Yize Li, and Li Ding. of Personalized Medicine, 14(4):410, April 2024. Moving pan-cancer studies from basic research toISSN 2075-4426. doi: 10.3390/jpm14040410. URL ward the clinic. Nature Cancer, 2(9):879–890, https://www.mdpi.com/2075-4426/14/4/410. September 2021. ISSN 2662-1347. doi: 10.1038/ s43018-021-00250-4. URL https://www.nature. Daniel I. Benjamin, Alyssa Cozzo, Xiaodan Ji, Lindcom/articles/s43018-021-00250-4. say S. Roberts, Sharon M. Louie, Melinda M. Mulvihill, Kunxin Luo, and Daniel K. Nomura. Ether lipid generating enzyme AGPS alters the Jiawei Chen, Hao Xu, Wanyu Tao, Zhaoxiong Chen, Yuxuan Zhao, and Jing-Dong J. Han. Transbalance of structural and signaling lipids to fuel former for one stop interpretable cell type annocancer pathogenicity. Proceedings of the Natation. Nature Communications, 14(1):223, Jantional Academy of Sciences, 110(37):14912–14917, uary 2023. ISSN 2041-1723. doi: 10.1038/ September 2013. ISSN 0027-8424, 1091-6490. doi: s41467-023-35923-4. URL https://www.nature. 10.1073/pnas.1310894110. URL https://pnas. com/articles/s41467-023-35923-4. org/doi/full/10.1073/pnas.1310894110. Ann M. Bode and Zigang Dong. Precision Yuhao Chen and Xiaowei Wang. miRDB: an ononcology- the future of personalized cancer line database for prediction of functional mimedicine? npj Precision Oncology, 1(1):1–2, croRNA targets. Nucleic Acids Research, 48 March 2017. ISSN 2397-768X. doi: 10.1038/ (D1):D127–D131, January 2020. ISSN 0305s41698-017-0010-5. URL https://www.nature. 1048, 1362-4962. doi: 10.1093/nar/gkz757. com/articles/s41698-017-0010-5. Number: 1 URL https://academic.oup.com/nar/article/ Publisher: Nature Publishing Group. 48/D1/D127/5557729. Le Chang, Guangyan Zhou, Othman Soufan, and Travers Ching, Xun Zhu, and Lana X. Garmire. Jianguo Xia. miRNet 2.0: network-based visual anCox-nnet: An artificial neural network method alytics for miRNA functional analysis and systems for prognosis prediction of high-throughput omics biology. Nucleic Acids Research, 48(W1):W244– data. PLOS Computational Biology, 14(4): W251, July 2020. ISSN 0305-1048, 1362-4962. doi: e1006076, April 2018. ISSN 1553-7358. doi: 10.1093/nar/gkaa467. URL https://academic. 10.1371/journal.pcbi.1006076. URL https://dx. oup.com/nar/article/48/W1/W244/5850315. plos.org/10.1371/journal.pcbi.1006076. Kumardeep Chaudhary, Olivier B. Poirion, Liangqun Lu, and Lana X. Garmire. Deep Learning–Based Paola Antonia Corsetto, Stefania Zava, Angela Maria Rizzo, and Irma Colombo. The Critical Impact Multi-Omics Integration Robustly Predicts Surof Sphingolipid Metabolism in Breast Cancer Provival in Liver Cancer. Clinical Cancer Research, 24 gression and Drug Response. International Journal (6):1248–1259, March 2018. ISSN 1078-0432. doi: of Molecular Sciences, 24(3):2107, January 2023. 10.1158/1078-0432.CCR-17-0853. URL https:// ISSN 1422-0067. doi: 10.3390/ijms24032107. URL doi.org/10.1158/1078-0432.CCR-17-0853. https://www.mdpi.com/1422-0067/24/3/2107. N. V. Chawla, K. W. Bowyer, L. O. Hall, and W. P. Kegelmeyer. SMOTE: Synthetic Minor- D. R. Cox. Regression Models and Lifeity Over-sampling Technique. Journal of ArTables. Journal of the Royal Statistical Sotificial Intelligence Research, 16:321–357, June ciety Series B: Statistical Methodology, 34(2): 19

187–202, January 1972. ISSN 1369-7412, 14679868. doi: 10.1111/j.2517-6161.1972.tb00899. x. URL https://academic.oup.com/jrsssb/ article/34/2/187/7027194.

(BIBM), pages 864–871, Lisbon, Portugal, December 2024. IEEE. ISBN 9798350386226. doi: 10.1109/BIBM62325.2024.10822388. URL https: //ieeexplore.ieee.org/document/10822388/.

Nello Cristianini, John Shawe-Taylor, André Elis- Kountay Dwivedi, Ankit Rajpal, Sheetal Rajpal, seeff, and Jaz S. Kandola. On Kernel-Target Manoj Agarwal, Virendra Kumar, and Naveen KuAlignment. In Thomas G. Dietterich, Suzanna mar. An explainable AI-driven biomarker disBecker, and Zoubin Ghahramani, editors, Adcovery framework for Non-Small Cell Lung Canvances in Neural Information Processing Systems cer classification. Computers in Biology and 14 [Neural Information Processing Systems: NatMedicine, 153:106544, February 2023. ISSN ural and Synthetic, NIPS 2001, December 3-8, 00104825. doi: 10.1016/j.compbiomed.2023. 2001, Vancouver, British Columbia, Canada], 106544. URL https://linkinghub.elsevier. pages 367–373. MIT Press, 2001. URL https: com/retrieve/pii/S0010482523000094. //proceedings.neurips.cc/paper/2001/hash/ 1f71e393b3809197ed66df836fe833e5-Abstract. Yannan Fan and Jianguo Xia. miRNet—Functional Analysis and Visual Exploration of miRNA–Target html. Interactions in a Network Context. In Louise Haotian Cui, Chloe Wang, Hassaan Maan, Kuan Von Stechow and Alberto Santos Delgado, ediPang, Fengning Luo, and Bo Wang. scGPT: Totors, Computational Cell Biology, volume 1819, wards Building a Foundation Model for Single-Cell pages 215–233. Springer New York, New York, Multi-omics Using Generative AI. preprint, BioinNY, 2018. ISBN 978-1-4939-8617-0 978-1formatics, May 2023. URL http://biorxiv.org/ 4939-8618-7. doi: 10.1007/978-1-4939-8618-7 10. URL http://link.springer.com/10.1007/ lookup/doi/10.1101/2023.04.30.538439. 978-1-4939-8618-7_10. Series Title: Methods in Shidong Cui, Sicong Yu, Hsi-Yuan Huang, Yang-ChiMolecular Biology. Dung Lin, Yixian Huang, Bojian Zhang, Jihan Xiao, Huali Zuo, Jiayi Wang, Zhuoran Li, Guang- Yannan Fan, Keith Siklenka, Simran K. Arora, Paula Ribeiro, Sarah Kimmins, and Jianguo hao Li, Jiajun Ma, Baiming Chen, Haoxuan Zhang, Xia. miRNet - dissecting miRNA-target inJiehui Fu, Liang Wang, and Hsien-Da Huang. miRteractions and functional associations through TarBase 2025: updates to the collection of exnetwork-based visual analysis. Nucleic Acids Reperimentally validated microRNA–target interacsearch, 44(W1):W135–W141, July 2016. ISSN tions. Nucleic Acids Research, 53(D1):D147–D156, 0305-1048, 1362-4962. doi: 10.1093/nar/ January 2025. ISSN 0305-1048, 1362-4962. doi: gkw288. URL https://academic.oup.com/nar/ 10.1093/nar/gkae1072. URL https://academic. article-lookup/doi/10.1093/nar/gkw288. oup.com/nar/article/53/D1/D147/7907368. Pedro Henrique da Costa Avelar, Min Wu, and Nikolaus Fortelny and Christoph Bock. Knowledgeprimed neural networks enable biologically Sophia Tsoka. Incorporating Prior Knowledge in interpretable deep learning on single-cell seDeep Learning Models via Pathway Activity Auquencing data. Genome Biology, 21(1): toencoders, June 2023. URL http://arxiv.org/ 190, December 2020. ISSN 1474-760X. abs/2306.05813. arXiv:2306.05813 [cs]. doi: 10.1186/s13059-020-02100-5. URL Pedro Henrique da Costa Avelar, Le Ou-Yang, Min https://genomebiology.biomedcentral.com/ Wu, and Sophia Tsoka. Pathway Activity Autoenarticles/10.1186/s13059-020-02100-5. coders for Enhanced Omics Analysis and Clinical Interpretability. In 2024 IEEE International Hemaasri-Neya Girithar, Ananda Staats Pires, Seong Beom Ahn, Gilles J. Guillemin, Laurence Conference on Bioinformatics and Biomedicine 20

Gluch, and Benjamin Heng. Involvement of https://linkinghub.elsevier.com/retrieve/ the kynurenine pathway in breast cancer: uppii/S1532046422000429. dates on clinical research and trials. British Journal of Cancer, 129(2):185–203, August 2023. M. G. Kendall. THE TREATMENT OF TIES IN RANKING PROBLEMS. Biometrika, 33 ISSN 0007-0920, 1532-1827. doi: 10.1038/ (3):239–251, 1945. ISSN 0006-3444, 1464-3510. s41416-023-02245-7. URL https://www.nature. doi: 10.1093/biomet/33.3.239. URL https: com/articles/s41416-023-02245-7. //academic.oup.com/biomet/article-lookup/ G. E. Hinton and R. R. Salakhutdinov. Reducdoi/10.1093/biomet/33.3.239. ing the Dimensionality of Data with Neural NetAutoworks. Science, 313(5786):504–507, July 2006. Diederik P. Kingma and Max Welling. Encoding Variational Bayes. In Yoshua Bengio ISSN 0036-8075, 1095-9203. doi: 10.1126/science. and Yann LeCun, editors, 2nd International Con1127647. URL https://www.science.org/doi/ ference on Learning Representations, ICLR 2014, 10.1126/science.1127647. Banff, AB, Canada, April 14-16, 2014, Conference Noah Hollmann, Samuel Müller, Katharina Track Proceedings, 2014. URL http://arxiv. Eggensperger, and Frank Hutter. TabPFN: org/abs/1312.6114. A Transformer That Solves Small Tabular Classification Problems in a Second. In The Eleventh Simon Kornblith, Mohammad Norouzi, Honglak Lee, and Geoffrey E. Hinton. Similarity of International Conference on Learning RepreNeural Network Representations Revisited. In sentations, ICLR 2023, Kigali, Rwanda, May Kamalika Chaudhuri and Ruslan Salakhutdi1-5, 2023. OpenReview.net, 2023. URL https: nov, editors, Proceedings of the 36th Interna//openreview.net/forum?id=cp5PvcI6w8_. tional Conference on Machine Learning, ICML Noah Hollmann, Samuel Müller, Lennart Purucker, 2019, 9-15 June 2019, Long Beach, CaliforArjun Krishnakumar, Max Körfer, Shi Bin Hoo, nia, USA, volume 97 of Proceedings of MaRobin Tibor Schirrmeister, and Frank Hutter. Acchine Learning Research, pages 3519–3529. PMLR, curate predictions on small data with a tabu2019. URL http://proceedings.mlr.press/ lar foundation model. Nature, 637(8045):319–326, v97/kornblith19a.html. January 2025. ISSN 0028-0836, 1476-4687. doi: 10.1038/s41586-024-08328-6. URL https://www. Yutong Li, Jonathan Cardoso-Silva, Lazaros G. Papageorgiou, and Sophia Tsoka. Optimisationnature.com/articles/s41586-024-08328-6. based modelling for drug discovery in malaria. Yi-Fan Huang, Shuji Mizumoto, and Morihisa preprint, Bioinformatics, February 2022. URL Fujita. Novel Insight Into Glycosaminoglyhttp://biorxiv.org/lookup/doi/10.1101/ can Biosynthesis Based on Gene Expression 2022.02.12.479469. Profiles. Frontiers in Cell and Developmental Biology, 9:709018, September 2021. ISSN Stan Lipovetsky and Michael Conklin. Analysis of regression in game theory approach. Ap2296-634X. doi: 10.3389/fcell.2021.709018. plied Stochastic Models in Business and InURL https://www.frontiersin.org/articles/ dustry, 17(4):319–330, October 2001. ISSN 10.3389/fcell.2021.709018/full. 1524-1904, 1526-4025. doi: 10.1002/asmb. Ishleen Kaur, M.N. Doja, and Tanvir Ahmad. 446. URL https://onlinelibrary.wiley.com/ Data mining and machine learning in cancer doi/10.1002/asmb.446. survival research: An overview and future recommendations. Journal of Biomedical Weijun Liu and Xiaowei Wang. Prediction of funcInformatics, 128:104026, April 2022. ISSN tional microRNA targets by integrative model15320464. doi: 10.1016/j.jbi.2022.104026. URL ing of microRNA binding and target expression 21

data. Genome Biology, 20(1):18, December 2019. ISSN 1474-760X. doi: 10.1186/s13059-019-1629-z. URL https://genomebiology.biomedcentral. com/articles/10.1186/s13059-019-1629-z.

Kelley, and David P. Bartel. The biochemical basis of microRNA targeting efficacy. Science, 366(6472):eaav1741, December 2019. ISSN 0036-8075, 1095-9203. doi: 10.1126/science. aav1741. URL https://www.science.org/doi/ 10.1126/science.aav1741.

Yanguang Liu, Hailong Yu, Xinya Duan, Xiaomin Zhang, Ting Cheng, Feng Jiang, Hao Tang, Yao Ruan, Miao Zhang, Hongyu Zhang, Wanjiru Muriithi, Lucy Wanjiku Macharia, Carlos and Qingye Zhang. TransGEM: a molecule Pilotto Heming, Juliana Lima Echevarria, Atunga generation model based on Transformer with Nyachieo, Paulo Niemeyer Filho, and Vivaldo gene expression data. Bioinformatics, 40 Moura Neto. ABC transporters and the hall(5):btae189, May 2024. ISSN 1367-4811. marks of cancer: roles in cancer aggressiveness doi: 10.1093/bioinformatics/btae189. URL beyond multidrug resistance. Cancer Biology and https://academic.oup.com/bioinformatics/ Medicine, 17(2):253–269, 2020. ISSN 2095-3941. article/doi/10.1093/bioinformatics/ doi: 10.20892/j.issn.2095-3941.2019.0284. URL btae189/7649318. http://www.cancerbiomed.org/lookup/doi/ 10.20892/j.issn.2095-3941.2019.0284. Scott M. Lundberg and Su-In Lee. A Unified Approach to Interpreting Model Predictions. In Marco Tulio Ribeiro, Sameer Singh, and Carlos Isabelle Guyon, Ulrike von Luxburg, Samy Bengio, Guestrin. ”Why Should I Trust You?”: Explaining Hanna M. Wallach, Rob Fergus, S. V. N. Vishthe Predictions of Any Classifier. In Proceedings wanathan, and Roman Garnett, editors, Advances of the 22nd ACM SIGKDD International Conferin Neural Information Processing Systems 30: Anence on Knowledge Discovery and Data Mining, nual Conference on Neural Information Processing pages 1135–1144, San Francisco California USA, Systems 2017, December 4-9, 2017, Long Beach, August 2016. ACM. ISBN 978-1-4503-4232-2. doi: CA, USA, pages 4765–4774, 2017. URL https: 10.1145/2939672.2939778. URL https://dl.acm. //proceedings.neurips.cc/paper/2017/hash/ org/doi/10.1145/2939672.2939778. 8a20a8621978632d76c43dfd28b67767-Abstract. Lindsay K. Ryland, Todd E. Fox, Xin Liu, Thomas P. html. Loughran, and Mark Kester. Dysregulation of Huanhuan Ma, Ran Ding, Junwen Wang, Guangying sphingolipid metabolism in cancer. Cancer BiDu, Yun Zhang, Qiuchen Lu, Yingyue Hou, ology & Therapy, 11(2):138–149, January 2011. Haosong Chen, and Hongguan Jiao. Global ISSN 1538-4047, 1555-8576. doi: 10.4161/cbt.11.2. research trends in tryptophan metabolism 14624. URL http://www.tandfonline.com/doi/ and cancer: a bibliometric and visualizaabs/10.4161/cbt.11.2.14624. tion analysis (2005–2024). Frontiers in Oncology, 15:1621666, July 2025. ISSN 2234- David A. Selby, Maximilian Sprang, Jan Ewald, and Sebastian J. Vollmer. Beyond the black 943X. doi: 10.3389/fonc.2025.1621666. URL box with biologically informed neural networks. https://www.frontiersin.org/articles/10. Nature Reviews Genetics, 26(6):371–372, June 3389/fonc.2025.1621666/full. 2025. ISSN 1471-0056, 1471-0064. doi: 10.1038/ Christopher D. Manning, Prabhakar Raghavan, and s41576-025-00826-1. URL https://www.nature. Hinrich Schütze. Introduction to information recom/articles/s41576-025-00826-1. trieval. Cambridge university press, Cambridge, David Antony Selby, Rashika Jakhmola, Maximil2008. ISBN 978-0-521-86571-5. ian Sprang, Gerrit Grossmann, Hind Raki, NilooSean E. McGeary, Kathy S. Lin, Charlie Y. far Maani, Daria Pavliuk, Jan Ewald, and SeShi, Thy M. Pham, Namita Bisaria, Gina M. bastian J Vollmer. Visible neural networks for 22

multi-omics integration: a critical review, December 2024. URL http://biorxiv.org/lookup/ doi/10.1101/2024.12.09.627465.

Lei Xiang, Peng Su, Shujun Xia, Zhiyan Liu, Yan Wang, Peng Gao, and Genyin Zhou. ABCG2 is associated with HER-2 Expression, lymph node metastasis and clinical stage in breast Lucas Seninge, Ioannis Anastopoulos, Hongxu invasive ductal carcinoma. Diagnostic PatholDing, and Joshua Stuart. VEGA is an inogy, 6(1):90, December 2011. ISSN 1746-1596. terpretable generative model for inferring biodoi: 10.1186/1746-1596-6-90. URL https: logical network activity in single-cell transcrip//diagnosticpathology.biomedcentral.com/ tomics. Nature Communications, 12(1):5684, articles/10.1186/1746-1596-6-90. September 2021. ISSN 2041-1723. doi: 10.1038/ s41467-021-26017-0. URL https://www.nature. Xinping Xie, Fengting Wang, Guanfu Wang, Weicom/articles/s41467-021-26017-0. wei Zhu, Xiaodong Du, and Hongqiang Wang. Learning the cellular activity representation based L. S. Shapley. A Value for n-Person Games. In on gene regulatory networks for prediction of Harold William Kuhn and Albert William Tucker, tumor response to drugs. Artificial Inteleditors, Contributions to the Theory of Games ligence in Medicine, 152:102864, June 2024. (AM-28), Volume II, pages 307–318. Princeton ISSN 09333657. doi: 10.1016/j.artmed.2024. University Press, December 1953. ISBN 978-1102864. URL https://linkinghub.elsevier. 4008-8197-0. doi: 10.1515/9781400881970-018. com/retrieve/pii/S0933365724001064. URL https://www.degruyter.com/document/ doi/10.1515/9781400881970-018/html. Linxuan Xue, Chao Wang, Yulu Qian, WenBora Uyar, Jonathan Ronen, Vedran Franke, Gaetano Gargiulo, and Altuna Akalin. Multi-omics and deep learning provide a multifaceted view of cancer. Technical report, bioRxiv, October 2021. URL https://www.biorxiv.org/content/ 10.1101/2021.09.29.462364v1. Section: New Results Type: article.

qiang Zhu, Lina Liu, Xiaohong Yang, Shuhua Zhang, and Daya Luo. Tryptophan metabolism regulates inflammatory macrophage polarization as a predictive factor for breast cancer immunotherapy. International Immunopharmacology, 125:111196, December 2023. ISSN 15675769. doi: 10.1016/j.intimp.2023.111196. URL https://linkinghub.elsevier.com/ retrieve/pii/S1567576923015230.

Arno Van Hilten, Steven A. Kushner, Manfred Kayser, M. Arfan Ikram, Hieab H. H. Adams, Caroline C. W. Klaver, Wiro J. Niessen, and Gen- Fan Yang, Wenchuan Wang, Fang Wang, Yuan Fang, Duyu Tang, Junzhou Huang, Hui Lu, and Jianhua nady V. Roshchupkin. GenNet framework: interYao. scBERT as a large-scale pretrained deep lanpretable deep learning for predicting phenotypes guage model for cell type annotation of single-cell from genetic data. Communications Biology, 4 RNA-seq data. Nature Machine Intelligence, 4(10): (1):1094, September 2021. ISSN 2399-3642. doi: 852–866, September 2022. ISSN 2522-5839. doi: 10.1038/s42003-021-02622-z. URL https://www. 10.1038/s42256-022-00534-z. URL https://www. nature.com/articles/s42003-021-02622-z. nature.com/articles/s42256-022-00534-z. David Wissel, Daniel Rowson, and Valentina Boeva. Hierarchical autoencoder-based integra- Melvyn Yap, Rebecca L. Johnston, Helena Foley, tion improves performance in multi-omics canSamual MacDonald, Olga Kondrashova, Khoa A. cer survival models through soft modality seTran, Katia Nones, Lambros T. Koufariotis, lection. Technical report, bioRxiv, January Cameron Bean, John V. Pearson, Maciej Trza2022. URL https://www.biorxiv.org/content/ skowski, and Nicola Waddell. Verifying explain10.1101/2021.09.16.460589v3. Section: New ability of a deep learning tissue classifier trained Results Type: article. on RNA-seq data. Scientific Reports, 11(1):2641, 23

January 2021. ISSN 2045-2322. doi: 10.1038/ s41598-021-81773-9. URL https://www.nature. com/articles/s41598-021-81773-9. Hung-Rong Yen, Wen-Chieh Liao, Chia-Hua Chen, Ying-Ai Su, Ying-Wei Huang, Chi Hsiao, Yu-Lun Chou, Yin-Hung Chu, Pin-Keng Shih, and ChiungHui Liu. Targeting chondroitin sulfate suppresses macropinocytosis of breast cancer cells by modulating syndecan-1 expression. Molecular Oncology, pages 1878–0261.13667, May 2024. ISSN 15747891, 1878-0261. doi: 10.1002/1878-0261.13667. URL https://febs.onlinelibrary.wiley.com/ doi/10.1002/1878-0261.13667. Shibo Yu, Xiaowen Wang, Lizhe Zhu, Peiling Xie, Yudong Zhou, Siyuan Jiang, Heyan Chen, Xiaoqin Liao, Shengyu Pu, Zhenzhen Lei, Bin Wang, and Yu Ren. A systematic analysis of a potential metabolism-related prognostic signature for breast cancer patients. Annals of Translational Medicine, 9(4):330–330, February 2021. ISSN 23055839, 23055847. doi: 10.21037/ atm-20-7600. URL https://atm.amegroups. com/article/view/62307/html. Oana A Zeleznik, Raji Balasubramanian, Yumeng Ren, Deirdre K Tobias, Bernard A Rosner, Cheng Peng, Alaina M Bever, Lisa Frueh, Sarah Jeanfavre, Julian Avila-Pacheco, Clary B Clish, Samia Mora, Frank B Hu, and A Heather Eliassen. Branched-Chain Amino Acids and Risk of Breast Cancer. JNCI Cancer Spectrum, 5(5):pkab059, September 2021. ISSN 2515-5091. doi: 10.1093/jncics/pkab059. URL https://academic.oup.com/jncics/article/ doi/10.1093/jncics/pkab059/6321892.

24

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

Alternatives for Omics-Pathway Mappings

We briefly considered other alternative databases for non-gene-based omics layers. While the dataset provided by Wissel et al. [2022] 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 µRNA, we briefly considered using predicted target genes instead of the experimentally validated genes of miRTarBase2025 miRTarBase SE WR [Cui et al., 2025]. That is, we considered using information from either miRDB [Liu and Wang, 2019, Chen and Wang, 2020], TargetScanHuman [McGeary et al., 2019], or miRNet [Fan et al., 2016, Fan and Xia, 2018, Chang et al., 2020]. However, we noticed that when using predicted targets the number of target genes for a significant portion of µRNA measurements in our dataset grew significantly, and compounded over each gene possibly being into more than one pathway. This resulted in every µRNA being associated, on average, with 2.3 pathways when following Chen and Wang [2020]’s recommendations. This translated in the memory requirement for training an early integration model almost doubling when adding the µRNA to a model with only Gene Expression, CNV, Mutation, and RPPA, despite the number of unique features in the µRNA layer being less than 5% of the total unique features.

B

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 Ts = 32 and with annealing duration 128 for Te = Ts + 128 in our annealed scheduling function STs ,Te (t). We use a gaussian distribution defined by µ and σ, from which the values of z are sampled using the reparametrization trick to maintain differentiability [Kingma and Welling, 2014]. 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) 25

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 µ), 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 × 64 units (64 for µ and 64 for σ); 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 × 64 units (64 for µ and 64 for σ); 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, α = 1, λ = 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 p kernel. For the Random Forest Classifier we build 100 Trees using Gini impurity, considering dinput 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 µ 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 26

27

PASL

PathME

PACL

t

m 2000

1 3

10−5 , 10−5 , 10−1

10−9 , 10−5 , 10−1 λ

Learning Rate

N/A N/A [],[32],[32,16] [],[32],[32,16]

Pathway Module Hidden Layer Sizes

λ

16,4,1

RBM k

N/A χt>32 (t), S32,160 (t) N/A χt>32 (t), S32,160 (t)

β Schedule

Table 1: Hyperparameters used during the internal validation grid search.

0.9

a1+a2 350+150

0,1

α

Layer Sizes [16],[32],[64]

16,64,1024,256

[64], [128,64], [256,128,64], [500,200,64]

Epochs

Layer Sizes

AE VAE PAAE PAVAE

β N/A 1,5,10,50,100 N/A 1,5,10,50,100

Encoder Layer Sizes

[128,64], [256,128,64], [512,256,128,64] [128,64], [256,128,64], [512,256,128,64] [64],[128,64] [128,64], [256,128,64]

Model

Classifier

LR,SVM,RF

LR,SVM,RF

LR,SVM,RF

LR,SVM,RF LR,SVM,RF LR,SVM,RF LR,SVM,RF

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.

B.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 (a) that our PAAE model and PathME have to reconstruct the full input; then (b) when we consider that our model only has to reconstruct the pathway set’s member genes; and, finally, (c) if we are only interested in dimensionality reduction, what would be our model size if we only kept the Encoding part and did not keep the reconstruction part. 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 (a), 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 (b), 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 (a) as we’d only need to prune the decoder, only needing these extra parameters for training. If we only need the encoder (c), 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 (a) since we’d only need to prune our model and discard these parameters after training.

B.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 28

Model

Configuration

Parameters

AE

[|G| , 128, 64, 128, |G|] [|G| , 256, 128, 64, 128, 256, |G|] [|G| , 512, 256, 128, 64, 128, 256, 512, |G|]

Space

Dimensionality

No Prior Knowledge 4,795,105 9,620,961 19,403,745

z z z

64 64 64

Table 2: Model parameter sizes for every pararameter configuration tested in our hyperparameter search for the Naive (no prior knowledge) models. Please see Tables 3 and 4 for other models and Section B.1 for discussion on these results.

a comparison of improved versions of these methods’ pipelines to future work due to computational budget constraints.

C

Interpreting the Pathway Activity Scores

As we can see in the original PAAE description [da Costa Avelar et al., 2024], 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 Epj , if it has k layers, can be described as in Eq. 11: aj = Epj (x) = hj,k hj,k = Wj,k × hj,k−1 + bj,k hj,i = fj,k (Wj,i × hj,i−1 + bj,i )

(11)

hj,1 = fj,k (Wj,1 × x:,pj + bj,1 ) Where Wj,i is a ℜdj,i−1 ×dj,i matrix, and fj,k is a nonlinear function, not present in the last layer so that our output is roughly normal-shaped. Since the weights Wj,i can be both positive and negative, one should not expect the final result aj to have a positive/negative directionality w.r.t. the pathway’s input x:,pj , 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 aj is the value that we use in both our clustermap (Figs. 12, 11) and featuremap (Fig. 3) 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 aj , one should not expect the values to behave linearly.

D

Clustermaps

For clustermaps, in Figs. 12 and 11, 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. 29

Model

Configuration

Parameters

Space

Dimensionality

Reconstruct Full Gene Matrix
|P | × [dp , 1] → [64, |G|]

1,219,041

|P | × [dp , 32, 1] → [64, |G|]

1,444,883

|P | × [dp , 32, 16, 1] → [64, |G|]

1,470,483

|P | × [dp , 1] → [128, 64, 128, |G|]

2,428,833

|P | × [dp , 32, 1] → [128, 64, 128, |G|]

2,654,675

PAAE

PathME

|P | × [dp , 32, 16, 1] → [128, 64, 128, |G|]

2,680,275

|P | × [dp , 16, 1, 16, G] |P | × [dp , 32, 1, 32, G] |P | × [dp , 64, 1, 64, G] |P | × [dp , 128, 1, 128, G] |P | × [dp , 256, 1, 256, G]

43,514,144 87,028,288 174,056,576 348,113,152 696,226,304

z a z a z a z a z a z a a a a a a

64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 |P | = 50

Reconstructing Only Pathway’s Input Genes
|P | × [dp , 1] → [64, |G|]

289,411

|P | × [dp , 32, 1] → [64, |G|]

515,253

|P | × [dp , 32, 16, 1] → [64, |G|]

540,853

|P | × [dp , 1] → [128, 64, 128, |G|]

583,875

|P | × [dp , 32, 1] → [128, 64, 128, |G|]

809,717

|P | × [dp , 32, 16, 1] → [128, 64, 128, |G|]

835,317

|P | × [dp , 16, 1, 16, dp ] |P | × [dp , 32, 1, 32, dp ] |P | × [dp , 64, 1, 64, dp ] |P | × [dp , 128, 1, 128, dp ] |P | × [dp , 256, 1, 256, dp ]

233,088 466,176 932,352 1,864,704 3,729,408

PAAE

PathME

z a z a z a z a z a z a a a a a a

64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 |P | = 50

Encoder Only

|P | × [dp , 1] → [500, 200, 64]

10,496 7,232 236,338 233,074 261,938 258,674 22,016 7,232 247,858 233,074 273,458 258,674 14,837 11,523 26,485 11,523 66,165 11,523 150,837 11,523

[D] = [D1 + D2 ] = [500] = [350 + 150]

354,197

|P | × [dp , 1] → [64, |G|]
|P | × [dp , 32, 1] → [64, |G|]
|P | × [dp , 32, 16, 1] → [64, |G|] PAAE
|P | × [dp , 1] → [128, 64, 128, |G|]
|P | × [dp , 32, 1] → [128, 64, 128, |G|]
|P | × [dp , 32, 16, 1] → [128, 64, 128, |G|]
|P | × [dp , 1] → [64]
|P | × [dp , 1] → [128, 64] PACL
|P | × [dp , 1] → [256, 128, 64]

PASL

z a z a z a z a z a z a z a z a z a z a a′ z′ L

64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 64 |P | = 50 D1 = 350 D2 = 150 D = 500

Table 3: Model parameter sizes for every pararameter configuration tested in our hyperparameter search for the models using the Hallmark Genes pathway set. Please see Tables 2 and 4 for other models and Section B.1 for discussion on these results.

30

Model

Configuration

Parameters

Space

Dimensionality

Reconstruct Full Gene Matrix
|P | × [dp , 1] → [64, |G|]

1,233,085

|P | × [dp , 32, 1] → [64, |G|]

1,628,955

|P | × [dp , 32, 16, 1] → [64, |G|]

1,724,187

|P | × [dp , 1] → [128, 64, 128, |G|]

2,451,581

|P | × [dp , 32, 1] → [128, 64, 128, |G|]

2,847,451

PAAE

|P | × [dp , 32, 16, 1] → [128, 64, 128, |G|]

PathME

|P | × [dp , 16, 1, 16, G] |P | × [dp , 32, 1, 32, G] |P | × [dp , 64, 1, 64, G] |P | × [dp , 128, 1, 128, G]

2,942,683 161,643,344 323,286,688 646,573,376 1,293,146,752

z a z a z a z a z a z a a a a a

64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 |P | = 186

Reconstruct Only Pathway’s Input Genes
|P | × [dp , 1] → [64, |G|]

350,385

|P | × [dp , 32, 1] → [64, |G|]

746,255

|P | × [dp , 32, 16, 1] → [64, |G|]

841,487

PAAE

PathME

|P | × [dp , 1] → [128, 64, 128, |G|]

699,761

|P | × [dp , 32, 1] → [128, 64, 128, |G|]

1,095,631

|P | × [dp , 32, 16, 1] → [128, 64, 128, |G|]

1,190,863

|P | × [dp , 16, 1, 16, dp ] |P | × [dp , 32, 1, 32, dp ] |P | × [dp , 64, 1, 64, dp ] |P | × [dp , 128, 1, 128, dp ]

408,544 817,088 1,634,176 3,268,352

z a z a z a z a z a z a a a a a

64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 |P | = 186

Encoder Only

|P | × [dp , 1] → [500, 200, 64]

24,540 12,572 420,410 408,442 515,642 503,674 44,764 12,572 440,634 408,442 535,866 503,674 29,739 17,585 50,091 17,585 107,179 17,585 225,035 17,585

[D] = [D1 + D2 ] = [500] = [350 + 150]

330,038

|P | × [dp , 1] → [64, |G|]
|P | × [dp , 32, 1] → [64, |G|]
|P | × [dp , 32, 16, 1] → [64, |G|] PAAE
|P | × [dp , 1] → [128, 64, 128, |G|]
|P | × [dp , 32, 1] → [128, 64, 128, |G|]
|P | × [dp , 32, 16, 1] → [128, 64, 128, |G|]
|P | × [dp , 1] → [64]
|P | × [dp , 1] → [128, 64] PACL
|P | × [dp , 1] → [256, 128, 64]

PASL

z a z a z a z a z a z a z a z a z a z a a′ z′ L

64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 64 |P | = 186 D1 = 350 D2 = 150 D = 500

Table 4: Model parameter sizes for every pararameter configuration tested in our hyperparameter search for the models using the KEGG pathway set. Please see Tables 2 and 3 for other models and Section B.1 for discussion on these results.

31

E

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.

F

Survival Analyses

For survival plots and statistical tests (Fig. 4) a time window limit of 5 years (365 × 5 = 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 F.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 F.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.

F.1

Preprocessing for Gene Expression Comparisons

For the analyses involving raw input features, such as those in Fig. 4, we re-normalise the genes from log(F P KM + 1) to log(T P M + 1) in the TCGA dataset, and from log(Intensity + 1) to log(IP M + 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 Ig ×106 intensity is normalised per sample as if it was a raw count value with IP Mg = P Ig . g

F.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) [Uyar et al., 2021]. NPWj =

k Y

Wj,i , NPWj ∈ ℜ|pj |

(12)

i=1

ANPWj [g] = | NPWj [g]| 32

(13)

PAM50 PEROXISOME IL2_STAT5_SIGNAL PI3K_AKT_MTOR_SI DNA_REPAIR PANCREAS_BETA_CE ESTROGEN_RESPONS KRAS_SIGNALING_D UNFOLDED_PROTEIN G2M_CHECKPOINT WNT_BETA_CATENIN NOTCH_SIGNALING HYPOXIA INTERFERON_GAMMA REACTIVE_OXYGEN_ BILE_ACID_METABO MITOTIC_SPINDLE COMPLEMENT EPITHELIAL_MESEN HEDGEHOG_SIGNALI INFLAMMATORY_RES APICAL_JUNCTION FATTY_ACID_METAB ESTROGEN_RESPONS XENOBIOTIC_METAB MYC_TARGETS_V1 MYOGENESIS ANDROGEN_RESPONS PROTEIN_SECRETIO CHOLESTEROL_HOME MTORC1_SIGNALING P53_PATHWAY UV_RESPONSE_DN UV_RESPONSE_UP ADIPOGENESIS E2F_TARGETS ALLOGRAFT_REJECT OXIDATIVE_PHOSPH MYC_TARGETS_V2 APOPTOSIS TNFA_SIGNALING_V APICAL_SURFACE IL6_JAK_STAT3_SI HEME_METABOLISM ANGIOGENESIS SPERMATOGENESIS INTERFERON_ALPHA TGF_BETA_SIGNALI KRAS_SIGNALING_U GLYCOLYSIS COAGULATION

PAM50Call_RNAseq TNFA_SIGNALING_V KRAS_SIGNALING_U MITOTIC_SPINDLE KRAS_SIGNALING_D ESTROGEN_RESPONS IL2_STAT5_SIGNAL EPITHELIAL_MESEN DNA_REPAIR BILE_ACID_METABO UNFOLDED_PROTEIN INTERFERON_GAMMA REACTIVE_OXYGEN_ PEROXISOME COMPLEMENT HEDGEHOG_SIGNALI HYPOXIA ESTROGEN_RESPONS WNT_BETA_CATENIN G2M_CHECKPOINT PI3K_AKT_MTOR_SI NOTCH_SIGNALING PANCREAS_BETA_CE GLYCOLYSIS XENOBIOTIC_METAB ADIPOGENESIS APOPTOSIS MYC_TARGETS_V1 MTORC1_SIGNALING CHOLESTEROL_HOME PROTEIN_SECRETIO ANDROGEN_RESPONS P53_PATHWAY ALLOGRAFT_REJECT HEME_METABOLISM MYOGENESIS E2F_TARGETS IL6_JAK_STAT3_SI COAGULATION UV_RESPONSE_UP APICAL_SURFACE TGF_BETA_SIGNALI INTERFERON_ALPHA SPERMATOGENESIS ANGIOGENESIS MYC_TARGETS_V2 OXIDATIVE_PHOSPH APICAL_JUNCTION FATTY_ACID_METAB UV_RESPONSE_DN INFLAMMATORY_RES SampleID

(a) TCGA (train)

(b) Metabric (test)

Figure 11: The clustermap using the cosine distance between samples’ inferred pathway activity vectors for Hallmark Genes PAAE’s pathway activity space with the colours marking BRCA’s 5-classes: Normal in blue, Luminal A in orange, Luminal B in green, Basal in red, Her2 in purple. Best seen in color. 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.

G

Supporting Figures

33

PAM50 DORSO_VENTRAL_AX MISMATCH_REPAIR PYRIMIDINE_METAB STEROID_HORMONE_ VASOPRESSIN_REGU GLYCOSPHINGOLIPI VALINE_LEUCINE_A GLIOMA ARRHYTHMOGENIC_R NUCLEOTIDE_EXCIS ETHER_LIPID_META PANCREATIC_CANCE ENDOMETRIAL_CANC GLYCOSAMINOGLYCA FATTY_ACID_METAB HYPERTROPHIC_CAR OOCYTE_MEIOSIS RENIN_ANGIOTENSI ECM_RECEPTOR_INT PRIMARY_IMMUNODE P53_SIGNALING_PA ONE_CARBON_POOL_ PENTOSE_PHOSPHAT TRYPTOPHAN_METAB N_GLYCAN_BIOSYNT NICOTINATE_AND_N DNA_REPLICATION PROSTATE_CANCER RNA_DEGRADATION SPHINGOLIPID_MET CHRONIC_MYELOID_ PROXIMAL_TUBULE_

PAM50Call_RNAseq HISTIDINE_METABO FRUCTOSE_AND_MAN GLYCOSPHINGOLIPI VALINE_LEUCINE_A GLIOMA LONG_TERM_DEPRES GLYCOSAMINOGLYCA AMINOACYL_TRNA_B PROXIMAL_TUBULE_ SPHINGOLIPID_MET CHRONIC_MYELOID_ GLYCEROPHOSPHOLI LEUKOCYTE_TRANSE MISMATCH_REPAIR FATTY_ACID_METAB LIMONENE_AND_PIN PANCREATIC_CANCE ETHER_LIPID_META GLYCOSPHINGOLIPI GLYCOLYSIS_GLUCO P53_SIGNALING_PA ONE_CARBON_POOL_ TRYPTOPHAN_METAB RNA_DEGRADATION EPITHELIAL_CELL_ HYPERTROPHIC_CAR PEROXISOME RENIN_ANGIOTENSI DORSO_VENTRAL_AX GLYCOSAMINOGLYCA PHENYLALANINE_ME ABC_TRANSPORTERS SampleID

(a) TCGA (train)

(b) Metabric (test)

Figure 12: The clustermap using the cosine distance between samples’ inferred pathway activity vectors for KEGG PAAE’s pathway activity space with the colours marking BRCA’s 5-classes: Normal in blue, Luminal A in orange, Luminal B in green, Basal in red, Her2 in purple. Best seen in color.

Figure 13: Average ROC AUC (left) from a one-vs-rest Logistic Regression and Concordance Index (right) of a Cox-PH model trained on the pathway activity vector of a PAAE model using the Hallmark Genes pathway set, trained on the MO-TCGA-BRCA dataset. Performance in shown across varying dropout rates applied on the pathway-activity layer. Both metrics show a gentle downward slope, with a sharper decrease between 70 and 80% dropout. Shaded regions indicate model variance.

34

Figure 14: The learned pathway space when applying, from left to right, 0%, 10%, 50%, 80%, and 90% dropout on the Pathway Activity Layer, clustered according to the result from 0% (top 5) and 90% dropout (bottom 5). One can see that, although dropout on the latent space has the ability to improve reproducibility, this comes at the cost of “bleaching” the learned representation.

Figure 15: The average Kendall τ -b coefficient between pathway rankings obtained from models trained with different random seeds on the MO-TCGA-BRCA dataset, based on mutual information. We compare nonconsensus (reg) models between themselves (top 3), non-consensus against aligned consensus models (cons) and ACMs among themselves. The correlation is taken between models with p dropout (x-axis) and models with p = 0 (left), p = 0.5 (middle) and p = 0.9 (right) dropout. We can see that representational similarity seems to be consistent until a phase transition happens between 70% and 80% dropout.

35

Figure 16: The average Kendall τ -b coefficient between pathway rankings obtained from models trained with different random seeds on the MO-TCGA-BRCA dataset, based on concordance index. We compare nonconsensus (reg) models between themselves (top 3), non-consensus against aligned consensus models (cons) and ACMs among themselves. The correlation is taken between models with p dropout (x-axis) and models with p = 0 (left), p = 0.5 (middle) and p = 0.9 (right) dropout. We can see that representational similarity seems to be consistent until a phase transition happens between 70% and 80% dropout.

36

Figure 17: Boxplots of the concordance index of a pathway obtained from models trained with different random seeds on the MO-TCGA-BRCA dataset. Models are labelled as either Aligned Consensus Model (bottom) or not (top). From left to right, we show the results for dropout p being 0.0, 0.1, 0.5, 0.8, and 0.9. Note that model order changes drastically at the last two dropout values.

37

Figure 18: Boxplots of the mutual information between a pathway and PAM50 labels, obtained from models trained with different random seeds on the MO-TCGA-BRCA dataset. Models are labelled as either Aligned Consensus Model (bottom) or not (top). From left to right, we show the results for dropout p being 0.0, 0.1, 0.5, 0.8, and 0.9. Note that model order changes drastically at the last two dropout values.

Figure 19: Average top-k agreement of Concordance-Index rankings between models models trained with different random seeds. Models are labelled as either Aligned Consensus Model (right) or not (left). Shaded regions indicate model variance.

38

Figure 20: Average top-k agreement the mutual information rankings between models models trained with different random seeds. Models are labelled as either Aligned Consensus Model (right) or not (left). Shaded regions indicate model variance.

39

0.12

0

LumA 0.016 0.028

0.51

0.38

0.072

LumB 0.042

0.16

0.61

0.016

True label

Basal

0.17

0

0.034 0.025 0.017

Basal

Her2 LumA LumB Normal Predicted label

0.72

0.14 0.0048 0.014

0.18

0.8

0.6 0.4

0.12

0.18

0.11

LumA 0.0029 0.054

0.35

0.28

0.31

LumB 0.048

0.25

0.14

Normal 0.047 0.047

0.14

0.095

0.082 0.67

0.4

Normal vs the rest (AUC = 0.94) LumA vs the rest (AUC = 0.83) LumB vs the rest (AUC = 0.73) Basal vs the rest (AUC = 0.96) Her2 vs the rest (AUC = 0.79)

0.0

0.0

0.0

0.7

0.2

0.4 0.6 False Positive Rate

0.8

1.0

PAAE-[32]-[64] (KEGG) - Pathway Space DecisionTreeClassifier (AUC = 0.74 on Metabric)

0.6

1.0

0.5

0.8

0.4 0.3

0.48

0.6

0.2

0.92

0.085

Basal

1.0

0.2

0.44

Her2

PAAE-[32]-[64] (KEGG) - Pathway Space DecisionTreeClassifier (AUC = 0.85 on TCGA)

0.8

True Positive Rate

0.18

Normal

0.12

0.057 0.0071 0.0071 0.071 0.58

Her2

True label

0.86

True Positive Rate

Basal

0.2

0.6 0.4

Normal vs the rest (AUC = 0.76) LumA vs the rest (AUC = 0.69) LumB vs the rest (AUC = 0.71) Basal vs the rest (AUC = 0.92) Her2 vs the rest (AUC = 0.65)

0.2

0.1

0.0

Her2 LumA LumB Normal Predicted label

0.0

0.2

0.4 0.6 False Positive Rate

0.8

1.0

Figure 21: Confusion matrices (left) and receiver operating characteristic curves (right) for a constrained decision tree model trained only on 5 KEGG pathways with highest OvR mutual information with each of PAM50 subtypes. Results are shown for both the TCGA dataset (top, training) and Metabric (bottom, external validation).

40

TCGA Metabric

PAM50 Normal LumA LumB Basal Her2

KEGG

Top 5 MI

Figure 22: 2-dimensional UMAP reprentations of the TCGA (top) and Metabric (bottom) datasets, for the pathway activity vector using all of the KEGG pathways (left) and using the top 5 pathways with highest OvR MI (right).

41

0.0

pway_def

gex rppa mirna meth cnv mutation gex,rppa rppa,cnv rppa,mirna gex,mutation gex,cnv gex,meth meth,rppa mutation,rppa meth,mirna mutation,mirna gex,mirna cnv,mirna meth,cnv meth,mutation mutation,cnv gex,mutation,rppa gex,rppa,cnv gex,rppa,mirna mutation,rppa,cnv mutation,rppa,mirna rppa,cnv,mirna meth,rppa,cnv gex,meth,rppa gex,mutation,cnv gex,meth,mirna meth,rppa,mirna gex,cnv,mirna meth,mutation,rppa gex,meth,mutation gex,meth,cnv meth,cnv,mirna gex,mutation,mirna meth,mutation,cnv meth,mutation,mirna mutation,cnv,mirna gex,meth,rppa,cnv gex,mutation,rppa,cnv gex,mutation,rppa,mirna mutation,rppa,cnv,mirna gex,rppa,cnv,mirna meth,mutation,rppa,cnv gex,meth,rppa,mirna gex,meth,mutation,rppa meth,mutation,rppa,mirna meth,rppa,cnv,mirna gex,meth,cnv,mirna gex,meth,mutation,cnv gex,meth,mutation,mirna meth,mutation,cnv,mirna gex,mutation,cnv,mirna gex,meth,mutation,rppa,cnv gex,mutation,rppa,cnv,mirna meth,mutation,rppa,cnv,mirna gex,meth,mutation,cnv,mirna gex,meth,mutation,rppa,mirna gex,meth,rppa,cnv,mirna gex,meth,mutation,rppa,cnv,mirna

accuracy 1.0 n_omics_layers = 1 n_omics_layers = 3

42

n_omics_layers = 5

0.8

0.6

0.4 integration late-concat late-mean early

0.2

Figure 23: Accuracy of one-vs-rest logistic regression for all combinations of omics layer in early, lateconcat and late-mean integration. The gray line indicates the highest median C-index among all tested configurations (dotted line marks accuracy of 0.2, equivalent to random chance).

0.0

Figure 24: Mutual information between the PAM50 Subtypes and the result of unsupervised clustering using pway_def the K-Means algorithm for all combinations of omics layer in early, late-concat and late-mean integration. The gray line indicates the highest median C-index among all tested configurations. cnv mirna rppa gex meth mutation rppa,cnv gex,mirna gex,rppa cnv,mirna rppa,mirna mutation,mirna mutation,cnv mutation,rppa gex,cnv gex,mutation meth,rppa gex,meth meth,mirna meth,cnv meth,mutation gex,meth,mutation mutation,rppa,cnv gex,rppa,mirna rppa,cnv,mirna gex,mutation,mirna meth,cnv,mirna mutation,rppa,mirna gex,mutation,cnv gex,mutation,rppa gex,cnv,mirna mutation,cnv,mirna meth,rppa,mirna gex,rppa,cnv gex,meth,cnv gex,meth,rppa gex,meth,mirna meth,mutation,cnv meth,mutation,rppa meth,rppa,cnv meth,mutation,mirna gex,mutation,rppa,mirna gex,mutation,rppa,cnv meth,rppa,cnv,mirna mutation,rppa,cnv,mirna gex,meth,rppa,mirna gex,mutation,cnv,mirna meth,mutation,rppa,mirna gex,meth,cnv,mirna gex,rppa,cnv,mirna gex,meth,mutation,cnv meth,mutation,cnv,mirna gex,meth,mutation,rppa gex,meth,mutation,mirna gex,meth,rppa,cnv meth,mutation,rppa,cnv gex,meth,mutation,cnv,mirna gex,meth,rppa,cnv,mirna meth,mutation,rppa,cnv,mirna gex,mutation,rppa,cnv,mirna gex,meth,mutation,rppa,mirna gex,meth,mutation,rppa,cnv gex,meth,mutation,rppa,cnv,mirna

mutual_info_score 1.0 n_omics_layers = 1 n_omics_layers = 3

43

n_omics_layers = 5

0.8

0.6

0.4 integration late-mean late-concat early

0.2

1.0

n_omics_layers = 1

n_omics_layers = 3

n_omics_layers = 5

0.8

0.6 rand_score

integration late-mean late-concat early

0.4

0.0

rppa cnv mirna gex meth mutation rppa,cnv gex,rppa mutation,rppa gex,mirna cnv,mirna rppa,mirna mutation,mirna mutation,cnv gex,mutation gex,cnv meth,mirna gex,meth meth,rppa meth,cnv meth,mutation gex,meth,mutation mutation,rppa,cnv gex,rppa,mirna gex,mutation,cnv gex,mutation,rppa gex,mutation,mirna meth,rppa,mirna rppa,cnv,mirna mutation,rppa,mirna meth,cnv,mirna mutation,cnv,mirna gex,cnv,mirna gex,rppa,cnv gex,meth,rppa gex,meth,mirna gex,meth,cnv meth,rppa,cnv meth,mutation,mirna meth,mutation,rppa meth,mutation,cnv gex,mutation,rppa,cnv gex,mutation,rppa,mirna gex,mutation,cnv,mirna meth,mutation,cnv,mirna mutation,rppa,cnv,mirna meth,rppa,cnv,mirna gex,meth,mutation,cnv gex,meth,rppa,mirna gex,rppa,cnv,mirna gex,meth,cnv,mirna gex,meth,mutation,rppa meth,mutation,rppa,mirna gex,meth,rppa,cnv gex,meth,mutation,mirna meth,mutation,rppa,cnv meth,mutation,rppa,cnv,mirna gex,meth,mutation,cnv,mirna gex,meth,rppa,cnv,mirna gex,mutation,rppa,cnv,mirna gex,meth,mutation,rppa,mirna gex,meth,mutation,rppa,cnv gex,meth,mutation,rppa,cnv,mirna

0.2

integration early late-concat late-mean

0.04 0.03 0.02 0.01 0.00 0.01 0.02

0.15

0.15

0.10

0.10

0.05 0.00 0.05 0.10 0.15

gex

meth mutation rppa omics_layer

cnv

mirna

concordance_index_effect

concordance_index_effect

0.05

concordance_index_effect

Figure 25: Rand Index between the PAM50 Subtypes and the result of unsupervised clustering using the pway_def K-Means algorithm for all combinations of omics layer in early, late-concat and late-mean integration. The gray line indicates the highest median C-index among all tested configurations.

integration early late-concat late-mean gex

0.05 0.00 0.05 0.10 0.15

meth mutation rppa omics_layer

cnv

mirna

integration early late-concat late-mean gex

meth mutation rppa omics_layer

cnv

mirna

Figure 26: Marginal contribution of each omics layer to survival prediction, measured by concordance index of a Cox proportional hazards model, across the three tested integration methods: early, late-concat and late-mean. The left panel shows the median marginal contribution, the right panel displays a swarm plot of values across all combinations excluding that layer, and in the middle the same information is shown as a boxplot.

44