A G RAPH F EATURE AUTO -E NCODER FOR THE P REDICTION OF

U NOBSERVED N ODE F EATURES ON B IOLOGICAL N ETWORKS

A P REPRINT

Ramin Hasibi Tom Michoel

arXiv:2005.03961v2 [[Link]] 23 Dec 2020

Department of Informatics Department of Informatics

University of Bergen University of Bergen
Bergen, 5008 Norway Bergen, 5008 Norway
[Link]@[Link] [Link]@[Link]

December 24, 2020

A BSTRACT
Motivation: Molecular interaction networks summarize complex biological processes as graphs,
whose structure is informative of biological function at multiple scales. Simultaneously, omics
technologies measure the variation or activity of genes, proteins, or metabolites across individuals or
experimental conditions. Integrating the complementary viewpoints of biological networks and omics
data is an important task in bioinformatics, but existing methods treat networks as discrete structures,
which are intrinsically difficult to integrate with continuous node features or activity measures. Graph
neural networks map graph nodes into a low-dimensional vector space representation, and can be
trained to preserve both the local graph structure and the similarity between node features.

Results: We studied the representation of transcriptional, protein-protein and genetic interaction

networks in E. coli and mouse using graph neural networks. We found that such representations
explain a large proportion of variation in gene expression data, and that using gene expression data as
node features improves the reconstruction of the graph from the embedding. We further proposed a
new end-to-end graph feature auto-encoder which is trained on the feature reconstruction task, and
showed that it performs better at predicting unobserved node features than auto-encoders that are
trained on the graph reconstruction task before learning to predict node features. When applied to
the problem of imputing missing data in single-cell RNAseq data, our graph feature auto-encoder
outperformed a state-of-the-art imputation method that does not use protein interaction information,
showing the benefit of integrating biological networks and omics data using graph representation
learning.

Availability: The source code of this project is written in python and is available on [Link]/
RaminHasibi/GraphFeatureAutoencoder

1 Introduction
Biological networks of genetic, transcriptional, protein-protein, or metabolic interactions summarize complex biological
processes as graphs, whose structure or topology is informative of biological function at multiple scales. For instance,
degree distributions reflect the relative importance of genes or proteins in a cell; 3-4 node network motifs have
well-defined information-processing roles; and network clusters or communities contain genes or proteins involved
in similar biological processes [1, 2, 3]. Simultaneously, genomics, transcriptomics, proteomics, and metabolomics
technologies measure the variation or activity of genes, proteins, or metabolites across individuals or experimental
conditions [4, 5]. There is a rich history of integrating the complementary viewpoints of biological networks and omics
data. For instance, “active subnetwork” identification methods treat omics data as features of network nodes in order
to identify well-connected subnetworks that are perturbed under different conditions [6]. Network propagation or
smoothing methods on the other hand use biological networks to extend partial information on some nodes (e.g., disease
 A PREPRINT - D ECEMBER 24, 2020

association labels, partially observed data) to other nodes (e.g., to discover new disease-associated genes or impute
missing data) [7, 8]. However, existing methods treat biological networks as discrete structures, which are intrinsically
difficult to integrate with continuous node features or activity measures.
Recently, with the advent of deep learning, the idea of representation learning on graphs has been introduced. In this
concept, nodes, subgraphs, or the entire graph are mapped into points in a low-dimensional vector space [9]. These
frameworks are known as Graph Neural Networks (GNNs), and use deep auto-encoders to preserve the local structure
of the graph around each node in the embedding, without having to specify in advance what “local” means. However,
not much attention has been paid so far to the representation of the node features in these embeddings [10, 11].
In this paper, we study whether graph representation learning of biological networks using GNNs results in embeddings
that are compatible with or informative for molecular profile data, concentrating for simplicity on gene expression data.
The three main contributions of this study are:

1. We introduce a method to systematically measure the relationship between the structure of a network and
the node feature (gene expression) values. This is done using the Graph Auto-Encoder approach of [12] and
measuring (i) the performance of reconstructing the network from the embedding, with and without node
feautures, and (ii) measuring the variance in feature values explained by the embedding matrix.
2. We investigate how well GNNs preserve node features (expression values) in their graph representation, and
introduce a new Graph Feature Auto-Encoder layer
that is tailored to reconstructing the representation of the node features rather than the graph structure.
3. We show that our new approach to graph representation learning has practical applications in tasks such as
imputation of missing values in single cell RNA-seq data and similar scenarios.

2 Approach

2.1 Problem Formulation

Assume that an undirected, unweighted graph G = (V, E) with N = |V| number of nodes (genes) has an adjacency
matrix A, where Aij = 1 if there is an edge between nodes i and j and zero otherwise, and degree matrix D, a diagonal
matrix with the degrees (number of neighbours) of each node on the diagonal. An N × Q matrix X and N × P matrix Y,
respectively called the feature matrix and the label matrix, denote the training and testing node feature values (molecular
profiles), referred to as features and labels, of the N nodes in Q + P different experiments. We consider two approaches
for mapping the graph into node feature values. In the first approach, the main purpose is to measure whether the graph
structure and feature values encode similar information. In the second approach, we aim to reconstruct the feature
matrix by leveraging the graph structure.

