Explainable Multilayer Graph Neural Network for Cancer Gene Prediction
The identification of cancer genes is a critical yet challenging problem in
cancer genomics research. Existing computational methods, including deep graph
neural networks, fail to exploit the multilayered gene-gene interactions or
provide limited explanation for their predictions. These methods are restricted
to a single biological network, which cannot capture the full complexity of
tumorigenesis. Models trained on different biological networks often yield
different and even opposite cancer gene predictions, hindering their
trustworthy adaptation. Here, we introduce an Explainable Multilayer Graph
Neural Network (EMGNN) approach to identify cancer genes by leveraging multiple
genegene interaction networks and pan-cancer multi-omics data. Unlike
conventional graph learning on a single biological network, EMGNN uses a
multilayered graph neural network to learn from multiple biological networks
for accurate cancer gene prediction. Our method consistently outperforms all
existing methods, with an average 7.15% improvement in area under the
precision-recall curve (AUPR) over the current state-of-the-art method.
Importantly, EMGNN integrated multiple graphs to prioritize newly predicted
cancer genes with conflicting predictions from single biological networks. For
each prediction, EMGNN provided valuable biological insights via both
model-level feature importance explanations and molecular-level gene set
enrichment analysis. Overall, EMGNN offers a powerful new paradigm of graph
learning through modeling the multilayered topological gene relationships and
provides a valuable tool for cancer genomics research.
Michail Chatzianastasis, Michalis Vazirgiannis, Zijun Zhang
# Explainable Multilayer Graph Neural Network for Cancer Gene Prediction

###### Abstract

The identification of cancer genes is a critical, yet challenging problem in cancer genomics research. Recently, several computational methods have been developed to address this issue, including deep neural networks. However, these methods fail to exploit the multilayered gene-gene interactions and provide little to no explanation for their predictions. Results: In this study, we propose an Explainable Multilayer Graph Neural Network (EMGNN) approach to identify cancer genes, by leveraging multiple gene-gene interaction networks and multi-omics data. Compared to conventional graph learning methods, EMGNN learned complementary information in multiple graphs to accurately predict cancer genes. Our method consistently outperforms existing approaches while providing valuable biological insights into its predictions. We further release our novel cancer gene predictions and connect them with known cancer patterns, aiming to accelerate the progress of cancer research.

## 1 Introduction

Understanding the precise function and disease pathogenicity of a gene is dependent on the target gene's properties, as well as its interaction partners in a disease-specific context [1; 2; 3]. High-throughput experiments, such as whole-genome sequencing and RNA sequencing of bulk and single-cell assays, have enabled unbiased profiling of genetic and molecular properties for all genes across the genome. Experimental methods to probe both physical [4; 5] and genetic interactions [6; 7] provide valuable insights of the functional relevance between a pair of genes. Based on these data, computational methods have been developed to predict gene functions for understudied and uncharacteristic genes by combining the gene's property with its network connectivity patterns [8; 9]. However, the prediction of gene pathogenicity in disease-specific contexts is challenging. Functional assays describing the gene and its gene network are relevant to disease only to the degree to which the measured property correlates with disease physiology [10]; while our understanding of complex disease physiology is poor, even for diseases with large sample size and data modalities, such as cancer [11].

As the completeness of known cancer genes is questioned, predicting novel cancer genes remains a crucial task in cancer genomics research. These genes, which are often mutated or over-expressed in cancer cells, play a key role in the development and progression of the disease [12]. Large-scale cancer sequencing consortia projects have generated genomic, molecular and histological profiling data for a variety of cancer types, providing an information-rich resource for identifying novel cancer genes. Building on the hypothesis that different genetic and molecular modalities provide complementary information to cancer gene pathogenicity, a pioneering work EMOGI [13] innovatively modeled the multi-omics features of cancer genes in Protein-Protein interaction (PPI) networks to predict novel cancer genes. To
address the challenge of functional properties irrelevant to cancer disease physiology, EMOGI featurized each gene by a vector summarizing multi-omics data levels across various cancer types in The Cancer Genome Atlas (TCGA) [14]. EMOGI then modeled the gene-gene interactions using pre-defined generic PPI networks using a Graph Convolution Neural network (GCN) and trained using a set of high-confidence cancer- and non-cancer genes. EMOGI identified \(165\) novel cancer genes without recurrent alterations, but interact with known cancer genes.

