Log In Sign Up

Predicting Biomedical Interactions with Higher-Order Graph Convolutional Networks

Biomedical interaction networks have incredible potential to be useful in the prediction of biologically meaningful interactions, identification of network biomarkers of disease, and the discovery of putative drug targets. Recently, graph neural networks have been proposed to effectively learn representations for biomedical entities and achieved state-of-the-art results in biomedical interaction prediction. These methods only consider information from immediate neighbors but cannot learn a general mixing of features from neighbors at various distances. In this paper, we present a higher-order graph convolutional network (HOGCN) to aggregate information from the higher-order neighborhood for biomedical interaction prediction. Specifically, HOGCN collects feature representations of neighbors at various distances and learns their linear mixing to obtain informative representations of biomedical entities. Experiments on four interaction networks, including protein-protein, drug-drug, drug-target, and gene-disease interactions, show that HOGCN achieves more accurate and calibrated predictions. HOGCN performs well on noisy, sparse interaction networks when feature representations of neighbors at various distances are considered. Moreover, a set of novel interaction predictions are validated by literature-based case studies.


page 1

page 2

page 3

page 4


MixHop: Higher-Order Graph Convolution Architectures via Sparsified Neighborhood Mixing

Existing popular methods for semi-supervised learning with Graph Neural ...

SkipGNN: Predicting Molecular Interactions with Skip-Graph Networks

Molecular interaction networks are powerful resources for the discovery....

Predicting Biomedical Interactions with Probabilistic Model Selection for Graph Neural Networks

A biological system is a complex network of heterogeneous molecular enti...

Graph Embedding on Biomedical Networks: Methods, Applications, and Evaluations

Motivation: Graph embedding learning which aims to automatically learn l...

Higher-order Weighted Graph Convolutional Networks

Graph Convolution Network (GCN) has been recognized as one of the most e...

Tri-graph Information Propagation for Polypharmacy Side Effect Prediction