2.2 Relationship between feature values and graph structure

We learn a N × H node embedding matrix Z of the graph adjacency matrix A using a graph auto-encoder neural
network:
Z = GNN(X, A), (1)
as depicted in Fig. 1(A). We use the metrics of Average Precision and area under the ROC curve related to the
reconstruction of A, and the Variance Explained of Y by Z quantify the relationship between the node features and
graph structure.

2.3 Prediction of unobserved node features

2.4 Neural network on graphs

One of the first attempts at learning neural networks over graph structures was the convolution operation on graphs. For
an input signal x ∈ RN , the spectral convolution is defined as

gθ ∗ x = Ugθ UT x, (2)
in which U is the matrix of eigenvectors of the symmetric Laplacian L = D − A = UΛUT . UT x is called the Fourier
transform of signal x and gθ is a matrix function of Λ, the diagonal matrix of eigenvalues of L.

2
 A PREPRINT - D ECEMBER 24, 2020

Graph
Convulotion
Embedding
Layer
Matrix

Graph
Decoder
(A) (B)
Graph
Structure

Feature
Matrix

Regression
Model

Powered by TCPDF ([Link])

(C)

Figure 1: Different graph auto-encoder schemes. (A) depicts the graph auto encoder scheme of [12] tailored to
reconstruction of the adjacency of the graph. In (B) we take the embedding matrix of graph auto-encoder and train an
indirect regression task for prediction of the expression values. (C) illustrates our proposed approach of graph feature
auto-encoder for end-to-end learning of graph node features.

Due to the high cost of calculating the eigenvalues in the case of large matrices, [13] proposed to use a Chebyshev series
expansion truncated after the K th term to approximate the graph convolution operation with a K th -order polynomial:

K K
X 0 X 0
gθ ∗ x ≈ U θk Tk (Λ̃)UT x = θk Tk (Λ̃)x, (3)
k=0 k=0
0
in which Tk (.) and θk are the k th -order Chebyshev polynomials and expansion coefficients, respectively, Λ̃ =
2 T 2
λmax Λ − 1N with λmax the largest eigenvalue of Λ, and L̃ = UΛ̃U = λmax L − 1N .

In Graph Convolution Networks (GCN) [14], further approximations were done by setting K = 1, λmax ≈ 2, and
0 0
θ = θ0 = −θ1 . As a result, (3) was transformed into
1 1
gθ ∗ x ≈ (1n + D− 2 AD− 2 )x. (4)
− 12 − 12
Repeated application of gθ resulting in high powers of D AD can cause numerical instabilities. [14] suggested to
set the diagonal elements of A to 1 (add self-loops) and to recompute D according to the updated adjacency matrix.
Therefore, they used the symmetrically normalized adjacency matrix à in their convolution layer, with
à = D−1/2 AD−1/2 . (5)

Thus, the forward operation in a Graph Convulotion Network (GCN) for Q input signals arranged in a N × Q matrix X
is computed as
Z = σ(Ã ReLU(ÃXW0 )W1 ) (6)
with weight matrices Wi containing the trainable weights for each input feature, and σ a non-linear task specific
function such as softmax for a node classification problem [15].
Additional studies on GNNs have show that a GNN can be viewed as a message-passing approach based on graph
structure, where every node’s message for its neighbours is the aggregation of the neighbourhood information, where
the aggregation is done through a trainable neural network [16]. This framework is also known as a Message Passing
Neural Network (MPNN). The forward pass in such a network consists of a message passing phase and a readout phate.
In the message passing phase, the hidden representation of each node is updated through an update function, which
aggregates the previous node representation and the messages of its neighbours according to :
X
hki = γ k (hk−1
i , M (hk−1
i , hk−1
j , eij )), (7)
j∈N (i)

3
 A PREPRINT - D ECEMBER 24, 2020

in which hki is the hidden state of node i in layer k, with h0i being the node i’s features and eij is the edge attribute
between nodes i and j. Additionally, γ and M are both differentiable functions called the update and message functions,
respectively. Furthermore, N (i) denotes the set of neighbouring nodes of node i.
In the readout phase, the feature vector of the graph is computed using some learnable, differentiable readout function
R according to
Y = R({hi |i ∈ V}). (8)

3 Methods
3.1 Graph representation learning for structural embedding

To learn the representation of the structure of a graoh, we used the Non-probabilistic Graph Auto-Encoder model of
[12]. The embedding matrix Z is calculated by
Z = GCN(X, A), (9)
where GCN(X,A) is a two layer Graph Convolutional Network on the input A and X (the features of the graph).
GCN(X,A) is obtained by setting σ in (6) to the identity function,
GCN(X, A) = ÃReLU(ÃXW0 )W1 . (10)
The weights are trained by measuring how well the embedding reconstructs the graph adjacency matrix, where the
reconstructed adjacancy matrix  is defined as
 = Sigmoid(ZZT ). (11)