A major limitation of EMOGI is that it didn't address the disease physiology relevance in the pre-defined graph topology and connectivity patterns. EMOGI employed six different pre-defined graphs, including genetic-focused networks such as Multinet [15], and generic protein interaction networks such as STRING-db [16]. Among EMOGI models trained on different PPI networks, we found an average standard deviation of \(25.2\%\) in unlabelled cancer gene predictions, demonstrating the predicted novel cancer genes were different when using different PPI networks. For these genes, a trustworthy adaptation of the method's output is challenging when conflicting prediction results are present. Because cancer disease physiology is complex, using a single predefined graph to represent the gene-gene relationships cannot fully capture its molecular landscape; therefore, more sophisticated, data-driven methods are needed to decipher the gene relationships in disease-specific contexts.

To alleviate this issue, we propose a novel graph learning framework, EMGNN (Explainable Multilayer Graph Neural Network), for predicting gene pathogenicity based on multiple input graphs. EMGNN maximizes the concordance of functional gene relationships with the unknown disease physiology by jointly modeling a multilayered graph structure. We evaluated the performance of EMGNN in predicting cancer genes using the same compiled datasets as EMOGI and showed that our proposed method achieves state-of-the-art performance by combining information from all six PPI networks. Furthermore, we explained EMGNN's prediction by both model-level integrated gradients and molecular-level gene pathways. By examining novel cancer genes predicted by EMGNN, we demonstrated novel biological insights by leveraging the complementary information in different types of biological networks. Overall, EMGNN provides a powerful new paradigm of graph learning through modeling the multilayered topological gene relationships. Our key contributions can be summarized as follows:

* We propose an Explainable Multilayer Graph Neural Network (EMGNN) approach to identify cancer genes by leveraging multiple protein-protein interaction networks and multi-omics data.
* Our method demonstrates superior performance compared to existing approaches as quantified by a significant increase in the AUPRC across six PPI networks. The average improvement in performance is 7.15% over the current state-of-the-art method, EMOGI.
* We identify the most important multi-omics features for the prediction of each cancer gene, as well as the most influential PPI networks, using model interpretation strategies.
* EMGNN identifies novel cancer genes by integrating multiple PPI networks, providing a unified and robust prediction for novel cancer genes discovery. Our code is publicly available on GitHub2