The use of drug combinations often leads to polypharmacy side effects (P...

SumGNN: Multi-typed Drug Interaction Prediction via Efficient Knowledge Graph Summarization

Thanks to the increasing availability of drug-drug interactions (DDI) da...

1 Introduction

Abiological system is a complex network of various molecular entities such as genes, proteins, and other biological molecules linked together by the interactions between these entities. The complex interplay between various molecular entities can be represented as interaction networks with molecular entities as nodes and their interactions as edges. Such a representation of a biological system as a network provides a conceptual and intuitive framework to investigate and understand direct or indirect interactions between different molecular entities in a biological system. Study of such networks lead to system-level understanding of biology [8] and discovery of novel interactions including protein-protein interactions (PPIs) [20], drug-drug interactions (DDIs) [41], drug-target interactions (DTIs) [21] and gene-disease associations (GDIs) [2].

Recently, the generalization of deep learning to the network-structured data 

[6] has shown great promise across various domains such as social networks [30], recommendation systems [37], chemistry [10], citation networks [17]. These approaches are under the umbrella of graph convolutional networks (GCNs). GCNs repeatedly aggregate feature representations of immediate neighbors to learn the informative representation of the nodes for link prediction. Although GCN based methods show great success in biomedical interaction prediction [39, 41], the issue with such methods is that they only consider information from immediate neighbors. SkipGNN [13] leverages skip graph to aggregate feature representations from direct and second-order neighbors and demonstrated improvements over GCN methods in biomedical interaction prediction. However, SkipGNN cannot be applied to aggregate information from higher-order neighbors and thus fail to capture information that resides farther away from a particular interaction [1].

To address the challenge, we propose an end-to-end deep graph representation learning framework named higher-order graph convolutional networks (HOGCN) for predicting interactions between pairs of biomedical entities. HOGCN learns a representation for every biomedical entity using an interaction network structure and/or features . In particular, we define a higher-order graph convolution (HOGC) layer where the feature representations from -order neighbors are considered to obtain the representation of biomedical entities. The layer can thus learn to mix feature representations of neighbors at various distances for interaction prediction. Furthermore, we define a bilinear decoder to reconstruct the edges in the input interaction network by relying on feature representations produced by HOGC layers. The encoder-decoder approach makes HOGCN an end-to-end trainable model for interaction prediction.

We compare HOGCN’s performance with that of state-of-the-art network similarity-based methods [19], network embedding methods [27, 12], and graph convolution-based methods [16, 17, 13] for biomedical link prediction. We experiment with various interaction datasets and show that our method makes accurate and calibrated predictions. HOGCN outperforms alternative methods based on network embedding by up to 30%. Furthermore, HOGCN outperforms graph convolution-based methods by up to 6%, alluding to the benefits of aggregating information from higher-order neighbors.

We perform a case study on the DDI network and observe that aggregating information from higher-order neighborhood allows HOGCN to learn meaningful representation for drugs. Moreover, literature-based case studies illustrate that the novel predictions are supported by evidence, suggesting that HOGCN can identify interactions that are highly likely to be a true positive.

In summary, our study demonstrates the ability of HOGCN to identify potential interactions between biomedical entities and opens up the opportunities to use the biological and physicochemical properties of biomedical entities for a follow-up analysis of these interactions.

2 Related works

With the increasing coverage of the interactome, various network-based approaches have been proposed to exploit already available interactions to predict missing interactions [3, 14, 18, 19]

. These methods can be broadly classified into (1) network similarity-based methods (2) network embedding methods (3) graph convolution-based methods. We next summarize these categories of methods for biomedical interaction prediction.

Given a network of known interactions, various similarity metrics are used to measure the similarity between the biomedical entities [3] with an assumption that higher similarity indicates interaction. Triadic closure principle (TCP) has been explored in biomedical interaction prediction with the hypothesis that biomedical entities with common interaction partners are likely to interact with each other [14]. TCP relies on a common neighbor algorithm to count the number of shared neighbors between the nodes and is quantified by where

is the adjacency matrix. Recently, L3 heuristic 

[19] shows the common neighbor hypothesis fails for most protein pairs in PPI prediction and proposes to consider nodes that are similar to the neighbors of the nodes and can be quantified by . This indicates that higher-order neighbors are important for interaction prediction.

Next, network embedding methods embed the existing networks to low-dimensional space that preserves the structural proximity such that the nodes in the original network can be represented as low-dimensional vectors. Deepwalk 

[27] is a popular approach that generates the truncated random walks in the network and defines a neighborhood for each node as a set of nodes within a window of size in each random walk. Similarly, node2vec performs a biased random walk by balancing the breadth-first and depth-first search in the network. The random walks generated by these methods can be considered as a combination of nodes from various order of neighborhoods such as 1-hop to -hop neighborhood. In other words, DeepWalk and node2vec learn the embeddings for the nodes in the network from the combination of where is the

power of the adjacency matrix. These embeddings can then be fed into a classifier to predict the interaction probability. These methods are only limited to the structure of the biomedical networks and cannot incorporate additional information about the biomedical entities. Also, they cannot learn the feature difference between nodes at various distances.

Furthermore, graph convolution-based methods use a message-passing mechanism to receive and aggregate information from neighbors to generate representations for the nodes in the network. Graph convolutional networks (GCNs) [16]

and variational graph convolutional autoencoder (VGAE) 

[17] aggregate feature representation from immediate neighbors to learn the representation of biomedical entities in an end-to-end manner using link prediction objective. These methods are only limited to the average pooling of the neighborhood features [1]. SkipGNN [13] therefore proposes to use skip similarity between the biomedical entities to aggregate information from second-order neighbors. However, these methods cannot aggregate feature representations from higher-order neighbors and also cannot learn feature differences between neighbors at various distances.

3 Preliminaries

A biomedical network is defined as where denotes the set of nodes representing biomedical entities (e.g. proteins, genes, drugs, diseases) and denotes the number of nodes. denotes a set of interactions between biomedical entities. is the features of biomedical entities where is the dimension of features.

Let denote the adjacency matrix of , where indicates an edge between nodes and . We assume the case of binary adjacency matrix where represents the existence of edge between the nodes and , indicating the presence of the experimental evidence for their interaction (i.e. ) or the absence of the experimental evidence for their interaction (i.e. ). Note that the same notation of adjacency matrix can be used to represent weighted graphs such that . Table I shows the notations and their definitions used in the paper.

Notation Definition
Graph with nodes , edges and features
Test edges
Adjacency matrix of graph
Degree matrix with
Identity matrix
Symmetrically normalized adjacency matrix
-dimensional feature matrix
Ground-truth interaction between nodes and
Probability of interaction between nodes and
Final node embeddings
Weight matrix for adjacency power for layer
Number of HOGC layers

Number of training epochs

The order of neighborhood
A set of integer adjacency powers
Representation of neighbors at distance in layer
TABLE I: Terms and notations

Problem Statement. (Biomedical interaction prediction) Given a biomedical interaction network and the set of potential biomedical interactions , we aim to learn a interaction prediction model to predict the interaction probabilities of , .

Fig. 1: Block diagram of proposed HOGCN model with 2 HOGC layers. Given a biomedical interaction network with initial features for biomedical entities, the encoder mixes the feature representation of neighbors at various distances and learn final representation . The decoder takes the representation and of nodes and to learn the representation for the edge (denoted by ) and predict probability of its existence.

3.1 Message Passing

Given a biomedical network , message passing algorithms learn the representation of biomedical entities in the network by aggregating information from immediate neighbors [10]. Additional information about biomedical entities can be used to initialize the feature matrix . These algorithms involve the message passing step in which each biomedical entity sends its current representation to, and aggregates incoming messages from its immediate neighbors. Representation for each biomedical entity can be obtained after steps of message passing and feature aggregation. However, such message passing operation is limited to average pooling of features from immediate neighbors and thus is unable to learn feature differences among neighbors at various distances [1].

Neighborhood nodes at various distances provide network structure information at different granularities [4, 7, 40, 32, 31]. Taking -hop neighborhoods into consideration, we aim at aggregating information from various distances at every message passing step. Different powers of adjacency matrices such as provide information about the network structure at different scales. Higher-order message passing operations can therefore learn to mix their representations using various powers of the adjacency matrix at each message passing step.

3.2 Graph Convolutional Networks (GCNs)

Graph convolutional networks (GCNs) are the generalization of convolution operation from regular grids such as images or texts to graph structured data [6, 36]. The key idea of GCNs is to learn the function to generate the node’s representation by repeatedly aggregating information from immediate neighbors. The graph convolutional layer is defined as:


where and are the input and output activations, is a trainable weight matrix of the layer , is the element-wise activation, and is a symmetrically normalized adjacency matrix with self-connections . A GCN model with layers is then defined as:

and can be used to predict the probability of interactions between biomedical entities.

4 Higher-order Graph Convolution Network (HOGCN)

In this work, we develop a higher-order graph convolutional network (HOGCN) that takes an interaction network as input and reconstruct the edges in the interaction network (Fig. 1). HOGCN has two main components:

  • Encoder: a higher-order graph convolution encoder that operates on an interaction graph and produces representations for biomedical entities by aggregating features from the neighborhood at various distances and

  • Decoder: a bilinear decoder that relies on these representations to reconstruct the interactions in .

4.1 Higher-Order Graph Encoder

We first describe the higher-order graph encoder, that operates on an undirected interaction graph and learns the representations for biomedical entities.

We develop an encoder with higher-order Graph Convolution (HOGC) layer to mix feature representations from neighbors at -distances. Specifically, HOGC layer considers the neighborhood information at different granularities and is defined as:


where is a set of integer adjacency powers, denotes the adjacency matrix multiplied times, and denotes column-wise concatenation [1]. Graph convolutional network [16] only considers the power of adjacency matrix and can be exactly recovered by setting in Equation (2). Similarly, SkipGNN [13] considers direct and skip similarity and can be recovered by setting in Equation (2).

Fig. 2: High-order graph convolution (HOGC) Layer with . The feature representation is a linear combination of the neighbors at multiple distances . represents feature representation of neighbors at distances for layer .

Fig. 2 shows a HOGC layer with where is maximum order of neighborhood considered in each HOGC layer. If , HOGC layer only considers the features of the biomedical entities and can capture the feature similarity between various biological entities. This is equivalent to a fully connected network with features of biomedical entities as input. For the HOGC layer, is the identity matrix where is the number of nodes in the network. This allows the HOGC layer to learn the transformation of node features separately and mix it with feature representations from neighbors.

The maximum order of neighborhood and the number of trainable weight matrices , one per each adjacency power, can vary across layers. However, we set the same for neighborhood aggregation and the same dimension for all the weight matrices across all layers.

Neighborhood features from different adjacency powers at layer are column-wise concatenated to obtain feature representation . As shown in Fig 2, weight at layer can learn the arbitrary linear combination of the concatenated features to obtain . Specifically, the layer can assign different coefficients to different columns in the concatenated matrix. For instance, the layer can assign positive coefficients to the columns produced by certain power of and assign negative coefficients to other powers. This allows the model to learn feature differences among neighbors at various distances. We apply HOGC layers to learn the latent representation for biomedical entities in the network, where and is the dimension of node’s representation for each adjacency power.

4.2 Interaction Decoder

We introduced the encoder based on HOGC layers that learns feature representation for biomedical entities by mixing neighborhood information at multiple distances. Next, we discuss the decoder that reconstructs the interactions in based on the representation .

We adopt a bilinear layer to fuse the representation of biomedical entities and and learn the edge representation . More precisely, we define a simple bilinear layer that takes the representation and as input:


where represents the learnable fusion matrix, is the representation of edge between nodes and , denotes the bias of the bilinear layer. is non-linearity.

The edge representation is then fed into 2-layered fully connected neural network to predict probability for edge :


where denotes fully connected layer with weight and bias , represents the probability that biomedical entities and interact.

So far, we have discussed the encoder and decoder of our proposed approach. Next, we describe the training procedure of our proposed HOGCN model. In particular, we explain how to optimize the trainable neural network weights in an end-to-end manner.

4.3 Training HOGCN

During HOGCN training, we employ binary cross entropy loss to optimize the model parameters


and encourage the model to assign higher probability to observed interactions than to randomly selected non-interactions. is the predicted interaction probability between and and

denotes the ground-truth interaction label between these nodes. The final loss function considering all interactions is


We follow an end-to-end approach to jointly optimize over all trainable parameters and backpropagate the gradients through encoder and decoder of HOGCN.

4.3.1 Algorithm

HOGCN leverages biomedical network structure along with additional information about biomedical entities as the initial feature representation . In this paper, we initialize the initial features

to be one-hot encoding

i.e. . The feature matrix can be initialized with properties of biomedical entities or pre-trained embeddings from other network-based approaches such as DeepWalk, node2vec.

1:  Inputs: , ,
3:  for  to  do
4:     Sample mini-batch of training edges and their interaction labels
5:     for  to  do
7:        for  to  do
10:        end for
12:     end for
15:     Compute loss in (6)
16:     Update model parameters via gradient descent
17:  end for
Algorithm 1 Training of HOGCN for biomedical interaction prediction

Given an adjacency matrix and the initial node representations , higher-order neighborhood indicated by the higher power of the adjacency matrix is iteratively computed that makes the model more efficient. By adopting right-to-left multiplication, for instance, we can calculate as (Line in Algorithm 1). Representation learned for the neighborhood at distances are concatenated to obtain the representation as shown in Fig. 2 (Line 11 in Algorithm 1). After passing through HOGC layers, we obtain the final representation for biomedical entities. With the final representations and the mini-batch of training edges, we retrieve the embeddings for the nodes in training edges and feed them into the interaction decoder to compute their interaction probabilities.

The parameters of HOGCN are optimized with a binary cross-entropy loss (Equation (6)) in an end-to-end manner. Given two biomedical entities and , the trained model can predict the probability of their interactions.

5 Experimental design

We view the problem of biomedical interaction prediction as solving a link prediction task on an interaction network. We consider various interaction datasets and compare our proposed method with the state-of-the-art methods.

5.1 Datasets

We conduct interaction prediction experiments on four publicly-available biomedical network datasets:

  • BioSNAP-DTI [42]: DTI network contains 15,139 drug-target interactions between 5,018 drugs and 2,325 proteins.

  • BioSNAP-DDI [42]: DDI network contains 48,514 drug-drug interactions between 1,514 drugs extracted from drug labels and biomedical literature.

  • HuRI-PPI [20]: HI-III human PPI network contains 5,604 proteins and 23,322 interactions generated by multple orthogonal high-throughput yeast two-hybrid screens.

  • DisGeNET-GDI [29]: GDI network consists of 81,746 interactions between 9,413 genes and 10,370 diseases curated from GWAS studies, animal models and scientific literature.

Table II provides summary of datasets used in our experiments. We provide the number of interactions used for training, validation, and testing for each interaction datasets. Also, the table includes the average number of interactions for each biomedical entity which can be computed as .

Dataset # Nodes # Edges Avg. node degree
Training (70%) Validation (10%) Testing (20%) Total
DTI 5,018 drugs, 2,325 proteins 10,597 1,514 3,028 15,139 4.12
DDI 1,514 drugs 33,960 4,852 9,702 48,514 64.09
PPI 5,604 proteins 16,326 2,332 4,664 23,322 8.32
GDI 9,413 genes, 10,370 diseases 57,222 8,175 16,349 81,746 8.26
TABLE II: Summary of the datasets used in our experiments.

5.2 Baselines

We compare our proposed model with the following network-based baselines for interaction prediction:

  • network similarity-based methods

    • [leftmargin=*]

    • L3 [19] counts the number of paths with length-3 normalized by the degree for all the node pairs.

  • Network embedding methods

    • [leftmargin=*]

    • DeepWalk [27] performs truncated random walk exploring the network neighborhood of nodes and applies skip-gram model to learn the

      -dimensional embedding for each node in the network. Node features are concatenated to form edge representation and train a logistic regression classifier.

    • node2vec [12] extends DeepWalk by running biased random walk based on breadth/depth-first search to capture both local and global network structure.

  • Graph convolution-based methods

    • [leftmargin=*]

    • VGAE [17] uses graph convolutional encoder with two GCN layers to learn representation for each node in the network and adopts inner product decoder to reconstruct adjacency matrix.

    • GCN [16] uses normalized adjacency matrix to learn node representations. The representation for nodes are concatenated to form feature representation for the edges and fully connected layer use these edge representation to reconstruct edges, similar to HOGCN. Setting in our proposed HOGCN is equivalent to GCN.

    • SkipGNN [13] learns the node embeddings by combining direct and skip similarity between nodes. Setting in our proposed HOGCN is equivalent to SkipGNN.

5.3 Experimental setup

We split the interaction dataset into training, validation, and testing interactions in a ratio of 7:1:2 as shown in Table II

. Since the available interactions are positive samples, the negative samples are generated by randomly sampling from the complement set of positive examples. Five independent runs of the experiments with different random splits of the dataset are conducted to report the prediction performance. We use (1) area under the precision-recall curve (AUPRC) and (2) area under the receiver operating characteristics (AUROC) as the evaluation metrics. With these evaluation metrics, we expect positive interactions to have higher interaction probability compared to negative interactions. So, the higher value of AUPRC and AUROC indicates better performance.

We implement HOGCN using PyTorch 

[26] and perform all experiments on a single NVIDIA GeForce RTX 2080Ti GPU. We construct a -layered HOGC network with for each layer. At each HOGC layer, the node mixes the feature representations from neighbors at distances and . The dimension of all weight matrices in HOGC layers is set to . All the weight matrices are initialized using Xavier initialization [11]. We train our model using mini-batch gradient descent with Adam optimizer [15] for a maximum of 50 epochs, with a fixed learning rate of . We set the mini-batch size to 256 and the dropout probability [33] to 0.1 for all layers. Early stopping is adopted to stop training if validation performance does not improve for 10 epochs. The dimension of the edge feature from the bilinear layer is followed by linear layers to project the edge features to edge probabilities. For baseline methods, we follow the same experimental settings discussed in [13].

6 Results

In this section, we investigate the performance and flexibility of HOGCN on interaction prediction using four different datasets. We further explore the robustness of HOGCN to sparse networks. Finally, we demonstrate the ability of HOGCN to make novel predictions with literature-based case studies.

6.1 Biomedical interaction prediction

We compare HOGCN against various baselines on biomedical interaction prediction tasks using four different types of interaction datasets including protein-protein interactions (PPIs), drug-target interactions (DTIs), drug-drug interactions (DDIs) and gene-disease associations (GDIs).

We randomly mask of interactions from the network as a test set and as a validation set. We train all models with

of interactions and evaluate their performances on test sets. The best set of hyperparameters is selected based on their performances on the validation dataset. Finally, the experiment is repeated for five independent random splits of the interaction dataset and the results with

one standard deviation are reported in Table 

III. All of our models used for the reported results are of same capacity (i.e. and ).

Dataset Method AUPRC AUROC
DTI DeepWalk 0.753 0.008 0.735 0.009
node2vec 0.771 0.005 0.720 0.010
L3 0.891 0.004 0.793 0.006
VGAE 0.853 0.010 0.800 0.010
GCN 0.904 0.011 0.899 0.010
SkipGNN 0.928 0.006 0.922 0.004
2-4 HOGCN 0.937 0.001 0.934 0.001
DDI DeepWalk 0.698 0.012 0.712 0.009
node2vec 0.801 0.004 0.809 0.002
L3 0.860 0.004 0.869 0.003
VGAE 0.844 0.076 0.878 0.008
GCN 0.856 0.005 0.875 0.004
SkipGNN 0.866 0.006 0.886 0.003
2-4 HOGCN 0.897 0.003 0.911 0.002
PPI DeepWalk 0.715 0.008 0.706 0.005
node2vec 0.773 0.010 0.766 0.005
L3 0.899 0.003 0.861 0.003
VGAE 0.875 0.004 0.844 0.006
GCN 0.909 0.002 0.907 0.006
SkipGNN 0.921 0.003 0.917 0.004
2-4 HOGCN 0.930 0.002 0.922 0.001
GDI DeepWalk 0.827 0.007 0.832 0.003
node2vec 0.828 0.006 0.834 0.003
L3 0.899 0.001 0.832 0.001
VGAE 0.902 0.006 0.873 0.009
GCN 0.909 0.002 0.906 0.006
SkipGNN 0.915 0.003 0.912 0.004
2-4 HOGCN 0.941 0.001 0.936 0.001
TABLE III: Average AUPRC and AUROC with one standard deviation on biomedical interaction prediction

Table III shows that HOGCN achieves huge improvement over network embedding methods such as DeepWalk and node2vec across all datasets. Specifically, HOGCN outperforms Deepwalk on AUPRC by 24.44% in DTI, 28.51% in DDI, 30.07% in PPI, and 13.79% in GDI. Although node2vec achieves better performance compared to DeepWalk by adopting a biased random walk, HOGCN still outperforms node2vec by a significant margin. DeepWalk and node2vec consider different orders of neighborhood defined by the window size and learns similar representations for the nodes in that window. In contrast, HOGCN learns feature differences between neighbors at various distances to obtain feature representation for the node and thus achieves superior performance. The improved performance suggests that feature differences between different order neighbors provide important information for interaction prediction.

A network similarity-based method, L3 [19] outperforms network embedding methods across four datasets but is limited to a single aspect of network similarity i.e. the number of paths of length connecting two nodes. So, L3 cannot be applied when other similarities between nodes such as similarity in features and common neighbors at various distances need to be considered. HOGCN overcomes these limitations and outperforms L3 across all interaction datasets with huge gain. In particular, HOGCN gains 3.5% AUPRC and 7.09% AUROC on PPI over L3 [19], which recently outperformed 20 network science methods in the PPI prediction problem.

Fig. 3: Visualization of learned representation for drugs with (a) GCN (b) SkipGNN (c) HOGCN. Drugs are mapped to the 2D space using t-SNE package [22] with learned drug representations. Drugs categories such as DBCAT002175, DBCAT000702 and DBCAT000592 are highlighted. The number of drugs in each categroy is reported in legend. Best viewed on screen.

Graph convolution-based methods such as GCN and VGAE achieves significant improvements over network embedding approaches but achieves comparable performance with L3. SkipGNN shows improvement over all other methods by incorporating skip similarity to aggregate information from second-order neighbors. Moreover, HOGCN with achieves an improvement over all graph convolution-based methods. Specifically, HOGCN achieves improvement in AUPRC over VGAE [17] by 6.3%, GCN [16] by 4.8% and SkipGNN [13] by 3.6% on DDI dataset.

As HOGCN can learn the linear combination of node features at multiple distances, it can extract meaningful representations from the interaction networks. The results in Table III demonstrate that our approach with higher-order neighborhood mixing outperforms the state-of-the-art methods on real interaction datasets.

6.2 Exploration of HOGCN’s drug representations

Next, we evaluate if HOGCN learns meaningful representation when feature representations of higher-order neighbors are aggregated. To this aim, we train GCN, SkipGNN, and HOGCN models on the DDI network to obtain the drug representations . The learned drug representations are mapped to 2D space using t-SNE [22] and visualize them in Fig. 3.

Drugbank [35] provides information about drugs and their categories based on different characteristics such as involved metabolic enzymes, class of drugs, side effects of drugs, and the like. For this experiment, we collect drug categories from Drugbank and limit the selection of drug categories such that the training set doesn’t contain any interactions between the drugs in the same category. The selected drug categories are ACE Inhibitors and Diuretics (DBCAT002175), Penicillins (DBCAT000702), and Antineoplastic Agents (DBCAT000592) with 10, 24, and 16 drugs respectively. Although these drugs don’t have direct interactions in the training set, we assume that these drugs share neighborhoods at various distances and can be explored accordingly with HOGCN.

Fig. 3 shows the clustering structure in drugs’ representations as neighborhood information at multiple distances are considered. Examining the figure, we observe that drugs in the same category are embedded close to each other in the 2D space when the model aggregates information from farther neighbors. For example, 24 drugs in the Penicillins (DBCAT000702) category (marked with blue triangles in Fig. 3) are scattered in the representation space learned by GCN that only considers feature aggregation from immediate neighbors (Fig. (a)a). Note that these drugs don’t have any direct interaction between themselves in the training set. Since GCN-based models can only average the representation from immediate neighbors, these drugs are mapped relatively farther to each other and closer to other interacting drugs. SkipGNN considers skip similarity to aggregate features from second-order neighbors and show relatively compact clusters compared to GCNs (Fig. (b)b). On the other hand, HOGCN considers the higher-order neighborhood and learns similar representations for drugs that belong to the same category demonstrated by compact clustering structure in Fig. (c)c even though no information about categorical similarity is provided to the model. This analysis demonstrates that HOGCN learns meaningful representation for drugs by aggregating feature representations from the neighborhood at various distances.

Next, we test if the clustering pattern in Fig. 3 holds across many drug categories. With this aim, we consider all drug categories in DrugBank and compute the average Euclidean distance between each drug’s representation and representations of other drugs within the same drug category. We then perform 2-sample Kolmogorov–Smirnov tests and found that HOGCN learns significantly more similar representations of drugs than expected by chance (p-value = ), GCNs (p-value = ) and SkipGNN (). Thus, this analysis indicates that HOGCN learns meaningful representations for drugs by aggregating neighborhood information at various distances.

6.3 Robustness to network sparsity

We next explore the robustness of network-based interaction prediction models to network sparsity. To this aim, we evaluate the performance with respect to the percentage of training edges varying from 10% to 70%. We make predictions on the rest of the interactions. We further use 10% of test edges for validation to select the best hyperparameter settings. For a fair comparison, we compare with graph convolution-based methods that aggregate information from direct and/or second-order neighbors.

Fig. 4: AUPRC comparison of HOGCN’s performance with that of alternative approaches with respect to network sparsity. HOGCN consistently achieves better performance in various fraction of training edges.

Fig. 4 shows the robustness of HOGCN to network sparsity. HOGCN achieves strong performance in all tasks with different network sparsity. The performance of HOGCN steadily improves with the increase in training edges. The mixing of features from a higher-order neighborhood in HOGCN and SkipGNN shows improvement over GCN and VGAE that only consider direct neighbors. Since HOGCN can learn the linear combination of features from a 3-hop neighborhood for this experiment, it shows improvement over SkipGNN in almost all cases. This demonstrates that features from farther distances are informative for interaction prediction in sparse networks.

6.4 Calibrating model’s prediction

All graph convolution-based model proposed for biomedical link prediction predicts the confidence estimate

for interaction between two biomedical entities and . We thus test if a predicted confidence represents the likelihood of being true interaction. In other words, we expect the confidence estimate to be calibrated, i.e. represents true interaction probability [24]. For example, given 100 predictions, each with the confidence of 0.9, we expect that 90 interactions should be correctly classified as true interactions.

Fig. 5: Reliability diagrams for different graph convolution-based methods. The calibration performance is evaluated with Brier score, reported in the legend (lower is better).

To evaluate the calibration performance, we use reliability diagrams [24] and Brier score [5]. In particular, the reliability diagram provides a visual representation of model calibration. These diagrams plot the expected fraction of positives as a function of predicted probability [24]. A model with perfectly calibrated predictions is represented by a diagonal in Fig. 5. In addition to reliability diagrams, it is more convenient to have a scalar summary statistics of calibration. Brier score [5] is a proper scoring rule for measuring the accuracy of predicted probabilities. Lower Brier score indicates better calibration of a set of predictions. It is computed as the mean squared error of a predicted probability and the ground-truth interaction label . Mathematically, the Brier score can be computed as:


where denotes the number of test edges.

Fig. 5 shows the calibration plots for GCN, SkipGNN and HOGCN (). For DTI dataset, SkipGNN show better calibration compared to GCN and HOGCN (Fig. (a)a), indicating that second-order neighborhood information is appropriate and aggregating features from farther away makes model overconfident. For other datasets, GCNs are relatively overconfident for all predicted confidence. For example, approximately of interactions are true positives among the interactions with high predicted confidence in PPI (Fig. (c)c) and GDI dataset (Fig. (d)d). In contrast, HOGCN achieves a lower Brier score in comparison to the GCN and SkipGNN across DDI, PPI, and GDI datasets, alluding to the benefits of aggregating higher-order neighborhood features for calibrated prediction. This analysis demonstrates that HOGCN with higher-order neighborhood mixing makes accurate and calibrated predictions for biomedical interaction.

6.5 Impact of higher-order neighborhood mixing

In Section 6.3, we contrast HOGCN’s performance with that of alternative graph convolution-based methods in varying fraction of edges. In this experiment, we aim to observe the performance of HOGCN when the order is increased to allow the model to aggregate neighborhood information from farther away. We follow a similar setup as discussed in 6.3.

Fig. 6: AUPRC comparison for higher-order message passing with different fractions of training edges. The values of for different HOGCN models are reported in the legend.

Fig. 6 shows the comparison of HOGCN with higher-order neighborhood mixing . The prediction performance of HOGCN improves with the increase in the number of training interactions for all cases. The results show that HOGCN’s performances are not sensitive to the hyperparameter settings of for all datasets since for settings , we achieve comparable performances across the datasets. This analysis indicates that the 3-hop neighborhood provides sufficient information for interaction prediction across all datasets and the performance remains stable even with a large value for .

6.6 Investigation of novel predictions

Next, we perform the literature-based validation of novel predictions. Our goal is to evaluate the quality of HOGCN’s novel predictions compared to that of GCN and SkipGNN and show that HOGCN predicts novel interactions with higher confidence. We consider GDIs and DDIs for this evaluation.

We first evaluate the potential of the HOGCN to make novel GDI predictions. We collect 1,134,942 GDIs and their scores from DisGeNET [29]. The score corresponds to the number and types of sources and the number of publications supporting the associations. With the score threshold of , we obtain 17,893 new GDIs that are not in the training set. We make predictions on these 17,893 GDIs with GCN, SkipGNN, and HOGCN (). Out of 17,893 GDIs, HOGCN predicts a higher probability than GCN for 17,356 (96.99%) GDIs and than SkipGNN for 11,418 (63.8%) GDIs. Table IV shows the top 5 GDIs with a significant increase in interaction probabilities when higher-order neighborhood mixing is considered. We also provide the number of evidence from DisGenNet [29] to support these predictions. Improvement in predicted probabilities by HOGCN models shows that aggregating feature representations from higher-order neighbors make HOGCN more confident about the potential interactions as discussed in Section 6.4.

Probability No. of
Gene Disease Evidence
PTGER1 Gastric ulcer 0.087 0.519 0.721 1
ANGPT1 Gastric ulcer 0.173 0.583 0.657 2
ABO Pancreatic carcinoma 0.265 0.615 0.770 26
VCAM1 Endotoxemia 0.294 0.529 0.639 5
GPC3 Hepatoblastoma 0.307 0.540 0.598 17
TABLE IV: Novel prediction of GDIs with the number of evidence from DisGenNet [29] supporting the interaction. GCN, SkipGNN and HOGCN are denoted by and respectively.

We select two predicted GDIs with a large number of supporting evidence and investigate the reason for the improvement in predicted confidence with HOGCN. Specifically, we choose gene-disease pairs (a) ABO and Pancreatic carcinoma (26 pieces of evidence) and (b) GPC3 and Hepatoblastoma ((17 pieces of evidence). To explain the prediction, the subnetworks containing all shortest paths between these pairs are selected. In particular, there are 49 shortest paths of length 3 between ABO and Pancreatic carcinoma including 6 diseases and 15 genes (Fig. (a)a). Similarly, there are 20 shortest paths of length 3 between GPC3 and Hepatoblastoma including 6 diseases and 9 genes (Fig. (b)b). Since these nodes are 3-hop away from each other and GCNs can only consider immediate neighbors, GCNs assign low confidence to these interactions.

Examining the subnetwork in Fig. (a)a, we found that most of the diseases are related to a cancerous tumor in the pancreas and the prostate. Furthermore, pancreatic carcinoma is associated with other diseases such as Pancreatic neoplasm, malignant neoplasm of pancreas, and malignant neoplasm of prostate [29]. Since ABO is linked with diseases that are related to pancreatic carcinoma and other genes are related to these diseases as well, HOGCN captures such association (Fig (a)a) even though they are farther away in the network. Similarly, HOGCN predicts association for GPC3 and Hepatoblastoma (Fig. (b)b).

Fig. 7: Subnetwork with predicted interactions (marked by bold dashed lines) between (a) ABO and Pancreatic carcinoma (b) GPC3 and Hepatoblastoma and all shortest paths between these pairs. The known interactions are presented as gray lines. Diseases are presented as dark circles and genes are presented as white circles.

Next, we perform a similar case study for DDIs and evaluate the predictions against DrugBank [35]. For this experiment, we make a prediction for every drug pair with GCN, SkipGNN and HOGCN and exclude the interactions that are already in the training set. Table V shows the top 5 interactions with an increase in interaction probabilities when highers orders of the neighborhood are considered. As discussed in Section 6.4, HOGCN makes predictions with higher confidence compared to GCN and SkipGNN for the interactions that are likely to be a true positive.

Drug 1 Drug 2 Probability
Nelfinavir Acenocoumarol 0.192 0.318 0.417 [9]
Praziquantel Itraconazole 0.609 0.721 0.811 [28]
Cisapride Droperidol 0.618 0.725 0.823 [23]
Dapsone Warfarin 0.632 0.720 0.885 [34]
Levofloxacin Tobramycin 0.663 0.760 0.823 [25]
TABLE V: Novel prediction of DDIs with the literature evidence supporting the interaction. GCN, SkipGNN and HOGCN are denoted by and respectively.

Moreover, we validate the false positive DDI predictions of GCNs and investigate the subnetwork for these drugs in DDI networks to reason the predictions. Table VI shows the top 5 interactions with a significant decrease in predicted confidence compared to GCN-based models. Since these DDIs are false positives [35], GCN-based models make overconfident predictions for such DDIs. In contrast, HOGCN significantly reduces the predicted confidence for these DDIs to be true positive, indicating that the higher-order neighborhood allows HOGCN to identify false positive predictions. In particular, HOGCN can identity false positive DDI between Belimumab and Estazolam even though they are 3-hop away from each other.

Drug 1 Drug 2 Probability with
Tranylcypromine Melphalan 0.925 0.478 0.065
Belimumab Estazolam 0.912 0.477 0.178
Methotrimeprazine Cloxacillin 0.907 0.406 0.065
Hydrocodone Melphalan 0.905 0.193 0.012
Ibrutinib Mecamylamine 0.899 0.398 0.353
TABLE VI: Predicted probability for negative DDIs. GCN, SkipGNN and HOGCN are denoted by and respectively.

We select subnetwork involving the drugs to investigate the reason for such predictions. Fig. 8 shows the subnetwork with all shortest paths between the drugs in Table VI. Examining the figure, we observe that the drugs in these false positive DDIs have common immediate neighbors for all cases. GCN makes wrong predictions for these DDIs with high confidence. However, SkipGNN becomes less confident about the interaction being true positive by considering the skip similarity. HOGCN further reduces the predicted confidence for Tranylcypromine and Melphalan to 0.065, indicating that there is no association between these drugs.

Fig. 8: Subnetwork containing false positive predictions (marked by dark dashed lines) and all shortest paths between (a) Tranylcypromine and Melphalan (b) Methotrimeprazine and Cloxacillin (c) Hydrocodone and Melphalan and (d) Ibrutinib and Mecamylamine. Other known interactions are presented as gray gray lines. Dark circles denotes drugs.

These case studies show that HOGCN with higher-order neighborhood mixing not only provide information for the identification of novel interactions but also help HOGCN to reduce false positive predictions.

7 Conclusion

We present a novel deep graph convolutional network (HOGCN) for biomedical interaction prediction. Our proposed model adopts a higher-order graph convolutional layer to learn to mix the feature representation of neighbors at various scales. Experimental results on four interaction datasets demonstrate the superior and robust performance of the proposed model. Furthermore, we show that HOGCN makes accurate and calibrated predictions by considering higher-order neighborhood information.

There are several directions for future study. Our approach only considers the known interactions to flag potential interactions. There are other sources of biomedical information such as various physicochemical and biological properties of biomedical entities that can provide additional information about the interaction and we plan to investigate the integration of such features into the model. As HOGCN aggregates the neighborhood information at various distances and can flag novel interactions, it would be interesting to provide interpretable explanations for the predictions in the form of a small subgraph of the input interaction network that are most influential for the predictions [38].


This work was supported by the NSF [NSF-1062422 to A.H.], [NSF-1850492 to R.L.] and the NIH [GM116102 to F.C.]


  • [1] S. Abu-El-Haija, B. Perozzi, A. Kapoor, N. Alipourfard, K. Lerman, H. Harutyunyan, G. Ver Steeg, and A. Galstyan (2019) MixHop: higher-order graph convolutional architectures via sparsified neighborhood mixing. In

    International Conference on Machine Learning

    pp. 21–29. Cited by: §1, §2, §3.1, §4.1.
  • [2] M. Agrawal, M. Zitnik, J. Leskovec, et al. (2018) Large-scale analysis of disease pathways in the human interactome.. In PSB, pp. 111–122. Cited by: §1.
  • [3] I. Ahmad, M. U. Akhtar, S. Noor, and A. Shahnaz (2020) Missing link prediction using common neighbor and centrality based parameterized algorithm. Scientific Reports 10 (1), pp. 1–9. Cited by: §2, §2.
  • [4] A. R. Benson, R. Abebe, M. T. Schaub, A. Jadbabaie, and J. Kleinberg (2018) Simplicial closure and higher-order link prediction. Proceedings of the National Academy of Sciences 115 (48), pp. E11221–E11230. Cited by: §3.1.
  • [5] G. W. Brier (1950) Verification of forecasts expressed in terms of probability. Monthly weather review 78 (1), pp. 1–3. Cited by: §6.4.
  • [6] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst (2017) Geometric deep learning: going beyond euclidean data. IEEE Signal Processing Magazine 34 (4), pp. 18–42. Cited by: §1, §3.2.
  • [7] S. Cao, W. Lu, and Q. Xu (2015) Grarep: learning graph representations with global structural information. In Proceedings of the 24th ACM international on conference on information and knowledge management, pp. 891–900. Cited by: §3.1.
  • [8] L. Cowen, T. Ideker, B. J. Raphael, and R. Sharan (2017) Network propagation: a universal amplifier of genetic associations. Nature Reviews Genetics 18 (9), pp. 551. Cited by: §1.
  • [9] B. Garcia, P. De Juana, T. Bermejo, J. Gómez, P. Rondon, I. Garcia Plaza, et al. (1999) Sequential interaction of ritonavir and nelfinavir with acenocoumarol (abstract 1069). In Seventh European Conference on Clinical Aspects and Treatment of HIV Infection. Lisbon, Portugal, Cited by: TABLE V.
  • [10] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl (2017) Neural message passing for quantum chemistry. In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pp. 1263–1272. Cited by: §1, §3.1.
  • [11] X. Glorot and Y. Bengio (2010) Understanding the difficulty of training deep feedforward neural networks. In

    Proceedings of the thirteenth international conference on artificial intelligence and statistics

    pp. 249–256. Cited by: §5.3.
  • [12] A. Grover and J. Leskovec (2016) Node2vec: scalable feature learning for networks. In Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 855–864. Cited by: §1, 2nd item.
  • [13] K. Huang, C. Xiao, L. Glass, M. Zitnik, and J. Sun (2020) SkipGNN: predicting molecular interactions with skip-graph networks. arXiv preprint arXiv:2004.14949. Cited by: §1, §1, §2, §4.1, 3rd item, §5.3, §6.1.
  • [14] O. Keskin, N. Tuncbag, and A. Gursoy (2016) Predicting protein–protein interactions from the molecular to the proteome level. Chemical reviews 116 (8), pp. 4884–4909. Cited by: §2, §2.
  • [15] D. P. Kingma and J. Ba (2015) Adam (2014), a method for stochastic optimization. In Proceedings of the 3rd International Conference on Learning Representations (ICLR), arXiv preprint arXiv, Vol. 1412. Cited by: §5.3.
  • [16] T. N. Kipf and M. Welling (2016) Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:1609.02907. Cited by: §1, §2, §4.1, 2nd item, §6.1.
  • [17] T. N. Kipf and M. Welling (2016) Variational graph auto-encoders. arXiv preprint arXiv:1611.07308. Cited by: §1, §1, §2, 1st item, §6.1.
  • [18] K. Kishan, R. Li, F. Cui, Q. Yu, and A. R. Haake (2019) GNE: a deep learning framework for gene network inference by aggregating biological information. BMC systems biology 13 (2), pp. 38. Cited by: §2.
  • [19] I. A. Kovács, K. Luck, K. Spirohn, Y. Wang, C. Pollis, S. Schlabach, W. Bian, D. Kim, N. Kishore, T. Hao, et al. (2019) Network-based prediction of protein interactions. Nature communications 10 (1), pp. 1–8. Cited by: §1, §2, §2, 1st item, §6.1.
  • [20] K. Luck, D. Kim, L. Lambourne, K. Spirohn, B. E. Begg, W. Bian, R. Brignall, T. Cafarelli, F. J. Campos-Laborie, B. Charloteaux, et al. (2020) A reference map of the human binary protein interactome. Nature 580 (7803), pp. 402–408. Cited by: §1, 3rd item.
  • [21] Y. Luo, X. Zhao, J. Zhou, J. Yang, Y. Zhang, W. Kuang, J. Peng, L. Chen, and J. Zeng (2017) A network integration approach for drug-target interaction prediction and computational drug repositioning from heterogeneous information. Nature communications 8 (1), pp. 1–13. Cited by: §1.
  • [22] L. v. d. Maaten and G. Hinton (2008) Visualizing data using t-sne. Journal of machine learning research 9 (Nov), pp. 2579–2605. Cited by: Fig. 3, §6.2.
  • [23] E. L. Michalets and C. R. Williams (2000) Drug interactions with cisapride. Clinical pharmacokinetics 39 (1), pp. 49–75. Cited by: TABLE V.
  • [24] A. Niculescu-Mizil and R. Caruana (2005)

    Predicting good probabilities with supervised learning

    In Proceedings of the 22nd international conference on Machine learning, pp. 625–632. Cited by: §6.4, §6.4.
  • [25] B. Ozbek and G. Otuk (2009) Post-antibiotic effect of levofloxacin and tobramycin alone or in combination with cefepime against pseudomonas aeruginosa. Chemotherapy 55 (6), pp. 446–450. Cited by: TABLE V.
  • [26] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al. (2019) Pytorch: an imperative style, high-performance deep learning library. In Advances in neural information processing systems, pp. 8026–8037. Cited by: §5.3.
  • [27] B. Perozzi, R. Al-Rfou, and S. Skiena (2014) Deepwalk: online learning of social representations. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 701–710. Cited by: §1, §2, 1st item.
  • [28] E. Perucca (2006) Clinically relevant drug interactions with antiepileptic drugs. British journal of clinical pharmacology 61 (3), pp. 246–255. Cited by: TABLE V.
  • [29] J. Piñero, J. M. Ramírez-Anguita, J. Saüch-Pitarch, F. Ronzano, E. Centeno, F. Sanz, and L. I. Furlong (2020) The disgenet knowledge platform for disease genomics: 2019 update. Nucleic acids research 48 (D1), pp. D845–D855. Cited by: 4th item, §6.6, §6.6, TABLE IV.
  • [30] J. Qiu, J. Tang, H. Ma, Y. Dong, K. Wang, and J. Tang (2018) DeepInf: social influence prediction with deep learning. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 2110–2119. External Links: ISBN 9781450355520 Cited by: §1.
  • [31] R. A. Rossi, N. K. Ahmed, E. Koh, S. Kim, A. Rao, and Y. A. Yadkori (2018) HONE: higher-order network embeddings. arXiv preprint arXiv:1801.09303. Cited by: §3.1.
  • [32] R. A. Rossi, N. K. Ahmed, and E. Koh (2018) Higher-order network representation learning. In Companion Proceedings of the The Web Conference 2018, WWW ’18, Republic and Canton of Geneva, CHE, pp. 3–4. External Links: ISBN 9781450356404, Link, Document Cited by: §3.1.
  • [33] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov (2014) Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research 15 (1), pp. 1929–1958. Cited by: §5.3.
  • [34] T. Truong and J. Haley (2012) Probable warfarin and dapsone interaction. Clinical and Applied Thrombosis/Hemostasis 18 (1), pp. 107–109. Cited by: TABLE V.
  • [35] D. S. Wishart, C. Knox, A. C. Guo, S. Shrivastava, M. Hassanali, P. Stothard, Z. Chang, and J. Woolsey (2006) DrugBank: a comprehensive resource for in silico drug discovery and exploration. Nucleic acids research 34 (suppl_1), pp. D668–D672. Cited by: §6.2, §6.6, §6.6.
  • [36] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and S. Y. Philip (2020) A comprehensive survey on graph neural networks. IEEE Transactions on Neural Networks and Learning Systems. Cited by: §3.2.
  • [37] R. Ying, R. He, K. Chen, P. Eksombatchai, W. L. Hamilton, and J. Leskovec (2018)

    Graph convolutional neural networks for web-scale recommender systems

    In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 974–983. External Links: ISBN 9781450355520 Cited by: §1.
  • [38] Z. Ying, D. Bourgeois, J. You, M. Zitnik, and J. Leskovec (2019) Gnnexplainer: generating explanations for graph neural networks. In Advances in neural information processing systems, pp. 9244–9255. Cited by: §7.
  • [39] X. Yue, Z. Wang, J. Huang, S. Parthasarathy, S. Moosavinasab, Y. Huang, S. M. Lin, W. Zhang, P. Zhang, and H. Sun (2020) Graph embedding on biomedical networks: Methods, applications and evaluations. Bioinformatics 36 (4), pp. 1241–1251. External Links: Document, 1906.05017, ISSN 14602059 Cited by: §1.
  • [40] D. Zhu, P. Cui, Z. Zhang, J. Pei, and W. Zhu (2018) High-order proximity preserved embedding for dynamic networks. IEEE Transactions on Knowledge and Data Engineering 30 (11), pp. 2134–2144. Cited by: §3.1.
  • [41] M. Zitnik, M. Agrawal, and J. Leskovec (2018) Modeling polypharmacy side effects with graph convolutional networks. Bioinformatics 34 (13), pp. i457–i466. Cited by: §1, §1.
  • [42] M. Zitnik, S. Rok Sosic, and J. Leskovec (2018) BioSNAP datasets: stanford biomedical network dataset collection. Cited by: 1st item, 2nd item.