The cross-entropy error over all the edges in the matrix is used as a loss function,
N
X
L=− An ln Ân , (12)
n=1

in which An is the adjacency row of the nth node in A and Ã. The training of the neural network is done by gradient
descent and stochasticity added by dropout rate.

3.2 Message passing neural network for end-to-end prediction of features

We used three popular message passing networks for finding the hidden representation of the nodes, as well as
introducing our own, for the task of predicting node feature values based on their network neighbourhood information.
These three methods are inductive GCN, GraphSAGE [17], and the GNN operator from [18] (from here on out refered
to as GraphConv). According to [16], the formula for inductive GCN is
X 1
hki = p p .(Whk−1
j ). (13)
j∈N (i)∪i
deg(i) ∗ deg(j)
The formula for GraphSAGE is
hki = W1 (hk−1
i ) + W2 Meanj∈N (i)∪i (hk−1
j ). (14)
The GraphConv operator is calculated through
X
hki = W1 hk−1
i + W2 .hk−1
j . (15)
j∈N (i)

Our version of MPNN predicts the features of a node based on its own features and those of its neighbours. Therefore,
we first obtain a representation of every node’s features by running them through a linear layer. Then, we aggregate the
neighbouring information by mean pooling of messages. For the update function we concatenate the node features with
its aggregated message representation, and run them through a shared weight network which determines how important
each of these values are in predicting the features of the node. Hence, the formulae for our FeatGraphConv operator are
as follows
gik = W1 ∗ hk−1
i (16)
hki = W2 (gik ||Meanj∈N (i)∪i (gjk )), (17)
in which (.||.) is the concatanation function.

4
 A PREPRINT - D ECEMBER 24, 2020

Table 1: Hyper-parameters of the graph neural network

Hyper-parameter Node Embedding MPNN
Epochs 500 20000
initial learning rate 0.001 0.001
first hidden layer size 64 64
second hidden layer size 32 32

3.3 Prediction of features from features

For comparison of the results obtained through GNNs, we also consider prediction of Y directly from X through simple
machine learning algorithms. These algorithms include:
• Multi Layered Perceptron (MLP): A simple form of neural networks which maps the input features into
output features through multiple layers of neurons (computing units).
• Linear Regression: A linear model for mapping the input to output.
• Random Forest: A set of decision tree models that determine the output value through the aggregation of
output of decision trees that each are trained on a subset of X
• Markov Affinity-based Graph Imputation of Cells (MAGIC): Uses signal-processing principles similar to
those used to clarify blurry and grainy images to recover missing values from already existing ones in a matrix.
[19]

3.4 Experimental setup

The hyper-parameters of all the experiments were determined after some initial experiments and were kept the same for
all the models, to measure the predictive performance of different approaches under the same set of initial circumstances.
These hyper-parameters are listed in Table 1.

3.4.1 Train-test splitting of features and nodes and edges

We divided different properties of the graph into train and test sets depending on the experiment. When experimenting
on the relationship between different graphs and the node features (cf. Section 2.2), for the task of calculating the
average percision and AUC of adjacency reconstruction we follow the work of [12], and split the E of each graph into 3
folds, 2 of which are for training and the remaining is used for testing. For calculating the Variance Explained and
Average Mean Square Error (MSE) of prediction of each dimension of Y, we split the nodes in 3 folds, 2 of which
are used for training and the remaining for testing. As a result, a regression model is trained on the embedding of the
training nodes for each of the features (experiments) and is used to predict or measure the Variance Explained on the
feature for the test nodes (Fig 1.a).
For assessing the feature auto-encoder scheme (cf. Section 2.3), we randomly divide all the features values into 3
separate folds. The aim is to reconstruct the training values and measure the error on the test values (Fig 1.b). After
every run the folds are exchanged in training and testing, such that the error measured is the mean of all the test errors
in 3 different runs. This procedure is also known as K-fold cross validation.

3.4.2 Structure embedding

We used the PytorchGeometry implementation of the graph auto-encoder provided by [20]. For our approach, the
normal auto-encoder provided in the package was used, and the variational auto-encoder was omitted. Four different
sets of input graphs and features to the model were tested:
1. Random Graph: In this approach Z is calculated by
Z = GCN(1N , Arand ), (18)
in which, 1N and Arand represent an identity matrix of size N × N and the adjacency matrix of a random
graph, respectively. For generating random graphs, we used the random graph generator of the Python3
package NetworkX, using the Erdős–Rényi model [21] (Appendix)
2. Expression + Random Graph: In this approach, the identity matrix is replaced with the actual expression
values of genes as input features:
Z = GCN(X, Arand ). (19)