Footnote 2: Code: [https://github.com/zhanglab-aim/EMGNN](https://github.com/zhanglab-aim/EMGNN)

## 2 Materials and Methods

### Datasets

We trained the proposed model with six PPI Networks: CPDB [17], Multinet [15], PCNet [18], STRING-db [16], Iref [19] and its newest version Iref(2015). As node features, we used mutation, copy number, DNA methylation and gene expression data of \(29,446\) samples from TCGA [14], from \(16\) different cancer types.

### Multilayer Graph Neural Network

**GNNs.** Let a graph be denoted by \(G=(V,E),\) where \(V=\{v_{1},\ldots,v_{N}\}\) is the set of vertices and \(E\) is the set of edges. Let \(\mathbf{A}\in\mathbb{R}^{N\times N}\) denote the adjacency matrix, \(\mathbf{X}=[x_{1},\ldots,x_{N}]^{T}\in\mathbb{R}^{N\times d_{I}}\) denote the node features and \(\mathbf{Y}=[y_{1},\ldots,y_{N}]^{T}\in\mathbb{N}^{N}\) denote the label vector. Graph neural networks have been successfully applied to many graph-structured problems[20; 21], as they can effectively leverage both the network structure and node features. They typically employ a message-passing scheme, which constitutes of the two following steps. In the first step, every node aggregates the representations of its neighbors using a permutation-invariant function. In the second step, each node updates its own representation by combining the aggregated message from the neighbors with its own previous
representation,

\[m_{u}^{(l)} =\text{Aggregate}^{(l)}\left(\left\{\mathbf{h}_{v}^{(l-1)}:v\in \mathcal{N}(u)\right\}\right), \tag{1}\] \[\mathbf{h}_{u}^{(l)} =\text{Combine}^{(l)}\left(\mathbf{h}_{u}^{(l-1)},\mathbf{m}_{u} ^{(l)}\right), \tag{2}\]

where \(h_{u}^{(l)}\) represents the hidden representation of node \(u\) at the \(l^{\text{th}}\) layer of the GNN architecture. Many choices for the _Aggregate_ and _Combine_ functions have been proposed in the recent years, as they have a huge impact in the representation power of the model[22]. Among the most popular architectures, are Graph Convolution Networks (GCNs) [23], and Graph Attention Networks (GAT) [24]. In GCN, each node aggregates the feature vectors of its neighbors with fixed weights inversely proportional to the central and neighbors' node degrees,\(\mathbf{h}_{u}^{\prime}=\mathbf{W}^{\top}\sum_{v\in\mathcal{N}(u)\cup\{u \}}\frac{\mathbf{h}_{v}}{\sqrt{d_{v}d_{u}}}\), with \(\hat{d}_{i}=1+\sum_{j\in\mathcal{N}(i)}1\). In GAT, each node aggregates the messages from its neighbor using learnable weighted scores: \(\mathbf{h}_{u}^{\prime}=\alpha_{u,u}\mathbf{W}\mathbf{h}_{u}+\sum_{v\in \mathcal{N}(u)}\alpha_{u,v}\mathbf{W}\mathbf{h}_{v}\), where the attention coefficients \(\alpha_{u,v}\) are computed as

\(\alpha_{u,v}=\frac{\exp\left(\text{LeakyReLU}\left(\mathbf{u}^{\top}\left[ \mathbf{W}\mathbf{h}_{u}\right.\left[\mathbf{W}\mathbf{h}_{u}\right]\right) \right)}{\sum_{k\in\mathcal{N}(u)\cup\{u\}}\exp\left(\text{LeakyReLU}\left( \mathbf{u}^{\top}\left[\mathbf{W}\mathbf{h}_{u}\right.\left[\mathbf{W} \mathbf{h}_{u}\right]\right)\right)})\).

**Multilayer Graph Construction.** Extending graph neural networks to handle multiple networks is not a trivial task, as they are designed to operate on a single graph. Next, we describe our method, which can accurately learn node representations using graph neural networks, from multilayer graphs.

Let \(N\) be the total number of genes, each associated with a feature vector \(x_{j}\in R^{d}\). Let also \(K\) be the number of gene-gene interaction networks. We represent each graph \(G^{(i)}\) with an adjacency matrix \(A^{(i)}\in\mathbf{Z}^{N_{i}\times N_{i}}\) and feature matrix \(X^{(i)}\in R^{N_{i}\times d}\), where \(n_{i}\) is the number of genes in the \(i\)-th network. Since some genes are not presented in all the graphs, the following equation holds \(N_{i}\leq N\), \(i\in\{0,1,\dots,K-1\}\).

In the first step, for each graph \(G^{i}\) we apply a graph neural network \(f_{1}\) that performs message-passing and updates the node representation matrix \(H^{(i)}=f_{1}(X^{(i)},A^{(i)})\) of each graph \(i\in\{0,1,\dots,K-1\}\). We set \(f_{1}\) to be shared across all graphs. This design allows us to handle a variable number of graphs while keeping the number of trainable parameters fixed. Next, we construct a meta graph \(G_{meta,j}\) for each gene/node \(j\), where genes \(j\) across all graphs are connected to a meta node \(v_{j}\). We initialize the features of the meta node \(v_{j}\) with the initial features of the corresponding gene \(j\).

Figure 1: An illustration of our proposed Explainable Multilayer Graph Neural Network (EMGNN) approach. The model consists of three main steps: (1) Apply a shared GNN to update the node representation matrix of each input graph, (2) Construct a meta graph for each gene, where the same genes across all graphs are connected to a meta node, and update the representation of the meta nodes with a second GNN (Meta GNN) (3) Use a multi-layer perceptron to predict the class of each meta node.
In the next step, we apply a second GNN \(f_{2}\) to update the representation of the meta node \(v_{j}\), \(H_{meta,j}=f_{2}(X_{meta,j},A_{meta,j})\), where \(X_{meta,j}\) contains the features of gene \(j\) from all the networks and \(A_{meta,j}\) is the adjacency matrix of the meta graph \(G_{meta,j}\), \(j\in\{0,1,\dots,N\}\). We set \(f_{2}\) to be shared across all genes. Therefore, in this stage, the model combines and exchanges information between the different networks. Finally, a multi-layer perceptron \(f_{3}\) predicts the class of the meta node \(j\), \(\hat{y_{j}}=f_{3}(H_{meta,j})\). An illustration of the proposed model can be found in Figure 1.

**Experimental Details.** To ensure a fair comparison with previous work, we adopted the experimental setup used in EMOGI. Specifically, we randomly divided the data for each testing graph into a 75% training set and a 25% testing set using stratified sampling, ensuring that the proportion of known cancer and non-cancer genes in both sets was equal. Given that our model uses multiple graphs as input for each experiment, we retained the test nodes of one graph as the test set, and added 90% of the remaining nodes from the other graphs to the training set and 10% to the validation set. The model was trained for \(2000\) epochs, using the cross-entropy loss function, and the ADAM optimizer [25] with a learning rate of \(0.001\). The initial GNN had three layers with a hidden dimension of \(64\), while the meta-GNN had a single layer with a hidden dimension of \(64\). The

### Model interpretation

Captum is a tool for understanding and interpreting the decision-making process of machine learning models [26]. It offers a range of interpretability methods that allow users to analyze the predictions made by their models and understand how different input variables contribute to these predictions. We used the integrated gradient (IG) module in Captum, to assign an importance score to each input feature. IG interprets the decisions of neural networks by estimating the contribution of each input feature to the final prediction. The integrated gradient approximates the integral of gradients of the model's output with respect to the inputs along a straight line path from a specific baseline input to the current input. The baseline input is typically chosen to be a neutral or a meaningless input, such as an all-zero vector or a random noise. Formally, let \(F(x)\) be the function of a neural network, where \(x\) is the input and \(\hat{x}\) is the baseline input. The integrated gradients for input \(x\) and baseline \(x0\) along the \(i\)-th dimension is defined as: \(\texttt{IntegratedGrads}_{i}(x)=(x_{i}-\hat{x_{i}})\int_{\alpha=0}^{1}\frac{ \partial F(\hat{x}+\alpha(x-\hat{x})}{\partial x_{i}}d\alpha\), where the integral is taken along the straight line path from \(\hat{x}\) to \(x\) and \(\partial F(x)/\partial x_{i}\) is the gradient of \(F(x)\) along the \(i\)-th dimension.

However, the traditional integrated gradient method, which is designed for single input models, is not directly applicable to graph neural networks as they have two distinct inputs, namely node features and network connectivity. This necessitates the development of a modified approach for computing integrated gradients in graph neural networks that considers both inputs. To this end, we propose a decomposition of the problem into two parts: identifying the most important node features and identifying the most crucial edges in the network separately. Since we predict the class of each gene by combining all the graphs, from the meta-node representations, we apply the interpretation analysis only to the meta-nodes.

**Node feature interpretation analysis.** We analyze the contribution of node features to the predictions of the GNN by using the traditional integrated gradient method while keeping the edges in the network fixed. Specifically, we interpolate between the current node features input and a baseline input where the node features are zero: \(Attribution_{x_{i}}=(x_{i}-\hat{x_{i}})\int_{\alpha=0}^{1}\frac{\partial F(\hat{ x}+\alpha(x-\hat{x},A)}{\partial x_{i}}d\alpha\), where \(A\) are the adjacency matrices of the graphs. Since the prediction for each gene is also based on the features of surrounding genes in the graphs, we extract attribution values for the \(k\)-hop neighbor genes as well, where \(k\) is equal to the number of message-passing layers in the first GNN. Therefore, the output of the attribution method for each node \(u\), is a matrix \(\mathbf{K}^{(u)}\in\mathbb{R}^{N\times d}\). Each entry \(K_{ij}\) of the matrix, corresponds to the attribution of the feature \(j\) of node \(i\) to the target node \(u\). From this matrix, we select the row that corresponds to the feature attributions of the corresponding meta node.

**Edge feature interpretation analysis.** To analyze the contribution of edges in the meta-graph to the predictions of the GNN, we use the integrated gradient method for the edges while keeping the node features fixed. Specifically, we interpolate between the current edge input and a baseline input where the weights of the edges are zero: \(Attribution_{e_{i}}=\int_{\alpha=0}^{1}\frac{\partial F(X,A_{\alpha})}{\partial w _{e_{i}}}d\alpha\), where \(A_{\alpha}\) corresponds to the graphs with the edge weights equal to \(\alpha\). We further normalize the attribution values of each meta node by dividing them by their maximum value, resulting in a range of [0, 1] for each edge. This explanation technique allows us to understand which edges in the meta-graph are crucial for the model's decision-making process, and therefore which input PPI networks are important for each gene prediction.

### Novel cancer gene discovery

We applied the trained EMBNN model that combined all six individual PPI networks to predict novel cancer genes in the \(n=14019\) unlabeled genes. We ranked these gene by their predicted cancer gene probability for potential novel
predicted cancer genes (NPCG) in this study. For each unlabeled gene, we also applied EMOGI models trained on individual PPI networks to predict the probability of it being a cancer gene. After inner join with all EMOGI models across six individual PPI networks, we analyzed the results of \(n=6591\) unlabelled genes whose predictions were available for all models.

### Gene set enrichment analysis

To understand the biological mechanisms of EMGNN's cancer gene prediction, we employed gene set enrichment analysis (GSEA) to analyze the functional enrichment of important gene features in curated cancer pathway annotations. Specifically, to determine the importance of neighboring gene nodes, we aggregated the maximum feature importance of each node using Captum's feature explanation results. Genes with zero importance were excluded in this analysis as they did not contribute to the prediction of this target gene. We then ranked the neighboring gene nodes based on their importance, and used this ranked gene list as input for GSEA. The enrichment p-value and multiple testing corrected FDR were computed by GSEA python package [27] against cancer hallmark gene sets [28].

## 3 Results

### Overview of EMGNN framework

To model multilayered graph structures, we developed a graph neural network model EMGNN (Figure 1). The input for EMGNN is a feature vector for each gene, and multiple graphs where each graph describes gene-gene relationships by an adjacency matrix. EMGNN first updates the graph representation within each graph layer by a shared message passing operator. As different graphs have distinct connectivity patterns, node representations will be updated differently in each layer. The shared message passing operator allows the model to incorporate new graphs while keeping the model's trainable parameters fixed.

EMGNN then introduces a meta graph layer that combines the layer-wise node representations of the same genes across the graphs, with a second GNN, referred to as the Meta GNN (Figure 1). Meta GNN enables directed message-passing to combine and exchange information from the different networks to the meta nodes, which will contain the final representations of the genes. A multi-layer perceptron (MLP) takes as input the meta node representations and performs the final node classification task.

Figure 2: Test AUPRC values of EMGNN(GCN) with respect to the number of input PPI networks. Each line represents a test set of cancer and non-cancer genes held-out in a specific PPI network. When there was only one input network, the corresponding PPI network was used to train an EMGNN model; then, the remaining PPI networks were added subsequently in the order of CPDB, Multinet, PCNet, STRING-db, Iref, and Iref(2015), excluding the first PPI network.
Notably, our EMGNN model is a generalized, multilayered form of single graph GNN. In the special cases where multilayered graphs have identical adjacency matrices or only one graph is provided as inputs, EMGNN reduces to a standard single graph GNN, where the shared message passing operators are standard GNN operators, and meta GNN reduces to an identical operation. Thus, EMGNN generalizes single graph GNN by capitalizing on the complementary information stored in multiple graphs.

### Multilayered graph improves EMGNN performance

We applied EMGNN to predict cancer genes using a previously compiled dataset (Methods). Briefly, this dataset consisted of a total of 887 labeled cancer genes, 7753 non-cancer genes and 14019 unlabeled genes. Six PPI networks were binarized to keep only high-confidence edges based on previously reported criteria.

We demonstrated that the integration of multiple graphs leads to an improvement in the performance of EMGNN. Specifically, we trained EMGNN models and evaluated the testing performance with respect to different numbers of PPI networks. As shown in Figure 2, the performance increased for each of the six PPI-network derived testing datasets, as the number of input networks increased. For CPDB, Multinet, STRING-db, IRefIndex networks, the incorporation of more graphs steadily increased the performance without reaching a plateau. For PCNet, EMGNN achieved the best testing performance by combining five networks. This behavior is largely consistent with previously reported benchmarking results, which suggest that the performance scale with the network size [18].

EMGNN trained by incorporating all six graphs achieved state-of-the-art performance for all test sets (Table 1). Each test set was an independent set of held-out labeled cancer and non-cancer genes from each network and we kept the set identical to previous reports (Methods). For most test sets, EMGNN outperformed EMGGI by a margin over 5% AURPC, with the largest gain of 11.1% in performance observed in the old version Iref. The smallest gain was observed in PCNet, likely because PCNet is already an expert-assembled graph combining the information from the other five graphs [18]. Nevertheless, for PCNet test set, EMGNN combining six graphs is significantly more accurate than EMGGI using PCNet (p-value=0.012, t-test). This aligns with our hypothesis that incorporating information from multiple networks leads to enhanced predictive power for gene pathogenicity prediction.

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline
**Method** & **CPDB** & **Multinet** & **PCNet** & **STRING-db** & **Iref** & **Iref(2015)** \\ \hline Random & 0.27 & 0.18 & 0.14 & 0.24 & 0.17 & 0.28 \\
20/20+ & 0.66 & 0.62 & 0.55 & 0.67 & 0.61 & 0.65 \\ MutSigCV & 0.38 & 0.33 & 0.27 & 0.41 & 0.35 & 0.43 \\ HotNet2 diffusion & 0.62 & 0.56 & 0.48 & 0.50 & 0.45 & 0.65 \\ DeepWalk+features RF & 0.74 & 0.71 & 0.72 & 0.71 & 0.66 & 0.71 \\ PageRank & 0.59 & 0.53 & 0.54 & 0.44 & 0.42 & 0.62 \\ GCN without omics & 0.57 & 0.53 & 0.47 & 0.39 & 0.37 & 0.64 \\ DeepWalk + SVM & 0.73 & 0.51 & 0.63 & 0.52 & 0.62 & 0.66 \\ RF & 0.60 & 0.59 & 0.51 & 0.61 & 0.54 & 0.62 \\ MLP & 0.58 & 0.63 & 0.47 & 0.63 & 0.55 & 0.64 \\ EMGGI[13] & 0.74 & 0.74 & 0.68 & 0.76 & 0.67 & 0.75 \\ EMGGI[29] & 0.775\(\pm\)0.003 & 0.732\(\pm\)0.003 & 0.745\(\pm\)0.002 & 0.763\(\pm\)0.003 & 0.701\(\pm\)0.004 & 0.757\(\pm\)0.001 \\ \hline
**EMGNN(GCN)** & **0.809**\(\pm\) 0.006 & **0.854**\(\pm\) 0.007 & **0.761**\(\pm\) 0.001 & **0.856**\(\pm\) 0.002 & **0.822**\(\pm\) 0.002 & **0.800**\(\pm\) 0.010 \\
**EMGNN(GAT)** & 0.776 \(\pm\) 0.018 & 0.796 \(\pm\) 0.034 & 0.730 \(\pm\) 0.031 & 0.805 \(\pm\) 0.307 & 0.739 \(\pm\) 0.033 & 0.773 \(\pm\) 0.049 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Test AUPRC values and standard deviation across different PPI networks across five different runs.

\begin{table}
\begin{tabular}{l c c c c c c} \hline \hline
**Method** & **CPDB** & **Multinet** & **PCNet** & **STRING-db** & **Iref** & **Iref(2015)** \\ \hline
**Random Features** & 0.703 \(\pm\) 0.001 & 0.727 \(\pm\) 0.002 & 0.615 \(\pm\) 0.009 & 0.745 \(\pm\) 0.002 & 0.674 \(\pm\) 0.001 & 0.697 \(\pm\) 0.005 \\
**All-one Features** & 0.726 \(\pm\) 0.002 & 0.769 \(\pm\) 0.001 & 0.657 \(\pm\) 0.010 & 0.779 \(\pm\) 0.010 & 0.710 \(\pm\) 0.015 & 0.725 \(\pm\) 0.013 \\ \hline
**Edge Removal(0.2)** & 0.800 \(\pm\) 0.007 & 0.841 \(\pm\) 0.016 & 0.746 \(\pm\) 0.017 & 0.841 \(\pm\) 0.009 & 0.796 \(\pm\) 0.005 & 0.786 \(\pm\) 0.011 \\
**Edge Removal(0.4)** & 0.795 \(\pm\) 0.004 & 0.834 \(\pm\) 0.009 & 0.743 \(\pm\) 0.003 & 0.828 \(\pm\) 0.004 & 0.790 \(\pm\) 0.012 & 0.802 \(\pm\) 0.006 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Test AUPRC and standard deviation of EMGNN(GCN) for different input perturbations methods across three different runs.
### Evaluating the performance of different GNN architectures and graph ablations.

Next, we performed an ablation study to assess the performance of our EMGNN model using different GNN architectures and input perturbations. Specifically, we compared Graph Convolutional Network (GCN), and Graph Attention Network (GAT)[24]. Using the best identified GNN architecture, we further permuted input features by random permutations and setting all node features to one, as well as randomly removing 20% and 40% of the edges in the input graphs.

As shown in Table 1, the GNN architecture played an essential role in EMGNN testing performance. We observed that GCN is the best-performing GNN architecture in all the datasets. Our findings demonstrated that the choice of GNN architecture has a significant impact on the performance of our model. Therefore, EMGNN refers to EMGNN(GCN) throughout the paper unless specified.

We then sought to answer the two following research questions: How robust is EMGNN to node feature and graph structure perturbations? Are the node features and the graph structure both crucial for the prediction of cancer genes? Biologically, the EMGNN node features and edges are determined using high-throughput assays that inherently have measurement errors. To this end, we examined the performance of EMGNN with GCN architecture under different types of input perturbations (Table 2). Specifically, we removed the node features and instead added uninformative vectors (random or constant).

We found that EMGNN decreased in performance for both random and all-one node features, suggesting that the node features derived from TCGA consortia were informative and highly relevant to cancer physiology and cancer gene pathogenicity. For edge ablations, we randomly removed 20% and 40% edges of each PPI network. The removal of edges slightly decreased EMGNN performance. This demonstrates the connectivity pattern robustness by jointly modeling a multilayered graph topology. Overall, EMGNN effectively leveraged both node features and edges to achieve accurate predictions.

### Explaining EMGNN reveals biological insights of cancer gene pathogenicity

Explainable and trustworthy models are essential for understanding the biological mechanisms of known cancer genes and facilitating the discovery of novel cancer genes. Therefore, we applied Captum [26] to explain the node and edge attributions of EMGNN (Methods). We focused our analysis on the relative contributions from each PPI network to the known cancer gene predictions, as measured by the importance of meta edges (Figure 3A). To examine if some PPI networks were statistically more important contributors than others, we performed an ANOVA test. We observed a significant difference in contributions from the different PPI networks (P-value=\(1.2e-65\)), suggesting certain PPI networks were more informative for cancer gene prediction. Among pairwise comparisons, Iref(2015) achieved a significantly higher contribution than Iref (P-value=\(1.3e-50\), t-test), which consolidated our observation that the incorporation of other PPI networks substantially improved upon the model using only Iref. We examined the relative

Figure 3: Explanation of each PPI network’s contribution to cancer gene predictions. A) Overall distribution of of meta-edge feature importance for all known cancer genes across six PPI networks. Meta-edge feature importance was normalized to 1 (see Methods for details). B) Representative PPI network contributions in known cancer genes and novel predicted cancer genes. CSMD3 and RSPO3 are known cancer genes; COL5A1 and MSLN are novel predicted cancer genes. A contribution of zero suggests the gene was missing in the corresponding PPI network.
contributions for two known cancer genes (CSMD3 and RSPO3) and two novel predicted cancer genes (COL5A1 and MSLN; see Figure 3B). Notably, different genes were predicted as cancer genes leveraging evidence from different PPI networks. For example, Multinet contributed to CSMD3, but not for MSLN; while COL5A1 combined all six PPI networks. Thus, EMGNN successfully learned complementary information from the connectivity patterns in each layer of the multilayer graph, as shown by the similar overall contributions from individual graphs and gene-specific variations.

As EMGNN's node features were derived by multi-omic cancer datasets, we next assessed if certain types of omic data were informative to cancer gene predictions. Our model explanation results of node features indicated that Single Nucleotide Variation (SNV) features were found to be significantly less informative than other types of features, which is consistent with previous reports that CENs were more detrimental to cancer progression than SNVs [30]. In contrast, we observed that DNA methylation was significantly more important for known cancer gene prediction than other omics data (P-value<0.01 for all three pairwise t-test of other omics against DNA methylation). We further examined the node feature importance of the same four genes as in Figure 3A. The omics feature contributions were highly gene specific. In general, DNA methylation had a moderate contribution across all four genes (Figure 4B). As DNA methylations are reversible epigenetic modifications, this may suggest potential novel therapeutic targets for certain cancer genes mediated by DNA methylation [31].

### EMGNN identifies novel cancer genes by integrating multilayer graphs

To gain biological insights for precision oncology, we applied the trained EMGNN model to predict cancer genes on the unlabeled genes. We discovered a non-trivial number of unlabeled genes with a high probability of predicted novel cancer genes (n=108 predicted cancer genes with over 95% predicted cancer gene probability), consistent with the previous report that machine learning predictions can augment the completeness of cancer gene catalogs [13].

Compared to models trained using single PPI networks, EMGNN achieved an accurate and unified novel cancer gene prediction by integrating multilayer graphs. Indeed, we observed a substantial divergence among the predictions of EMOGI models trained on individual PPI networks, with an average of \(29\%\) and \(63\%\) difference between the highest and lowest predictions for the top-100 and all unlabelled nodes, respectively. Furthermore, we found an average standard deviation of 25.2% in unlabelled cancer gene predictions of EMOGI, demonstrating the predicted novel cancer genes were different when using different PPI networks.

As a case study, we analyzed the predictions of a NPCG, COL5A1 (Figure 5). For this gene, EMOGI model trained on STRINGdb predicted a non-cancer gene with high confidence, EMOGI models trained on IRefIndex, CPDB and PCNet predicted a cancer gene with high confidence, while the models trained on Multinet and IRefIndex2015 predicted a cancer gene with moderate likelihood. The fact that STRINGdb was the best performer among models trained on individual PPI networks further complicated the decision making whether COL5A1 should be considered as a cancer/non-cancer gene. This level of divergence in predictions hinders a trustworthy adaptation of model predictions in clinical and pragmatic settings. In contrast, EMGNN integrated the information from each individual PPI networks in a data-driven approach and provided more accurate, unified predictions of cancer genes (Table 1). EMGNN predicted

Figure 4: Explanations of multi-omic node features importance in cancer gene predictions. A) Overall distribution of of node feature importance grouped by omic feature types, including single-nucleotide variants (MF), DNA methylation (METH), gene expression (GE) and copy number aberrations (CNA), for known cancer genes. B) Detailed node feature importance for the four genes analyzed in Figure 3B. X-axis labels were color-coded to match the omic feature types in panel A. Individual tumor types were coded according to TCGA study abbreviations [14].
COL5A1 as a cancer gene with high confidence (Figure 5A). Importantly, we also found all individual PPI networks were contributing similarly to the final EMGNN predictions (Figure 4B).