5
 A PREPRINT - D ECEMBER 24, 2020

Train Nodes Train Values

Test Values

Test Nodes

(a) Separate model for each dimension (b) One model for values

Figure 2: The two different approaches of prediction of expression values. In (a), each dimension of Y is predicted
using a separate model trained on X. In (b), a single model is used to predict all the values of the feature matrix.

3. Real Graph: In this approach, the embedding matrix Z is calculated by

Z = GCN(1N , A), (20)
where A is the adjacency matrix of the input graph.
4. Expression + Real Graph: The embedding Z in this case is calculated through (9),
Z = GCN(X, A).
5. Expression: The network in this model, is inferred from the correlation between the expression values of
different genes. In this approach, the correlation directly outputs the probability of the edge between two
nodes.
By choosing the identity matrix as input features in approach 1 and 3, each one of the nodes has a distinct set of features
which do not have any indication to the functionality of the node. This way the model will only pay attention to the
graph structure when producing the embedding matrix.

3.4.3 Node feature representation learning

Since Y has P number of dimensions, we can either train a different model for prediction of each of the dimensions, or
predict Y in an auto-encoder manner. In the latter, the values of Y are considered to be missing values in the feature
matrix and the model is trained on reconstruction of the X values. In the first approach, the evaluation of the methods
is done on the prediction of each of the Y dimensions of the test genes. In the second approach, the Y values are
scattered in the feature matrix, similar to the case of imputation problems. The difference between these two methods is
illustrated in Fig. 2.
For the MPNN modules of inductive GCN, GraphSage, and GraphConv we used the already available implementations
from the package PytorchGeo. Furthermore, we used MessagePassing class of this package to implement our proposed
FeatGraphConv. The hyper-parameters of all the MPNNs are available in Table 1. The architectures of all the models
consist of two message passing layers followed by a fully connected layer as the readout function, with Relu as the
activation function of the message passing layers. As before, both the cases of with node features and without node
features are used in the training of the MPNN. In the case of nodes without features, 1N is used, turning (3) into
Y = GNN(1N , A). (21)

3.4.4 Datasets
We evaluated the performance of our method on data for the organisms Escherichia Coli (E. coli) and Mus Musculus.
• Network datasets: For E. coli, we used transcriptional, protein-protein, and genetic interaction networks. The
transcription network was obtained from RegulonDB [22]. All the positive and negative regulatory effects from
the TF-gene interactions dataset file were included regardless of their degree of evidence (strong or weak) to
construct the adjacency matrix. A PPI and genetic interaction network were obtained from BioGRID [23, 24].
We extracted the interactions from the file ”BIOGRID-ORGANISM-Escherichia_coli_K12_W3110-3.5.180”,
and considered the ”physical” and ”genetic” values of ’Experimental System Type’ column for constructing
the PPI and genetic networks, respectively. For Mus Musculus, we used a protein-protein interaction network
extracted in the same way from the file ”BIOGRID-ORGANISM-Mus_musculus-3.5.182”.

6
 A PREPRINT - D ECEMBER 24, 2020

Table 2: Summary description of benchmark datasets

TF_net PPI Genetic
Input Expressions
Nodes Edges Nodes Edges Nodes Edges
Ecoli 466 1559 3184 1929 11592 3688 147475
Mus Musculus 1468 - - 9951 75587 - -

• Expression level dataset: For E. coli, we used the Many Microbes Microarray Database (M3DB) [25, 26].
All the experiments from the file ”avg_E_coli_v4_Build_6_exps466probes4297” were used to construct the
feature matrix. For Mus Musculus, the single cell RNA-Seq data from [27] were obtained from the Gene
Expression Omnibus.
For each of the networks, the common genes between the network and the expression data were extracted, and an
adjacency matrix and a matrix of features were constructed from the network and expression level datasets, respectively.
A detailed description of each of the networks is available in table 2.

3.5 Computational resources and Source Code

All the experiments were done on Tesla V100 with Python 3. The source code of the experiments is availabe at
[Link]/RaminHasibi/GraphFeatureAutoencoder

4 Results
4.1 Measuring the relationship between biological networks and gene expression values

To explore whether the structure of biological networks relates to the expression values of genes, and whether such a
relationship can be expressed quantitatively using graph neural networks, we learned low-dimensional embeddings
of transcriptional regulatory (TF_net), protein-protein interaction (PPI) and genetic interaction networks in E. coli
and a protein-protein interaction network in mouse, with and without expression data as node features, and trained to
optimize reconstruction of the original graph from the node embedding (see Methods Section 3 for details). Table 3
shows the area under the ROC curve (AUC) and Average Precision results for the graph reconstruction task for various
embeddings. Embeddings learned from the structure of the real graph alone (“Graph” row) performed considerably
better than embeddings learned from random graphs (“Random Graph” row), as expected. The same was true for a
standard Pearson correlation coexpression network inferred from the expression data alone (“Expression” row), showing
that graph embeddings and gene expression data independently predict graph structure. Interestingly, when gene
expression data were used as node features during the graph representation learning (see Section 3.1 and 3.4.2), graph
reconstruction performance further increased (“Expression + Graph” row), but this was not the case when expression
data was combined with random graphs (“Expression + Random Graph” row). In other words, graph embeddings where
the distance between nodes respects both their graph topological distance and their expression similarity results in better
graph reconstruction than embeddings that are based on topological information alone. This shows that expression
profiles are informative of graph structure in a way that is consistent with, but different from, the traditional view where
networks are inferred directly from expression data using expression similarity measures.
Next we computed the variance of the expression data explained by the different embeddings (see Appendix A.1).
Despite being trained on the graph reconstruction task, graph embeddings learned with and without expression data as
node features explained a high percentage of variation in the expression data, but not when random graphs were used
(Table 4).
In summary, graph representation learning results in low-dimensional node embeddings that faithfully reconstruct
the original graph as well as explain a high percentage of variation in the expression data, suggesting that graph
representation learning may aid the prediction of unobserved expression data.

4.2 Predicting unobserved features for a fixed set of nodes

To investigate the use of graph neural networks in predicting unobserved node features, we first considered the problem
where feature data is available for a constant set of training or predictor nodes in all experiments, in addition a set of
training experiments is available with feature data for all graph nodes, and the task is to predict the values of the test
nodes from test data for the training nodes (see Fig. 2(a)). As in Section 2, we call the training data (for the train and
test nodes) X the feature matrix and the testing data Y the label matrix, and hence the task is to predict the values of the
test nodes in the label matrix.

7
 A PREPRINT - D ECEMBER 24, 2020

Table 3: Area under the ROC curve (AUC) and average precision (AP) for reconstructing biological network in E. coli
and mouse from five graph structure embedding approaches (detailed Section 3.4.2)
AUC AP
Input Ecoli Mus Musculus Ecoli Mus Musculus
TF_net PPI Genetics PPI TF_net PPI Genetics PPI
Expression + Graph 0.868±0.017 0.8502±0.002 0.894±0.01 0.803±0.000 0.918±0.007 0.872±0.003 0.904±0.006 0.762±0.000
Graph 0.574±0.020 0.8316±0.018 0.882±0.014 0.6674±0.0000 0.727±0.012 0.860±0.010 0.904±0.006 0.653±0.000
Expression + Random Graph 0.529±0.01 0.4958±0.002 0.485±0.005 0.5108±0.0000 0.525±0.011 0.488±0.001 0.489±0.004 0.506±0.000
Random Graph 0.492±0.006 0.4892±0.006 0.490±0.007 0.5025±0.0000 0.496±0.005 0.486±0.005 0.5002±0.003 0.505±0.000
Expression 0.579±0.02 0.580±0.018 0.557±0.02 0.610±0.000 0.611±0.002 0.624±0.018 0.590±0.02 0.627±0.002

Table 4: The average Variance Explained on Gene Expression from the Regulatory Networks graph embedding

Input Ecoli
TF_net PPI Genetics
Expression + Graph 0.883±0.005 0.805±0.003 0.7824±0.006
Graph 0.716±0.010 0.752±0.008 0.6697±0.010
Expression + Random Graph 0.660±0.010 0.667±0.013 0.6618±0.014
Random Graph 0.647±0.017 0.652±0.018 0.6483±0.010

We considered three categories of prediction methods: (i) standard baseline methods that don’t use graph information
(linear regression, random forest regression, and a multi-layer perceptron, see Section 3.3), (ii) standard regression
methods trained on graph embeddings instead of directly on the training data (linear regression (LR) and random forest
(RF) regression), and (iii) graph feature auto-encoder methods for end-to-end learning of features (GCN, GraphSage,
GraphConv, and FeatGraphConv, see Section 3.2), see also Fig. 1.
For the baseline models, we trained only on the expression values in the training data X and predicted the values of
the test nodes for each dimension of Y separately. For the indirect learning frameworks, first the graph embedding is
obtained as in Section 4.1 using the graph auto-encoder, with and without using the training data X as node features,
and then linear regression and random forest prediction models trained on the embedding in order to predict the values
of the test nodes in Y, as in the baseline models. For each of the end-to-end graph convolution models, we trained the
GNNs using the specified network and the training data X as node features and again predicted the values of the test
nodes in Y in each dimension. Table 5 shows the performance (average mean squared error) of all methods on the E.
coli data. For this experiment the mouse single-cell RNA-seq data was omitted due to sparsity of the expression values.
As shown by the results, the newly proposed graph feature auto-encoder FeatGraphConv is able to predict the unobserved
expression values better than the other graph convolution models, due to the fact that this layer is tailored to the prediction
of features rather than the reconstruction of the graph. As expected, all end-to-end methods perform considerably better
when training data is included as node features. The end-to-end methods, with the exception of GCN, also perform
better than the indirect methods where regression models are trained on graph embeddings.
On the other hand, when the graph structure alone is used, the indirect embedding-based methods get less error. This
could be due to the fact that these models better capture the structure of the graph, since their loss function is defined on
the reconstruction of the adjacency matrix. Hence when the graph structure is the only information provided to the
model, they are able to better capture this information and therefore obtain an embedding that better predicts expression
data (on the basis of the results in Section 4.1), compared to end-to-end models which try to predict the features directly
and are operating blindly when no training features are provided.
We also observe that the lowest MSE overall is in fact obtained by baseline linear regression on the training data
alone. However, since LR performs also better than MLP and RF in this context, we speculate that this result is due to
the relative simplicity of the expression data with highly linear correlations among genes, and would not necessarily
generalize to other data. Moreover, experiments on FeatGraphConv with 20,000 iterations (as opposed to the default
of 500 used for all end-to-end methods in Table 5) showed that this model can decrease the MSE to 0.204 ± 0.12,
0.133 ± 0.089, and 0.107 ± 0.083 for each of the TF_net, PPI, and Genetic networks, respectively. However, due
to the high number of experiments and the need to train a different model for each of the experiments of Y, it is not
computationally efficient to train the more complex GNN models with a higher number of iterations by default for this
prediction task.