Leveraging the explanation results of node contributions for COL5A1 (Figure 4), we further illustrated the potential biological mechanisms of COL5A1 by a gene set enrichment analysis (Methods). We discovered that three cancer hallmark gene sets, i.e. apical junction, coagulation, and complement system (part of the innate immune system), were significantly enriched in COL5A1 neighboring genes (Figure 5B). These neighboring gene's contributions were important for EMGNN's prediction of COL5A1 as cancer gene or not. For example, the apical junction cancer hallmark genes contained genes annotated to function in cell-cell adhesion among epithelial cells, many of those were enriched in the top contributors (Figure 5C). This was further supported by mouse studies where loss of COL5A1 triggered skin cancer [32], the type of cancer with a strong epithelial cell origin. Therefore, we demonstrated how molecular mechanisms of novel predicted cancer genes could be interpreted and discovered using the explainable EMGNN framework.

## 4 Discussion

The biomedical and biological domain contains a wealth of information, which is often represented and analyzed using graph structures to reveal relationships and patterns in complex data sets. Indeed, various gene interaction and protein-protein interaction networks describe the functional relationships of genes and proteins, respectively. The gene-gene relationships were often described in generic cellular contexts and/or by integrating different, heterogeneous sources of information. Therefore, a single graph often struggles to best match disease-specific conditions, and different graph construction and integration methods render distinct predictive powers [33; 34]. Substantial efforts have been devoted to develop integrated [33; 34] and tissue-specific graphs [35]. Here, we took a complementary approach and developed a new graph learning framework, EMGNN, to jointly model multilayered graphs. Applying EMGNN to predict cancer driver genes demonstrated its superior performance over previous graph neural networks trained on single graphs. We also employed model explanation techniques to assess both node and edge feature importance. Our results showed that EMGNN leveraged the complementary information from different graph layers and omics features to predict cancer genes. Importantly, we found that cancer genes that have conflicting predictions based on different single graphs, or are missed by previous state-of-the-art predictors, can be recovered effectively using EMGNN. This demonstrates the robustness of EMGNN predictions by joint modeling the multilayered graphs.