8
 A PREPRINT - D ECEMBER 24, 2020

Table 5: The average MSE of predicting test features of the test nodes using different models.
TF_Net PPI Genetics
Method
Features Graph Features Graph Features Graph
GCN 7.791±3.550 15.127±2.280 6.208±0.607 11.106±0.52198 5.988±0.696 4.560±0.351
GraphSAGE 0.332±0.160 8.078±2.592 0.265±0.135 2.844±0.349 0.233±0.127 4.466±1.605
GraphConv 0.318±0.154 13.812±4.534 0.308±0.139 3.094±0.431 0.234±0.116 5.226±1.054
FeatGraphConv 0.285±0.135 7.525±2.941 0.244±0.130 5.207±1.476 0.201±0.112 3.414±0.691
Multilayer perceptron 0.424±0.170 - 0.354±0.153 - 0.332±0.134 -
Linear Regression 0.21526+-0.12626 - 0.147±0.105 - 0.108±0.084 -
Random Forest 0.50969+-0.14335 - 0.194 ±0.103 - 1.882±0.343 -
LR-embedding 1.583±0.200 2.279±0.403 1.091±0.166 1.453±0.264 1.863±0.332 1.653±0.271
RF-embedding 1.945±0.31759 2.150±0.363 1.472±0.267 1.452±0.267 1.883±0.343 1.897±0.351

Table 6: The average MSE of predicting randomly distributed test values using different auto-encoder models.
Ecoli Mus Musclus
Model TF_Net PPI Genetics PPI
GCN 0.441±0.045 0.314±0.019 0.515±0.038 0.083±0.313
GraphSAGE 0.071±0.001 0.056±0.0003 0.059±0.001 0.015±0.003
GraphConv 0.144±0.003 0.425±0.187 0.487+-0.245 0.015±0.001
FeatGraphConv 0.059±0.002 0.049±0.002 0.044±0.001 0.010±0.0003
MLP Auto-encoder 0.062±0.004 0.053±0.004 0.049±0.006 0.013±0.002
MAGIC 3.505±0.006 3.661±0.017 3.215±0.003 0.050±0.0002

4.3 Predicting randomly distributed missing node features

Based on the results in the previous section, we next considered the more challenging prediction task where unobserved
node features are randomly distributed over the nodes and differ between experiments, that is, the task of imputing
missing data (Fig. 2(b)). Since there are no fixed sets of training and test nodes, neither the baseline regression methods
nor the indirect GNN learning frameworks are applicable in this case. In contrast, the end-to-end graph convolution
methods allow to train a single model in an auto-encoder manner for the prediction of all the Y values, which may be
placed in any possible order inside the feature matrix. We used these GNN models for the prediction of expression
values in E. coli and of non-zero values of the single-cell RNAseq data in mouse, and benchmarked them against two
methods that don’t use graph information, namely a normal MLP used in an auto-encoder scheme and MAGIC, a
method designed specifically to impute missing data in single-cell RNA-seq data [19] (see Section 3.3).
As shown in Table 6, our FeatGraphConv convolution layer is able to predict missing features more accuratly compared
to all other methods. It is interesting to note that the GNN models, with the exception of GCN, outperform MAGIC on
the single-cell RNAseq imputation task, although the MLP, which does not use graph information, also performs well in
this case.

5 Discussion