The EMGNN model can be viewed as a data-driven, gradient-enabled integration method for multiple graphs. By providing multiple PPI networks as input, EMGNN learns from the different connectivity patterns that represent complementary information to predict cancer genes. Since all PPI networks share the same type of nodes and edges, EMGNN currently integrates homogeneous, undirected graphs; however, the EMGNN framework can be extended to various types of graphs and to perform cross-data modality integration. In biology and biomedicine, hierarchical graphs

Figure 5: EMGNN predicts COL5A1 as a novel cancer gene and reveals biological insights. A) A comparison of predicted cancer gene probability from EMGNN and EMOGI models trained on single PPI networks. As a probability of 50% equalled random guessing between cancer vs non-cancer gene, the bar heights reflected the prediction confidence. B) Three cancer hallmark genes were significantly enriched in the important neighboring genes of COL5A1 as revealed by interpreting EMGNN model. C) Enrichment of apical junction cancer hallmark genes in COL5A1 neighboring genes. The neighboring genes of COL5A1 were ranked by their EMGNN node importance on the x-axis, where each blue bar represented a gene in the apical junction geneset. A strong left-shifted curve demonstrates enrichment of apical junction geneset in the top important genes to predict COL5A1 as a cancer gene.
and heterogeneous graphs are particularly prevalent, such as Gene Ontology[36]. For example, biomedical data is often organized in hierarchical levels, starting with genes and molecules, moving on to cells and tissues, and finally reaching the level of individual patients and populations. Therefore, an interesting future direction is to apply EMGNN to model multiple graphs with more heterogeneous node and edge types, and with more complex inter-graph structures.

DNA methylation and gene expression aberrations are major contributors to EMGNN's cancer gene predictions, which are considered important features when explaining its omics node features. Unlike single nucleotide variation and copy number alteration that introduced permanent mutations to DNAs, epigenetic and transcriptomic alterations of cancer genes are potentially reversible by targeted therapies. Model explanations for EMGNN revealed these molecular aberrations that may be leveraged for screening and re-purposing of drugs, especially for previously less well-characterized, novel predicted cancer genes. This highlights the importance of model explanations to gain biological and biomedical insights when developing deep learning models to predict gene pathogenicity.

In summary, we present a novel deep learning approach for the prediction of cancer genes by integrating multiple gene-gene interaction networks. By applying graph neural networks to each individual network and then combining the representations of the same genes across networks through a meta-graph, our model is able to effectively integrate information from multiple sources. We demonstrate the effectiveness of our approach through experiments on benchmark datasets, achieving state-of-the-art performance. Furthermore, the ability to interpret the model's decision-making process through the use of integrated gradients allows for a better understanding of the contribution of different multi-omic features and PPI networks. Overall, our approach presents a promising avenue for the prediction of novel cancer genes.

## Acknowledgements

We thank all members of the Zhang laboratory for helpful discussions. This work was performed at the high-performance computing resources at Cedars-Sinai Medical Center and the Simons Foundation.

## Funding

This work has been supported by an institutional commitment fund from Cedars-Sinai Medical Center to ZZ.