In this paper we studied whether graph neural networks, which learn embeddings of nodes of a graph in a low-
dimensional space, can be used to integrate discrete structures such as biological interaction networks with information
on the activity of genes or proteins in certain experimental conditions. Traditionally, this is achieved by for instance
network propagation methods, but these methods do not extract quantitative information from a graph that could be
used for downstream modelling or prediction tasks. Graph neural networks on the other hand can include node features
(gene or protein expression levels) in the learning process, and thus in theory can learn a representation that better
respects the information contained in both data types. Thus far the integration of node features in graph representation
learning has mainly been pursued for the task of link prediction. Here instead we focused on the task of predicting
unobserved or missing node feature values.
We showed that representations learned from a graph and a set of expression profiles simultaneously result in better
reconstruction of the original graph and higher expression variance explained than using either data type alone, even
when the representations are trained on the graph reconstruction task. We further proposed a new end-to-end graph
feature auto-encoder which is trained on the feature reconstruction task, and showed that it performs better at predicting

9
 A PREPRINT - D ECEMBER 24, 2020

unobserved node features than auto-encoders that are trained on the graph reconstruction task before learning to predict
node features.
Predicting or imputing unobserved node features is a common task in bioinformatics. In this paper we demonstrated
the value of our graph feature auto-encoder on the problem of imputing missing data in single-cell RNAseq data,
where it performs better than a state-of-the-art method that does not include protein interaction data. Other potential
application areas are the prediction of new disease-associated genes from a seed list of known disease genes on the basis
of network proximity [7], or the prediction of non-measured transcripts or proteins from new low-cost, high-throughput
transcriptomics and proteomics technologies that only measure a select panel of genes or proteins [28, 29].
A potential drawback of our method is that it assumes that the interaction graph is known and of high-quality. Future
work could investigate whether it is feasible to learn graph representations that can do link prediction and node feature
prediction simultaneously, or whether network inference followed by graph representations learning for one type of
omics data (e.g. bulk RNAseq data) can aid the prediction of another type of omics data (e.g. single-cell RNAseq).
In summary, graph representation learning using graph neural networks is a powerful approach for integrating and
exploiting the close relation between molecular interaction networks and functional genomics data, not only for network
link prediction, but also for the prediction of unobserved functional data.

References
[1] A-L Barabási and Z N Oltvai. Network biology: understanding the cell’s functional organization. Nat Rev Genet,
5:101–113, 2004.
[2] X. Zhu, M. Gerstein, and M. Snyder. Getting connected: analysis and principles of biological networks. Genes &
Development, 21(9):1010–1024, April 2007.
[3] Uri Alon. An introduction to systems biology: design principles of biological circuits. CRC Press, 2nd edition,
2020.
[4] Marylyn D Ritchie, Emily R Holzinger, Ruowang Li, Sarah A Pendergrass, and Dokyoon Kim. Methods of
integrating data to uncover genotype-phenotype interactions. Nature Reviews Genetics, 16(2):85–97, 2015.
[5] Yehudit Hasin, Marcus Seldin, and Aldons Lusis. Multi-omics approaches to disease. Genome Biology, 18(1):1–15,
2017.
[6] Hung Nguyen, Sangam Shrestha, Duc Tran, Adib Shafi, Sorin Draghici, and Tin Nguyen. A comprehensive survey
of tools and software for active subnetwork identification. Frontiers in genetics, 10:155, 2019.
[7] Lenore Cowen, Trey Ideker, Benjamin J Raphael, and Roded Sharan. Network propagation: a universal amplifier
of genetic associations. Nature Reviews Genetics, 18(9):551, 2017.
[8] Jonathan Ronen and Altuna Akalin. netSmooth: Network-smoothing based imputation for single cell RNA-seq.
F1000Research, 7, 2018.
[9] William L. Hamilton, Rex Ying, and Jure Leskovec. Representation learning on graphs: Methods and applications.
CoRR, abs/1709.05584, 2017.
[10] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu. A comprehensive
survey on graph neural networks. CoRR, abs/1901.00596, 2019.
[11] Ines Chami, Sami Abu-El-Haija, Bryan Perozzi, Christopher Ré, and Kevin Murphy. Machine learning on graphs:
A model and comprehensive taxonomy. arXiv preprint arXiv:2005.03675, 2020.
[12] Thomas N. Kipf and Max Welling. Variational graph auto-encoders, 2016.
[13] David K. Hammond, Pierre Vandergheynst, and Rémi Gribonval. Wavelets on graphs via spectral graph theory.
Applied and Computational Harmonic Analysis, 30(2):129 – 150, 2011.
[14] Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. CoRR,
abs/1609.02907, 2016.
[15] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu. A comprehensive
survey on graph neural networks. CoRR, abs/1901.00596, 2019.
[16] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing
for quantum chemistry. CoRR, abs/1704.01212, 2017.
[17] William L. Hamilton, Rex Ying, and Jure Leskovec. Inductive representation learning on large graphs. CoRR,
abs/1706.02216, 2017.

10
 A PREPRINT - D ECEMBER 24, 2020

[18] Christopher Morris, Martin Ritzert, Matthias Fey, William L. Hamilton, Jan Eric Lenssen, Gaurav Rattan, and
Martin Grohe. Weisfeiler and leman go neural: Higher-order graph neural networks. CoRR, abs/1810.02244,
2018.
[19] David van Dijk, Roshan Sharma, Juozas Nainys, Kristina Yim, Pooja Kathail, Ambrose J. Carr, Cassandra
Burdziak, Kevin R. Moon, Christine L. Chaffer, Diwakar Pattabiraman, Brian Bierie, Linas Mazutis, Guy Wolf,
Smita Krishnaswamy, and Dana Pe’er. Recovering gene interactions from single-cell data using data diffusion.
Cell, 174(3):716–729.e27, July 2018.
[20] Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR Workshop
on Representation Learning on Graphs and Manifolds, 2019.
[21] P. Erdös and A. Rényi. On random graphs i. Publicationes Mathematicae Debrecen, 6:290, 1959.
[22] RegulonDB. Regulondb network interactions - gene interaction, 2019. data retrieved from RegulonDB Download-
able Experimental Datasets, [Link]
[23] BioGRID. Biogrid ppi interaction network ecoli, 2019. data retrieved from BioGRID download page available at,
[Link]
[24] Rose Oughtred, Chris Stark, Bobby-Joe Breitkreutz, Jennifer Rust, Lorrie Boucher, Christie Chang, Nadine Kolas,
Lara O’Donnell, Genie Leung, Rochelle McAdam, Frederick Zhang, Sonam Dolma, Andrew Willems, Jasmin
Coulombe-Huntington, Andrew Chatr-aryamontri, Kara Dolinski, and Mike Tyers. The BioGRID interaction
database: 2019 update. Nucleic Acids Research, 47(D1):D529–D541, 11 2018.
[25] Jeremiah J. Faith, Michael E. Driscoll, Vincent A. Fusaro, Elissa J. Cosgrove, Boris Hayete, Frank S. Juhn,
Stephen J. Schneider, and Timothy S. Gardner. Many microbe microarrays database. 2007. data retrieved from
M3DB download page available at,[Link]
[26] Jeremiah J. Faith, Michael E. Driscoll, Vincent A. Fusaro, Elissa J. Cosgrove, Boris Hayete, Frank S. Juhn,
Stephen J. Schneider, and Timothy S. Gardner. Many microbe microarrays database: uniformly normalized
affymetrix compendia with structured experimental metadata. Nucleic Acids Research, 36(Database-Issue):866–
870, 2008.
[27] Shangli Cheng, Yu Pei, Liqun He, Guangdun Peng, Björn Reinius, Patrick P.L. Tam, Naihe Jing, and Qiaolin
Deng. Single-cell RNA-seq reveals cellular heterogeneity of pluripotency transition and x chromosome dynamics
during early mouse development. Cell Reports, 26(10):2593–2607.e3, March 2019.
[28] Aravind Subramanian, Rajiv Narayan, Steven M Corsello, David D Peck, Ted E Natoli, Xiaodong Lu, Joshua
Gould, John F Davis, Andrew A Tubelli, Jacob K Asiedu, et al. A next generation connectivity map: L1000
platform and the first 1,000,000 profiles. Cell, 171(6):1437–1452, 2017.
[29] Karsten Suhre, Mark I McCarthy, and Jochen M Schwenk. Genetics meets proteomics: perspectives for large
population-based studies. Nature Reviews Genetics, pages 1–19, 2020.

11
 A PREPRINT - D ECEMBER 24, 2020

Appendix
A.1 Variance explained on features

For an embedding matrix Z ∈ RN ×H and feature matrix X ∈ RN ×M , we calculated the amount of variance in X
explained by Z as

tr PZ XT X
VZ = (A1)
tr(XT X)
where PZ is the projection matrix onto the subspace of RN spanned by the columns of Z,
PZ = Z(ZT Z)−1 ZT .
Note that if we write the eigendecomposition of XT X as XT X = VT ∆V, then the columns of V corresponding to
th
the nonzero eigenvalues in ∆ are the principal components of X. If Z consistP of the i principal component, then
eq. (A1) reduces to the familiar variance explained by this component, ∆i /( j ∆j ). If Z consists of a single vector
z ∈ RN with unit length, kzk = 1, then eq. (A1) reduces to
N
∆
P i (z T ui )2 ,
X
Vz =
i=1 j ∆j

a weighted sum of the variances explained by each principal component, with weights determined by the extent of
overlap between z and each principal component. Eq. (A1) generalizes this to summing the variances explained by
multiple vectors simultaneously that need not be mutually orthogonal.

A.2 Random Graph Generation using Erdő s–Rényi model

According to Erdő s–Rényi, a random graph G(n, p) has n2 p edges placed at random. The degree distribution of each


node is calculated through:  
n−1 k
P (deg(v) = k) = p (1 − p)n−1−k , (A2)
k
which is the binomial distribution. We adjusted the value of p so that the random graphs would approximately have the
same number of edges as the real networks.

12