# Biological Random Walks: multi-omics integration for disease gene
prioritization

Michele Gentili  Department of Computer, Control, and Management Engineering
Antonio Ruberti, Sapienza University of Rome, Rome, Italy  Leonardo Martini
Department of Computer, Control, and Management Engineering Antonio Ruberti,
Sapienza University of Rome, Rome, Italy  Marialuisa Sponziello  Translational
and Precision Medicine Department Sapienza University of Rome, Rome, Italy
Luca Becchetti  Department of Computer, Control, and Management Engineering
Antonio Ruberti, Sapienza University of Rome, Rome, Italy

###### Abstract

Motivation: Over the past decade, network-based approaches have proven useful
in identifying disease modules within the human interactome, often providing
insights into key mechanisms and guiding the quest for therapeutic targets.
This is all the more important, since experimental investigation of potential
gene candidates is an expensive task, thus not always a feasible option. On
the other hand, many sources of biological information exist beyond the
interactome and an important research direction is the design of effective
techniques for their integration.
Results: In this work, we introduce the Biological Random Walks (BRW) approach
for disease gene prioritization in the human interactome. The proposed
framework leverages multiple biological sources within an integrated
framework. We perform an extensive, comparative study of BRW’s performance
against well-established baselines.
Availability and implementation: All code is publicly available and can be
downloaded at https://github.com/LeoM93/BiologicalRandomWalks. We used
publicly available datasets, details on their retrieval and preprocessing are
provided in the supplementary material.
Supplementary material: Supplementary material available.

##  1 Introduction

In recent years, through the advent of big data, genomics, and quantitative in
silico methodologies, medicine is witnessing tremendous advancements towards
the understanding of the human pathophysiology. [SSA+20]. Gene-disease
associations have been identified by genome-wide association studies (GWAS)
[HS09] and more recently by whole exome or whole genome sequencing studies
[KSL+13]. While many of the mechanisms underlying these associations remain
largely unclear, a growing body of research highlights associations between
groups of interacting proteins and diseases within the so-called “human
interactome”, representing the cellular network of all physical molecular
interactions [BGL11]. A key feature is that disease proteins do not appear to
be uniformly scattered across the interactome [MSK+15], but they are prone to
participation in common biological activities such as, for example, genome
maintenance, cell differentiation or growth signaling, which are the most
relevant pathways in carcinogenesis [ODC+18]. For these reasons, while
traditional single protein (i.e., magic bullet) approaches have limited
effectiveness in addressing complex diseases, network-based ones can prove
useful in identifying disease modules within the interactome, hopefully
providing insights into key mechanisms and guiding the quest for therapeutic
targets. Moreover, experimental investigation of potential gene candidates is
an expensive task, thus not always a feasible option.

The human interactome refers to all protein-protein interactions within a
cell, including regulatory interaction of transcription factors, metabolic
enzyme-coupled interactions, protein complexes, and kinase/substrate
interactions. This network is largely incomplete. Currently, more than 140,000
interactions involving over 13,000 proteins are known (e.g., see [KSSF17,
GNZ+14]). The interactome-based approach to network medicine [BGL11] proved
effective for a number of diseases, e.g., by identifying putative biomarkers
and subtypes, thus allowing a principled approach to drug targeting [BGL11,
ODC+18]. The need for new disease genes (or disease proteins) as putative
candidates for diagnosis, prognosis or treatment, motivated the development of
a number of algorithms for disease genes and modules prediction [GMB15].

Two related classes of methodologies have emerged over the last decade as the
most promising: module-based [GMB15, BGL11] and network propagation [CIRS17,
KBHR08] algorithms. Module-based algorithms find topological, functional or
disease modules in the interactome network, on the hypothesis that these
represent cellular components likely involved in the same disease. Network
propagation (or diffusion-based) algorithms leverage the information flow
through nearby proteins in the network from initial (known) disease genes as
their main ingredient.

While important and effective in many cases, the interactome is only one of
many sources of biological information. An important research direction is the
quest for effective techniques that allow the seamless integration of rich and
heterogeneous biological sources into methods that were originally designed to
leverage the topological features of the interactome [DHH+18, SL20, DBTVOM07].

##### Our contribution

In this work, we introduce the Biological Random Walks (BRW) framework for
disease gene prioritization. The proposed framework leverages the integration
of multiple biological sources within a propagation-based approach. We compare
BRW’s performance against well-established baselines, such as: RWR [NK10],
DIAMOnD [GMB15], DADA [EBEK11] and RWR-M [VTN+19]. In particular, we
investigate BRW’s performance along different axes: i) an in-depth,
comparative analysis on four cancer phenotypes (i.e., breast cancer, lung
adenocarcinoma, papillary thyroid cancer, colorectal adenocarcinoma); ii) a
broad comparative analysis of BRW’s prioritization performance over a wide
spectrum of Mendelian diseases with different characteristics and prior
information available; iii) an external validation using FDA approved drugs
for Breast Cancer treatment, along with an evaluation of the algorithm
results’ stability across multiple population studies. 111Some of the ideas
presented in this submission, albeit in preliminary form and accompanied by a
minimal, internal validation, appeared in [GMP+19]. For the sake of
completeness, Supplementary Section C.4 presents a comparison of our approach
with [GMP+19].

##  2 Materials and Methods

Figure 1: Network aggregation and network flow of Biological Random Walks. a)
shows the preprocessing step that combines PPI and CO-Expression topology to
derive a combined network. The transition matrix of PPI and CO-Expression
network are combined using a convex combination to compute an aggregated
transition matrix. Note that what is termed “aggregated network” in tab a) is
actually a _weighted_ network that corresponds to the aggregated transition
matrix. b) Shows the flow of the random walker when the personalization vector
is biased with disease specific information. In this way the flow can
propagate from seed nodes (i.e., known disease genes), but also from node that
are not part of the seed set but biologically similar to it (i.e., biological
teleporting probability (BTP)).

###  2.1 Biological Random Walks

Biological Random Walks (BRW) build on the hypothesis that integrating
different biological information sources may better reflect the complexity of
protein interactions in a cell’s process. In light of this insight, our
approach integrates information on pairwise protein interaction of the
Protein-Protein Interaction network (PPI) [BGL11] with other biological data
in a unified framework. Our approach is to some extent agnostic to the
particular biological data source, as long as it affords a principled notion
of similarity between proteins. In the remainder, we use bold lowercase to
denote vectors and capital, non-bold letters to denote matrices. Given a
vector $\mathbf{x}$, $x_{i}$ denotes its $i$-th entry.

##### Notation.

We represent the PPI as an undirected graph $G=(V,E)$, with genes as vertices.
Any edge $(i,j)$ represents a known protein-protein interaction recorded in
the PPI. We assume $|V|=n$ in the following. For a given gene $i$, $N_{h}(i)$
denotes the subset of $G$’s vertices whose shortest path distance from $i$ is
exactly $h$.

####  2.1.1 Random Walks with Restart

Random Walk with Restart (RWR) [KBHR08] is a diffusion-based method, whose
purpose is identifying disease modules that are topologically “close” to known
disease genes in the interactome. It was shown to outperform other
prioritization algorithms in many cases [NK10]. In a nutshell, this algorithm
can be seen as performing multiple random walks over the PPI network, each
starting from a seed node associated to a known disease gene, iteratively
moving from one node to a random neighbour, thus simulating the diffusion of
the disease phenotype across the interactome. More formally, the random walk
with restart is defined as:

$\mathbf{p}^{(t+1)}=(1-r)W\mathbf{p}^{(t)}+r\mathbf{q}.$ (1)

Here, $W$ is the column-normalized adjacency matrix of the graph
andm$\mathbf{p}^{(t)}$ is a vector, whose $i$-th entry $p_{i}^{(t)}$ is the
probability of the random walk being at node $i$ at the end of the $t$-th
step. $r\in(0,1)$ is the restart probability, i.e., the probability that the
random walk is restarted from one of the (disease-associated) seed nodes in
the next step. Upon a restart, the probability of restarting the random walk
from some seed node $j$ is $q_{j}$. Vector $\mathbf{q}$ is normally called a
_personalization vector_ in the Data Science literature. This random walk
corresponds to an ergodic Markov chain [LP17] that admits a stationary
distribution (i.e., a fixed point) $\mathbf{p}$. Nodes of the PPI are simply
ranked by considering the corresponding entries of $\mathbf{p}$ in descending
order of magnitude.

####  2.1.2 Biological Information-Aware Random Walks

For the sake of exposition, in the remainder we refer to the biological
information associated to a gene (e.g., the set of its annotations) as the set
of its features. These can include (more precisely, be derived from)
annotations from the Gene Ontology database [Con19] (GO in the remainder) or
gene expression levels. We remark that, in principle, any reliable information
source on gene biology can be integrated. BRW ranks genes according to the
main steps outlined below.

Unlike [KBHR08] and similar approaches, our method consists of two main steps:
i) extracting statistically significant features from biological data, using
them to compute a personalization vector and a transition matrix used by the
BRW algorithm; ii) using the stationary distribution of the corresponding
random walk to rank genes. Our approach to computing aggregated
personalization vectors and transition matrices is outlined below, with
further details given in Sections B.3 and B.4 of the supplementary material.

##### Computing a Personalization Vector.

Both gene annotation data and gene expression levels allow derivation of
personalization vectors that reflect some notion of similarity between a gene
and a disease, with the latter represented by the set of its seed genes. In
the first case, this similarity is defined in terms of knowledge about genes’
involvement in various biological functions and/or diseases. In the second
case, similarity is defined in terms of co-expression levels of different
genes in subject cases as opposed to expression levels in a control group,
using data about a population of patients involved in a clinical trial.

For annotation data, assume we have $\ell$ sources of biological information
(e.g., GO, KEGG pathways etc.).Let $S$ denote the seed set. For every
$j=1,\ldots,\ell$, we use $\mathcal{F}^{j}$ to denote the subset of
annotations from the $j$-th source that are associated with at least one gene
in $S$. Then, for every $j=1,\ldots,\ell$, we select a subset of annotations
$\widehat{\mathcal{F}^{j}}$, filtering out annotations that are not
statistically significant, as shown in Supplementary Figure 2 (i.e.,
$p-value>10^{-5}$, using Fisher Exact Test and FDR correction), so that
$\mathcal{F}=\cup_{j=1}^{\ell}\widehat{\mathcal{F}^{j}}$ denotes the set of
all statistically significant annotations that are associated with genes in
$S$.Likewise, for every gene $i$ (not necessarily belonging to $S$), we denote
by $\mathcal{A}(i)$ the set of its annotations, possibly extracted from
multiple biological information sources. We assign each gene $i$ a weight
$\theta_{i}$, which reflects the extent to which $i$ shares annotations that
are statistically significant for genes that belong to the seed set of the
disease under consideration. While other choices are possible, the definition
we adopted reflects the extent of the _inclusion_ of $\mathcal{A}_{i}$ in each
of the $\ell$ sources:

$\theta_{i}=\sum_{j=1}^{\ell}\frac{|\mathcal{A}(i)\cap\widehat{\mathcal{F}^{j}}|}{|\widehat{\mathcal{F}^{j}}|}.$
(2)

At this point, the components of a personalization vector $\mathbf{q}$ can be
computed as follows:

$q_{i}=\frac{\theta_{i}}{\sum_{i=0}^{N}\theta_{i}},$ (3)

with $N$ the number of genes we consider. Note that $q_{i}$ denotes the
probability that, upon teleportation, the random walkers jumps to gene $i$.
Further details are given in Section B.3.1 of the Supplementary material.

We next discuss how to compute personalization vectors from gene expression
data. An important goal is the identification of Differentially Expressed (DE)
genes, whose expression levels systematically differ between case (Breast
Cancer Cells) and a control group (Normal Breast Cells). We follow the
approach proposed in [MGS+17], in which subjects of the control group are
assigned personalized perturbation profiles (PEEPs), from which gene
expression-aware personalization vectors can be derived.Succinctly put, for
each gene $i$ and for each subject $j$, the expression level $l_{i}^{j}$ of
gene $i$ in subject $j$ is compared with the distribution of the expression
level of gene $i$ within the control group by taking the corresponding z-score
$z_{ij}$. This approach allows association of a “bar code” to each subject.
Following [MGS+17], we set $|z_{ij}|>2.5$ as the threshold to declare gene $i$
differentially expressed in subject $j$. Following the general intuition
stated in the introduction that disease genes generally are not scattered
across the interactome, we also bias our choice towards differentially
expressed genes that are closer to disease genes in the PPI. Eventually, we
obtain a personalization vector $\mathbf{q}$ that reflects both genes’
differential expressions and vicinity to disease genes. For full details on
computing PEEPs and deriving gene expression-aware personalization vectors, we
refer the reader to Supplementary Section B.3.1.

##### Computing a Transition Matrix.

Similar approaches can be used to derive a transition matrix for therandom
walk with restart. Both approaches rely on the PPI, differing on the way PPI’s
edges are assigned weights that reflect genes’ similarity and determine the
probabilties of edge traversals. We leverage categorical, biological
information (e.g., gene annotations)by defining a weighted transition matrix
$W$, in which each entry $W_{ij}$ depends on the extent to which nodes/genes
$i$ and $j$ of the PPI share common annotations (i.e., they are involved in
common biological processes) that are also significant for the disease. In
more detail, considered genes $i$ and $j$, we define the following _Disease
Specific Interaction Function_ (DSI function in the remainder):

$DSI(i,j)=|\mathcal{A}(i)\cap\mathcal{A}(j)\cap\mathcal{F}|,$ (4)

where we remind that $\mathcal{A}(k)$ denotes the set of gene $k$’s
annotations, while $\mathcal{F}$ denotes the overall set of annotations that
are statistically significant for disease genes.Intuitively, the higher
$DSI(i,j)$, the more $i$ and $j$ share annotations that are also statistically
significant for the disease under consideration.In the end, edges in the PPI
will be assigned weights depending on the DIS function as follows:

$W_{ij}=\begin{cases}c+DSI(i,j)&\text{if }(i,j)\in E,\\\
0&otherwise.\end{cases}$ (5)

Here, $c$ is a positive constant that accounts for usual sparsity of the
available datasets, so that no biological information may be available for the
end-points of a link in the PPI. In this case, the link receives a minimum
weight $c$.As with personalization vectors, gene expression information about
a population of patients can also be used to define a tissue-specific,
population-dependent transition matrix. Specifically, gene expression
information is used to assign weights to edges of the underlying PPI network,
this time reflecting similarities between genes in terms of co-expression with
respect to the subject population.Consider a CO-Expression network in which
each pair of genes $(i,j)$ is assigned a score equal to the Pearson’s
correlation coefficient $pc_{i,j}$ between the expression levels of $i$ and
$j$ within the population. We can define a transition matrix $W$ by assigning
each edge $(i,j)$ of the PPI network a probability as follows:

$W_{ij}=\frac{|pc_{ij}|}{\sum_{k\in N(i)}|pc_{ik}|}$ (6)

where N(i) is the set of i’s neighbors in the PPI network.Note that i) for
every node/gene $i$ we consistently have a probability distribution over its
incident edges and ii) the importance of edges reflects the absolute value of
the correlation between the expression levels of two genes within the
population of interest.

##### Integrating Biological Information and Gene Expression.

We discussed above two orthogonal approaches to the design of personalization
vectors. The first one leverages similarities between the biological processes
associated to known disease proteins and those to be prioritized. Hence,
teleporting probabilities depend on the seed set through association with
common biological processes.In the second case, teleporting probabilities
depend on information that is tissue-specific (the level of expression in a
population of subjects affected by a certain disease) and partly on the seed
set, but this time through the PPI’s topology. Hence, these two approaches
largely rely on complementary sources of information. In order to integrate
these complementary sources into a unique personalization vector that
leverages both, we follow a simple, yet mathematically principled approach,
whereby we take a _convex combination_ of the corresponding personalization
vectors. Namely, assume we have computed two personalization vectors
$\mathbf{q}_{1}$ and $\mathbf{q}_{2}$, the former using biological information
only, the latter using gene expression data and the PPI. We obtain a
personalization vector as follows:

$\mathbf{q}=\alpha\mathbf{q}_{1}+(1-\alpha)\mathbf{q}_{2},$ (7)

where $\alpha\in[0,1]$. Parameter $\alpha$ allows to weigh in the relative
importance of the different information sources we are using. Intuitively,
this type of aggregation amounts to considering a gene a potential candidate
if it is statistically significant in terms of its involvement in biological
processes, of its gene expression levels in the subject group, or both.333Note
that this approach seamlessly extends to an arbitrary number of information
sources. The choice of $\alpha$ (and other parameters of the model), its
impact on performance and dependence of the optimal choice on the scenario at
hand are discussed in detail in Sections 3, 4 and in Section C.2 of the
supplementary material. Remark. Despite its simplicity, this is a
mathematically principled choice. In particular, it is well-known [JW03] and
easy to show that the stationary distribution corresponding to the convex
combination of two personalization vectors $\mathbf{q}_{1}$ and
$\mathbf{q}_{2}$ is itself the linear combination of the stationary
distributions corresponding to $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$
respectively. Briefly put, the parameter $\alpha$ allows us to tune the
relative importance of the information providedby gene annotations and gene
expression levels respectively. Moreover, this approach extends seamlessly to
an arbitrary number of heterogeneous sources of information(and corresponding
personalization vectors). Other aggregation methods were also considered, yet
they provided worse or at most comparable results. Some are presented in
Section B.3.2 of the Supplementary material for the sake of completeness. In
the previous paragraphs, we have seen how we can derive random walk transition
matrices using only information about biological processes (i.e., annotations)
or gene expression data from a population of subjects. As with personalization
vectors, multiple transition matrices derived from complementary biological
sources can be integrated into a single _aggregate transition matrix_. For
example, assume we computed transition matrices $W_{1}$ and $W_{2}$ using
biological and gene expression information respectively. Any convex
combination of $W_{1}$ and $W_{2}$ is a feasible transition matrix. Namely,
for some $\beta\in[0,1]$, we consider the transition matrix

$W=\beta W_{1}+(1-\beta)W_{2}.$

It is easy to see that i) $W$ is still a transition matrix, i.e., the entries
of each row sum to one and that ii) the approach seamlessly extends to any
number of complementary sources. In the case of biological annotations and
gene expression data, $\beta=1$ corresponds to only considering biological
annotations, whereas $\beta=0$ corresponds to only leveraging gene expression
data. So, $\beta$ is a parameter, whose tuning allows us to weigh the
importance of one source of information with respect to the other.

##  3 Results

This section investigates BRW’s performance in prioritizing gene candidates
for genetic diseases.

###  3.1 Experimental Setup

We used a number of biological data sources. Some of them were used as inputs,
to define key parameters of our algorithms, while others were used to
biologically validate the results of the algorithms we considered. They are
briefly described here and more extensively in Section A of the supplementary
material.

##### Data sources.

The experiments discussed in Section 3 were conducted on the HIPPIE Protein-
Protein Interaction network (PPI) [ALANS16] and on the same Protein-Protein
Interaction network (PPI) as in [GMB15] for the sake of comparison.

We used three different sources of gene biological information: Gene Ontology
Consortium444http://geneontology.org/. where, for each gene, we downloaded its
biological processes, KEGG [Kan19, KSF+18, KG00] and Reactome [JMV+19], which
we used to download pathways’ annotations. Gene expression datasets were
downloaded from The Cancer Genome Atlas
database555https://portal.gdc.cancer.gov/.

We validated the methods considered in this study both via an internal
validation on known disease genes and through an external validation using
drug target associations. In more detail, we used disease-gene associations as
in [GMB15] that describe a corpus of 70 Mendelian diseases. From [PRASP+20],
we further derived known disease-gene associations for the four different
cancer types that we investigate in Section 3 (i.e, breast cancer, lung
adenocarcinoma, papillary thyroid cancer, colorectal adenocarcinoma). Finally,
we used Drug-Gene Target associations from
DrugBank666https://go.drugbank.com/. [WDG+17]. We selected only drugs approved
for Breast Cancer treatment from Food and Drug Administration
777https://www.cancer.gov/about-cancer/treatment/drugs/breast. (Supplementary
Table 2).

All data sources mentioned above are more extensively discussed in Section A
of the supplementary material.

Figure 2: Algorithm comparison: a) and b) compare respectively the Recall@K
and nDCG of some of well known network based approaches with BRW on four types
of cancer using a 100-fold Monte Carlo random sampling validation choosing
uniformly at random the 70% of known disease genes and considering the
remaining ones as the test set (30%). c) shows the average shortest path
length of top 200 candidates predicted by the analyzed algorithms. d) shows
the comparison between BRW with the state of the art on a corpus of 70
Mendelian diseases downloaded from [GMB15]. e) compares multi-omics frameworks
(BRW and RWR-M) on average Recall@k on cancer phenotypes when the PPI network
is randomized. The average Recall@K is computed as described in the previous
experiment. f) illustrate how multi-omics integration affect the bias induced
by curated ontologies such as GO, Kegg, and Reactome and the PPI network. It
compares the distribution of p-values (negative log scale) computed by the
testsuit using KEGG pathways and gseapy on gene sets predicted on randomized
PPI networks with those predicted on the original PPI. As hypothesized, BRW,
having in input statistical significant pathways from KEGG, Reactome, and GO,
returns biologically meaningful candidate sets. However, their significance is
not comparable with the set predicted using the original PPI.

##### Baselines.

To provide a robust performance assessment, we compared BRW with a number of
well-known, state-of-art baselines for disease gene prioritization, namely,
RWR [KBHR08], DIAMOnD [GMB15], DADA [EBEK11] and RWR-M [VTN+19]. For the sake
of completeness, supplementary Section C.4 also compares our approach against
an embryonic (and underperforming) version of our framework that was presented
in [GMP+19].

##### Performance indices.

We used the widely adopted indices Recall@K and nDCG in some of the
experiments described in the remainder. To keep presentation self-contained,
they are briefly described in Section C.1 of the supplementary material.

###  3.2 Multi-Omics Integration Improves Algorithm performance

In a first round of experiments, we performed two (internal) validation steps:
i) we first compared the algorithms with respect to four cancer phenotypes
(i.e, breast cancer, lung adenocarcinoma, papillary thyroid cancer, colorectal
adenocarcinoma), for which both biological annotations and gene expression
data are available; ii) as a further step, we performed a broader, yet less
specific validation on a corpus of 70 manually curated mendelian diseases
[GMB15], for which only biological annotations were used. For each disease, we
performed a mean 100-fold validation, by sampling 70% of known disease genes
uniformly at random and using the rest to test the algorithms. For each tested
disorder, we computed Recall@K and the nDCG of each algorithm.

##### Four cancer phenotypes.

We performed a grid search to identify the best combination of hyperparameters
$(\alpha,\beta,r)$. The benchmark, discussed in Section C.2 of the
supplementary material and illustrated in Supplementary Fig. 6 and
Supplementary Fig. 7, highlights interesting, mixed results. For the task of
disease gene prioritization888We remark that the settings discussed here are
not optimal for other tasks in general, e.g., drug target discovery (see
Section 3.4). and to the purpose of computing personalization vectors (i.e.,
teleporting probabilities), the signal contained in statistically significant
annotations derived from seed genes is definitely stronger than the signal
carried by gene expression, so that the best choice for disease gene
prioritization on cancer phenotypes is $\alpha=1.0$, thus completely removing
information about differentially expressed genes. However, gene expression
information is crucial in determining the transition matrix of the random
walk, with values of $\beta\in[0.25,0.5]$ achieving best predictive
performance, indicating that both information sources provide crucial
information to BRW for disease gene ranking.

Results for the four cancer phenotypes show that algorithms that only leverage
the PPI tend to perform worse in terms of both Recall@k and nDCG, as shown in
figure 2 a) and 2 b). Conversely, multi-Omics methods that, like BRW and
RWR-M, rely on multiple biological information sources typically perform
better in terms of the aforementioned indices. Improvement of these methods
over single-source baselines at least in part stems from the well-known fact
that disease-associated genes tend to be involved in similar pathways and
biological processes [BGL11], see also Section 3.3. Interestingly, BRW’s
ability to bias the random walk towards related genes that do not necessary
belong to the seed set seems to play an important role in improving
prioritization of the test set, at least in terms of Recall@k. At the same
time, RWR-M achieves similar performance (slightly better or worse, depending
on the dataset) if one considers nDCG as a global measure of rank, as shown in
figure 2 b). As previously remarked, biased teleporting makes BRW explore
areas of the PPI network that could be relatively far from the seed set, a
fact that is reflected in its candidate genes in the top 200 positions having
higher average shortest path distance than other RWR-based methods that only
teleport to genes of the seed set, as shown in Figure 2 c).

##### Mendelian diseases.

As a further internal validation, we provided a less specific yet broader,
comparative assessment of BRW, by performing a Monte Carlo cross-validation
[DGB07] on “gold standard” gene sets. These sets contain known genes
associated with 70 diseases, which were previously selected in [GMB15] from
OMIM and PheGenI databases. As discussed more in detail at the end of Section
A.1, in this experiment (and only in this one) we used the PPI considered in
[GMB15], for the sake of consistency.999All other experiments use the more
recent HIPPIE-v2.2 PPI network of [ALANS16]. In this second round of
experiments, we set $\alpha=1.0$ and $\beta=1.0$, since gene expression is not
used. Figure 2 d) compares the performances of RwR, DIAMOnD, DADA, RWR-M and
BRW in terms of Recall@k. Our heuristic ranks more known disease proteins than
DIAMOnD in the top 200 positions for 95 percent of disorders analyzed and
respectively for 72, 73, and 63 percent of the disorders for RwR, DaDA, and
RWR-M.

###  3.3 Randomization and Bias

Figure 3: Algorithm Comparison on Breast Cancer Disease: a) percentage of
Breast Cancer Drug targets found by each framework in the top K positions. b)
drugs that are enriched (corrected p-value ¡ $10^{-5}$) in at least one of the
predicted candidate gene sets. c) number of candidate genes, predicted by each
algorithm, that are frequently mutated in at least j of the remaining studies,
with $j\in\\{1,2,3,4\\}$.

Biological Random Walks and RWR-M leverage multiple data sources. As a result,
we expect their results to be less affected by random noise in the PPI with
respect to other heuristics. To test this hypothesis, we performed a first
experiment, by performing an internal validation as done in Section 3.2, but
this time running the random walk-based heuristics using degree-preserving
randomized version of the PPI. To this purpose, we implemented the degree-
preserving randomization algorithm of [MKI+03], also described in Section B.5
of the supplementary material. As hypothesized (Fig. 3 e)), by leveraging
multiple biological sources, BRW and RWR-M are less affected by randomization
of the PPI. This effect is stronger in BRW, whose teleporting probabilities
also depend on phenotype information (statistically significant ontologies,
pathways, and differentially expressed genes). This is further shown in the
more detailed Supplementary Fig. 8 a) of the supplementary material,
highlighting a positive correlation between the number of test genes ranked in
the first 200 positions (Recall@K) and the value of BRW’s restart probability
$r$. As remarked above, this effect is also present in RWR-M, in which
heterogeneous biological sources are summarized in different, layered
networks. In particular, as shown in Supplementary Fig. 8 a) of the
supplementary material, RWR-M’s performance degrades significantly if one also
randomizes the Pathways network used by the algorithm. Altogether, these
results suggest that phenotypical information associated to the nodes plays an
important role in prioritization for these heuristics, owing to the fact that
known disease proteins tend to be involved in the same pathways, a fact that
becomes apparent in an internal-only validations as the one considered in this
section.

The results above are part of a more general phenomenon. As shown in [LBLB21],
several heuristics inherit biases present in up-to-date PPI networks and
manually curated ontologies such as GO, KEGG, and Reactome, in some cases with
results that improve when the PPI is replaced by a randomized one. To further
investigate this issue, in particular, the bias inherited by BRW when
ontologies from manually curated data sources are used, we considered the AMIM
(Active Module Identification Methods) test suite proposed in [LBLB21], which
allows systematic comparison of candidate gene sets predicted using an
original PPI network and perturbed versions thereof. To this purpose, we
considered the Lung Cancer phenotype, which was both considered in Section 3.2
and is used as a benchmark in [LBLB21]. We considered randomized versions of
the original PPI obtained using all randomization algorithms implemented in
the test suite (i.e., expected degree, uniform, shuffled and scale-free).
Candidate genes ranked by BRW, RWR-M, and RWR using the original PPI as input
were compared with the results obtained using each of the aforementioned
randomized counterparts of the PPI. Statistical significance was computed on
KEGG pathways enrichments computed by the test suite on predefined KEGG
pathway.

Fig. 2 f) shows that BRW, RWR-M, and RWR extract statistically significant
candidate gene sets on the original PPI network. Relying only on the PPI, RWR
provides no statistically significant results when the PPI is randomized,
except when expected degree is preserved.101010This phenomenon arises because
the stationary distribution of a random walk on an undirected network follows
exactly the degree distribution. As a result, the stationary distribution of a
random walk with restart over an undirected network is positively correlated
with degree distribution. On the other hand, BRW and RWR-M inherit PPI and
ontology’s biases, so that their results are still statistically significant
when the PPI is randomized using expected degree, scale-free and uniform
randomization algorithms. Still, BRW’s results are considerably more
significant when the original network is used, indicating that topology of the
PPI is crucial to (more) effectively propagate information extracted from
other biological sources. BRW predicts more significant outcomes on the actual
PPI than on the randomized versions. While all ontologies (GO, KEGG, and
Reactome) were considered in the experiment summarized in Fig. 2 f), we
further investigated BRW’s behaviour when KEGG (which is used in the
enrichment) is not used. In this case, BRW’s results are no longer
statistically significant if one randomizes the network, with the exception of
a mild significance when expected degree is preserved (see again footnote 10
for an explanation).

###  3.4 A Case Study: Breast Cancer Phenotype

In this section, we discuss the results of an in-depth analysis of BRW’s
performance on the Breast Cancer phenotype, a global health concern [SMJ+14,
CRA+21], with 284,200 new cases and more than 44,000 deaths in the USA in 2021
[SMJ+14]. In particular, we present i) an external validation using drugs FDA
approved drugs; ii) Drug enrichment; iii) an assessment of the algorithm(s)
stability across multiple populations.

####  3.4.1 Drug Target Discovery

We validated the top candidate gene sets prioritized by each algorithm along
different axes. To this purpose, we created a test set of target genes for
Breast Cancer drugs approved by the FDA as described in Section 3.1. Results
are summarized in Figure 3 a), showing Recall@K for the algorithms we
considered. Compared to other baselines, BRW prioritizes the highest number
(25%) of drug targets in the top 200 candidates. Looking at the specific gene
targets predicted by the algorithms, BRW predicts the highest number of genes
(14 out of the 20 predicted genes), making four unique predictions, namely,
Cyclin Dependent Kinases 4 and 6 (CDK4, CDK6), the Protein Kinase C Zeta
(PRKCZ) and Caspase 3 (CASP3). On the other hand, the DNA Topoisomerase II
Alpha (TOP2A) and the progesterone receptor (PGR) genes are uniquely predicted
by DIAMOnD. The Protein Kinase C Theta gene (PRKCQ) is only predicted by RWR-M
algorithm. More details are given in Supplementary Table 5.

Furthermore, we validated the prioritized genes from a different perspective:
we considered the group of drug targets returned by each framework and showed
the drugs that target them. We filtered out drugs in Drug Bank that were not
annotated with the ”Breast Cancer” or related associated condition.
Supplementary Table 3 shows how the number of drugs ranges from 6 to 17 as we
change the combination of hyper-parameters. As expected, the number of drugs
correlates positively with the predicted drug target percentage (Recall@200).
Supplementary Table 4 shows the drugs prioritized (i.e., a drug with at least
one drug target prioritized in the top 200 positions) by BRW and the other
algorithms. BRW prioritizes the highest number of drugs. Interestingly, CDK4/6
inhibitors palbociclib, ribociclib, and abemaciclib, currently used to treat
hormone receptor-positive/HER2-negative metastatic breast cancer [DFF+21], are
only predicted by BRW.

Finally, we identified drugs that have an enrichment with the top candidates
genes predicted by each algorithm. To identify enriched drugs, we used the
gseapy package for gene set enrichment analysis [STM+05], and we chose drugs
that were enriched by at least one algorithm with a corrected p-value lower
than $10^{-5}$. In particular, Supplementary Table 3 shows how the enrichment
is affected by various hyper-parameters combinations. In this case, drug
enrichment correlates positively with gene expression, and the best
combination is obtained when $\alpha$ and $\beta$ are equal to $0.25$. Indeed,
while drugs target and inhibit genes involved in disease-specific pathways,
the effect of the drugs can be measured by the differential expression and the
co-expression between targets and nearby genes [CMP+17]. In general, results
highlight the following trends: i) different biological sources (in our case,
gene expression and ontologies) provide complementary information, with
different subsets of FDA approved drug targets significantly enriched for both
low (gene expression bias) and high (ontologies bias) values of the parameters
$\alpha$ and $\beta$; ii) best performance is achieved for lower values of the
restart probability $r$ ($0.25$), confirming that BRW’s prioritization depends
on information that is not necessarily confined to the immediate neighbourhood
of the seed set.

Fig. 3 b) highlights drugs that are enriched in BRW and in the other baselines
considered in this manuscript, with details for each drug reported in
Supplementary Table 6. Interestingly, 11 FDA-approved drugs for breast cancer
treatment are enriched in BRW gene candidates (i.e., fulvestrant, doxorubicin,
paclitaxel, tamoxifen, methotrexate, letrozole, cyclophosphamide, trastuzumab,
fluorouracil, gemcitabine). The mTOR inhibitor sirolimus showed the best
significance. While it is not currently used to treat breast cancer,
preclinical in vivo studies demonstrated its potent antiangiogenic activity on
breast cancer models [MSTDJ+22]

####  3.4.2 Breast Cancer - Multi Population Study

A desirable property of disease gene prioritization should be a certain
stability in the set of proposed disease gene candidates across different
populations. In other words, to some extent (e.g., up to intrinsic biases or
qualitative differences in the datasets used), results should not overly
depend on the specific dataset the algorithm is analyzing. To investigate this
aspect, we characterized the behaviors of BRW and the other baselines we
considered when applied to different population studies on Breast Cancer. We
considered Invasive Breast Cancer population studies retrieved from cBioPortal
datasets [STD+12, SRG+12, BCRE+12, KDK+18, CGB+15] (see Section A of
supplementary material), which allowed us to identify five different seed sets
for the algorithms we considered, one per study. Figure 3 d) plots the number
of candidate genes, predicted by each algorithm, that are seeds in at least j
other studies (frequently mutated genes in the associated populations, i.e.
mutation frequency ¿ 1%), for $j\in\\{1,2,3,4\\}$. Notably, BRW has the
highest number of candidates that are frequently mutated in the other 4
studies (the far right blue bar). The 10 genes identified by BRW, that are
frequently mutated in 4 out of the 5 breast cancer populations are AR, JUN,
STAT3, NOTHC1, JAK2, HDAC1, SMAD4, YAP1, and CHD4. While the Histone
Deacetylase 1 (HDAC1) is also retrieved by all the other algorithms, SMAD4 and
CHD4 are only predicted by BRW. The remaining genes are returned by BRW and at
least another heuristic. Details are reported in Supplementary Table 7.

##  4 Discussion

Guided by the hypothesis that disease causing genes often share important
common biological processes and pathways, we extended the random walk with
restart approach to disease gene prioritization, proposing a framework that
allows seamless integration of multiple biological information sources. The
proposed approach consists of two main steps: i) extracting significant
disease features from disease-term association data, such as statistically
significant biological processes and pathways, and ii) using these features to
bias the random walk with restart in a way that is consistent with the
biological sources used. These two aspects are discussed in Section 2.

In general, BRW outperforms, in terms of standard indices of predictive
accuracy such as recall and nDCG, previous frameworks that only rely on a s
single biological source, such as Random Walks with Restart [NK10], DaDa
[EBEK11] and DIAMOnD [GMB15] that rely on Protein-Protein Interaction network.
This is also true for an extension of RWR, namely, RWR-M [VTN+19], that
performs a random walk on a multi-layer network. Using a Monte Carlo random
sampling validation, we showed that prioritization results returned by BRW
frequently outperform other baselines on four different cancer types: Breast,
Colorectal, Lung, and Thyroid cancer. These results were further supported by
a broader, computational validation on a corpus of seventy Mendelian disease
manually curated by [GMB15].

A main aim of precision medicine is to use disease genes to enable tailored
treatments. In this perspective, we investigated how the candidate genes
prioritized by each framework are related to Breast Cancer drug targets.
Results show that BRW prioritizes the highest number of drug targets in the
top 200 candidates (Figure 3a). As a further assessment, we considered the
group of drug targets returned by each algorithm and identified the drugs that
target them. We found that BRW prioritizes the highest number of FDA-approved
drugs for Breast Cancer treatment (Supplementary Table 4). In general, we
noticed that drug targets are more correlated with phenotypic pathways, while
drug enrichment highlights how drugs affect gene expression in terms of co-
expression and differential expression. This is not surprising: indeed,
although drugs often target and inhibit genes involved in disease-specific
pathways, the effect of the drugs can be measured by the differential
expression and the co-expression between targets and nearby genes [CMP+17].

To investigate the stability of the proposed disease gene candidates, we
characterized the behaviors of BRW and the previous frameworks across
different population studies. We selected “invasive breast cancer” as a
phenotype and we retrieved data for 5 different populations from cBioPortal
(https://www.cbioportal.org/). Consistently, BRW showed excellent stability,
with the highest number of gene candidates in one study that are frequently
mutated in the other four (Figure 3 c), Supplementary Table 7).

The BRW framework is not without limitations. On one hand, inducing a bias in
teleporting probability and the transition matrix through the use of
ontologies improves predictive accuracy of the algorithm. On the other hand,
this bias can hinder BRW’s ability to identify new disease-related pathways.
For this reason, we believe it is important to exploit the framework’s ability
to integrate heterogeneous, hopefully complementary, data sources.
Furthermore, experiments summarized in Fig. 7 and Supplementary Table 3
quantitatively showed that every data source comes with its bias. As a result,
the relative weights attributed to different biological information sources
(the combined choice of $\alpha$ and $\beta$ in our case) can significantly
affect predictive accuracy, with a magnitude that in general depends on the
validation test used. For example, while ontologies showed most effective in a
computational validation and for drug target discovery, gene expression proved
particularly useful in identifying a candidate gene set that is highly
enriched in breast cancer and cancer related drugs. In our opinion, these
results provide support to the use of (at least partially) complementary
sources of biological information, even though these sometimes present non-
negligible correlations, as discussed in Section 3.3.

We also emphasize that, in this study, we only leveraged a limited set of
disease-specific data sources (co-expression and differential expression).
Potential performance improvements might be achieved by integrating further
disease information sources, such as methylation data, microRNA expression, or
microRNA-target gene associations. Doing this might provide further insights
into key biological mechanisms and provide new prospective gene targets,
though, of course, only functional studies can provide the ultimate answer as
to their biological role.

##  5 Funding.

Partially supported by the ERC Advanced Grant 788893 AMDROMA “Algorithmic and
Mechanism Design Research in Online Markets”, the EC H2020RIA project
“SoBigData++” (871042), and the MIUR PRIN project ALGADIMAR “Algorithms,
Games, and Digital Markets”.

## References

  * [ABG+09] I. Armean, A. Bridge, A. T. Ghanbarian, B. Aranda, B. Roechert, C. Derow, C. Leroy, H. Hermjakob, J. Kerssemakers, J. Khadake, K. van Eijk, L. Montecchi-Palazzi, M. Feuermann, M. Menden, M. Michaut, P. Achuthan, S. Kerrien, S. Orchard, S. N. Neuhauser, V. Perreau, and Y. Alam-Faruque.  The IntAct molecular interaction database in 2010.  Nucleic Acids Research, 38(1):D525–D531, 2009.
  * [ABSH08] Joanna Amberger, Carol A Bocchini, Alan F Scott, and Ada Hamosh.  Mckusick’s online mendelian inheritance in man (omim®).  Nucleic acids research, 37(1):D793–D796, 2008.
  * [ALANS16] Gregorio Alanis-Lobato, Miguel A Andrade-Navarro, and Martin H Schaefer.  Hippie v2. 0: enhancing meaningfulness and reliability of protein–protein interaction networks.  Nucleic acids research, page gkw985, 2016.
  * [BCRE+12] Shantanu Banerji, Kristian Cibulskis, Claudia Rangel-Escareno, Kristin K Brown, Scott L Carter, Abbie M Frederick, Michael S Lawrence, Andrey Y Sivachenko, Carrie Sougnez, Lihua Zou, et al.  Sequence analysis of mutations and translocations across breast cancer subtypes.  Nature, 486(7403):405–409, 2012.
  * [BGL11] Albert-László Barabási, Natali Gulbahce, and Joseph Loscalzo.  Network medicine: a network-based approach to human disease.  Nature reviews genetics, 12(1):56, 2011.
  * [CACP+09] Andrew Chatr Aryamontri, Arnaud Ceol, Daniele Peluso, Gianni Cesareni, Leonardo Briganti, Livia Perfetto, Luana Licata, and Luisa Castagnoli.  Mint, the molecular interaction database: 2009 update.  Nucleic Acids Research, 38(1):D532–D539, 2009.
  * [CGB+15] Giovanni Ciriello, Michael L Gatza, Andrew H Beck, Matthew D Wilkerson, Suhn K Rhie, Alessandro Pastore, Hailei Zhang, Michael McLellan, Christina Yau, Cyriac Kandoth, et al.  Comprehensive molecular portraits of invasive lobular breast cancer.  Cell, 163(2):506–519, 2015.
  * [CGP+17] Debyani Chakravarty, Jianjiong Gao, Sarah Phillips, Ritika Kundra, Hongxin Zhang, Jiaojiao Wang, Julia E Rudolph, Rona Yaeger, Tara Soumerai, Moriah H Nissan, et al.  Oncokb: a precision oncology knowledge base.  JCO precision oncology, 1:1–16, 2017.
  * [CIRS17] 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.
  * [CMP+17] Bin Chen, Li Ma, Hyojung Paik, Marina Sirota, Wei Wei, Mei-Sze Chua, Samuel So, and Atul J Butte.  Reversal of cancer gene expression correlates with drug efficacy and reveals therapeutic targets.  Nature communications, 8(1):1–12, 2017.
  * [Con19] The Gene Ontology Consortium.  The gene ontology resource: 20 years and still going strong.  Nucleic Acids Research, 47(Database-Issue):D330–D338, 2019.
  * [CRA+21] Edoardo Crimini, Matteo Repetto, Philippe Aftimos, Andrea Botticelli, Paolo Marchetti, and Giuseppe Curigliano.  Precision medicine in breast cancer: From clinical trials to clinical practice.  Cancer Treatment Reviews, 98, 2021.
  * [CS04] Bond AT Carter SL, Griffin M.  Gene co-expression network topology provides a framework for molecular characterization of cellular state.  Bioinformatics (Oxford, England), pages 2242–50, 09 2004.
  * [CWB+10] Andrew Chatraryamontri, Andrew Winter, Bobby-Joe Breitkreutz, Chris Stark, Jennifer M. Rust, Julie Nixon, Kara Dolinski, Kimberly Van Auken, Lorrie Boucher, Michael S. Livstone, Mike Tyers, Rose Oughtred, Teresa Reguly, Xiaodong Wang, and Xiaoqi Shi.  The biogrid interaction database: 2011 update.  Nucleic Acids Research, 39(1):D698–D704, 2010.
  * [DBTVOM07] Tijl De Bie, Leon-Charles Tranchevent, Liesbeth MM Van Oeffelen, and Yves Moreau.  Kernel-based data fusion for gene prioritization.  Bioinformatics, 23(13):i125–i132, 2007.
  * [DFF+21] Simona Duranti, Alessandra Fabi, Marco Filetti, Rosa Falcone, Pasquale Lombardi, Gennaro Daniele, Gianluca Franceschini, Luisa Carbognin, Antonella Palazzo, Giorgia Garganese, et al.  Breast cancer drug approvals issued by ema: A review of clinical trials.  Cancers, 13(20):5198, 2021.
  * [DGB07] Werner Dubitzky, Martin Granzow, and Daniel P Berrar.  Fundamentals of data mining in genomics and proteomics.  Springer Science & Business Media, 2007.
  * [DHH+18] Christos Dimitrakopoulos, Sravanth Kumar Hindupur, Luca Häfliger, Jonas Behr, Hesam Montazeri, Michael N Hall, and Niko Beerenwinkel.  Network-based integration of multi-omics data for prioritizing cancer genes.  Bioinformatics, 34(14):2441–2448, 2018.
  * [ea08] Abhilash et al.  Human protein reference database—2009 update.  Nucleic Acids Research, 37(1):D767–D772, 2008.
  * [EBEK11] Sinan Erten, Gurkan Bebek, Rob M Ewing, and Mehmet Koyutürk.  D a d a: Degree-aware algorithms for network-based disease gene prioritization.  BioData mining, 4(1):19, 2011.
  * [GMB15] Susan Dina Ghiassian, Jörg Menche, and Albert-László Barabási.  A disease module detection (diamond) algorithm derived from a systematic analysis of connectivity patterns of disease proteins in the human interactome.  PLoS computational biology, 11(4):e1004120, 2015.
  * [GMP+19] Michele Gentili, Leonardo Martini, Manuela Petti, Lorenzo Farina, and Luca Becchetti.  Biological random walks: Integrating heterogeneous data in disease gene prioritization.  In 2019 IEEE Conference on Computational Intelligence in Bioinformatics and Computational Biology (CIBCB), pages 1–8. IEEE, 2019.
  * [GNZ+14] Mika Gustafsson, Colm E Nestor, Huan Zhang, Albert-László Barabási, Sergio Baranzini, Sören Brunak, Kian Fan Chung, Howard J Federoff, Anne-Claude Gavin, Richard R Meehan, et al.  Modules, networks and systems medicine for understanding disease and aiding diagnosis.  Genome medicine, 6(10):82, 2014.
  * [HS09] John Hardy and Andrew Singleton.  Genomewide association studies and human disease.  New England Journal of Medicine, 360(17):1759–1768, 2009.
  * [JMV+19] Bijay Jassal, Lisa Matthews, Guilherme Viteri, Chuqiao Gong, Pascual Lorente, Antonio Fabregat, Konstantinos Sidiropoulos, Justin Cook, Marc Gillespie, Robin Haw, Fred Loney, Bruce May, Marija Milacic, Karen Rothfels, Cristoffer Sevilla, Veronica Shamovsky, Solomon Shorser, Thawfeek Varusai, Joel Weiser, Guanming Wu, Lincoln Stein, Henning Hermjakob, and Peter D’Eustachio.  The reactome pathway knowledgebase.  Nucleic Acids Research, 48(D1):D498–D503, 11 2019.
  * [JW03] Glen Jeh and Jennifer Widom.  Scaling personalized web search.  In Proceedings of the 12th international conference on World Wide Web, pages 271–279, 2003.
  * [Kan19] Minoru Kanehisa.  Toward understanding the origin and evolution of cellular organisms.  Protein Science, 28(11):1947–1951, 2019.
  * [KBHR08] Sebastian Köhler, Sebastian Bauer, Denise Horn, and Peter N Robinson.  Walking the interactome for prioritization of candidate disease genes.  The American Journal of Human Genetics, 82(4):949–958, 2008.
  * [KDK+18] Zhengyan Kan, Ying Ding, Jinho Kim, Hae Hyun Jung, Woosung Chung, Samir Lal, Soonweng Cho, Julio Fernandez-Banet, Se Kyung Lee, Seok Won Kim, et al.  Multi-omics profiling of younger asian breast cancers reveals distinctive molecular signatures.  Nature communications, 9(1):1–13, 2018.
  * [KG00] Minoru Kanehisa and Susumu Goto.  KEGG: Kyoto Encyclopedia of Genes and Genomes.  Nucleic Acids Research, 28(1):27–30, 01 2000.
  * [KLPK+03] A. E. Kel, B. Lewicki-Potapov, D. Karas, D.-U. Kloos, E. Fricke, E. Gößling, E. Wingender, H. Michael, H. Saxel, I. Reuter, K. Hornischer, M. Haubrock, M. Scheer, O. V. Kel-Margoulis, R. Geffers, R. Hehl, R. Münch, S. Land, S. Rotert, S. Thiele, and V. Matys.  Transfac ® : transcriptional regulation, from patterns to profiles.  Nucleic Acids Research, 31(1):374–378, 2003.
  * [KSF+18] Minoru Kanehisa, Yoko Sato, Miho Furumichi, Kanae Morishima, and Mao Tanabe.  New approach for understanding genome variations in KEGG.  Nucleic Acids Research, 47(D1):D590–D595, 10 2018.
  * [KSL+13] Daniel C. Koboldt, Karyn Meltz Steinberg, David E. Larson, Richard K. Wilson, and Elaine R. Mardis.  The next-generation sequencing revolution and its impact on genomics.  Cell, 155(1):27 – 38, 2013.
  * [KSSF17] Tamas Korcsmaros, Maria Victoria Schneider, and Giulio Superti-Furga.  Next generation of network medicine: interdisciplinary signaling approaches.  Integrative Biology, 9(2):97–108, 2017.
  * [LBLB21] Olga Lazareva, Jan Baumbach, Markus List, and David B Blumenthal.  On the limits of active module identification.  Briefings in Bioinformatics, 22(5):bbab066, 2021.
  * [LP17] David A Levin and Yuval Peres.  Markov chains and mixing times, volume 107.  American Mathematical Soc., 2017.
  * [MGS+17] Jörg Menche, Emre Guney, Amitabh Sharma, Patrick J Branigan, Matthew J Loza, Frédéric Baribaud, Radu Dobrin, and Albert-László Barabási.  Integrating personalized gene expression profiles into predictive disease-associated gene pools.  NPJ systems biology and applications, 3(1):1–10, 2017.
  * [MKI+03] Ron Milo, Nadav Kashtan, Shalev Itzkovitz, Mark EJ Newman, and Uri Alon.  On the uniform generation of random graphs with prescribed degree sequences.  arXiv preprint cond-mat/0312028, 2003.
  * [MSK+15] Jörg Menche, Amitabh Sharma, Maksim Kitsak, Susan Dina Ghiassian, Marc Vidal, Joseph Loscalzo, and Albert-László Barabási.  Uncovering disease-disease relationships through the incomplete interactome.  Science, 347(6224):1257601, 2015.
  * [MSTDJ+22] Muhammad Shahidan Muhammad Sakri, Tengku Ahmad Damitri Al-Astani Tengku Din, Hasnan Jaafar, Vinod Gopalan, and Wan Faiziah Wan Abdul Rahman.  Rapamycin as a potent and selective inhibitor of vascular endothelial growth factor receptor in breast carcinoma.  International Journal of Immunopathology and Pharmacology, 36:20587384211059673, 2022.
  * [NK10] Saket Navlakha and Carl Kingsford.  The power of protein interaction networks for associating genes with diseases.  Bioinformatics, 26(8):1057–1063, 2010.
  * [ODC+18] Kivilcim Ozturk, Michelle Dow, Daniel E Carlin, Rafael Bejar, and Hannah Carter.  The emerging potential for network analysis to inform precision cancer medicine.  Journal of molecular biology, 2018.
  * [PRASP+20] Janet Piñero, Juan Manuel Ramírez-Anguita, Josep Saüch-Pitarch, Francesco Ronzano, Emilio Centeno, Ferran Sanz, and Laura I Furlong.  The disgenet knowledge platform for disease genomics: 2019 update.  Nucleic acids research, 48(D1):D845–D855, 2020.
  * [RBW+09] Andreas Ruepp, Barbara Brauner, Brigitte Waegele, Corinna Montrone, Gisela Fobo, Goar Frishman, H.-Werner Mewes, Irmtraud Dunger-Kaltenbach, and Martin Lechner.  CORUM: the comprehensive resource of mammalian protein complexes—2009.  Nucleic Acids Research, 38(1):D497–D501, 2009.
  * [RHJ+14] Erin M Ramos, Douglas Hoffman, Heather A Junkins, Donna Maglott, Lon Phan, Stephen T Sherry, Mike Feolo, and Lucia A Hindorff.  Phenotype–genotype integrator (phegeni): synthesizing genome-wide association study (gwas) data with existing genomic resources.  European Journal of Human Genetics, 01 2014.
  * [SCCR12] Gabriele Sales, Enrica Calura, Duccio Cavalieri, and Chiara Romualdi.  g raphite-a bioconductor package to convert pathway topology to gene network.  BMC bioinformatics, 13(1):1–12, 2012.
  * [SL20] Haixia Shang and Zhi-Ping Liu.  Network-based prioritization of cancer genes by integrative ranks from multi-omics data.  Computers in biology and medicine, 119:103692, 2020.
  * [SMJ+14] Rebecca Siegal, KD Miller, Ahmddin Jemal, et al.  Cancer statistics, 2012.  Ca Cancer J Clin, 64(1):9–29, 2014.
  * [SPCP10] Jan Schellenberger, Junyoung O. Park, Tom M. Conrad, and Bernhard Palsson.  Bigg: a biochemical genetic and genomic knowledgebase of large scale metabolic reconstructions.  BMC Bioinformatics, 11(1):213, 2010.
  * [SRG+12] Sohrab P Shah, Andrew Roth, Rodrigo Goya, Arusha Oloumi, Gavin Ha, Yongjun Zhao, Gulisa Turashvili, Jiarui Ding, Kane Tse, Gholamreza Haffari, et al.  The clonal and mutational evolution spectrum of primary triple-negative breast cancers.  Nature, 486(7403):395–399, 2012.
  * [SSA+20] Edwin K Silverman, Harald HHW Schmidt, Eleni Anastasiadou, Lucia Altucci, Marco Angelini, Lina Badimon, Jean-Luc Balligand, Giuditta Benincasa, Giovambattista Capasso, Federica Conte, et al.  Molecular networks in network medicine: Development and applications.  Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 12(6):e1489, 2020.
  * [STD+12] PJ Stephens, PS Tarpey, H Davies, P Van Loo, C Greenman, DC Wedge, S Nik-Zainal, S Martin, I Varela, GR Bignell, et al.  Oslo breast cancer consortium (osbreac) the landscape of cancer genes and mutational processes in breast cancer.  Nature, 486(7403):400–404, 2012.
  * [STM+05] Aravind Subramanian, Pablo Tamayo, Vamsi K. Mootha, Sayan Mukherjee, Benjamin L. Ebert, Michael A. Gillette, Amanda Paulovich, Scott L. Pomeroy, Todd R. Golub, Eric S. Lander, and Jill P. Mesirov.  Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles.  Proceedings of the National Academy of Sciences, 102(43):15545–15550, 2005.
  * [VTN+19] Alberto Valdeolivas, Laurent Tichit, Claire Navarro, Sophie Perrin, Gaelle Odelin, Nicolas Levy, Pierre Cau, Elisabeth Remy, and Anaïs Baudot.  Random walk with restart on multiplex and heterogeneous biological networks.  Bioinformatics, 35(3):497–505, 2019.
  * [WDG+17] David Wishart, Yannick Djoumbou, An Chi Guo, Elvis Lo, Ana Marcu, Jason Grant, Tanvir Sajed, Daniel Johnson, Carin Li, Zinat Sayeeda, Nazanin Assempour, Ithayavani Iynkkaran, Yifeng Liu, Adam Maciejewski, Nicola Gale, Alex Wilson, Lucy Chin, Ryan Cummings, Diana Le, and Michael Wilson.  Drugbank 5.0: A major update to the drugbank database for 2018.  Nucleic acids research, 46, 11 2017.
  * [XL06] Jianzhen Xu and Yongjin Li.  Discovering disease-genes by topological features in human protein–protein interaction network.  Bioinformatics, 22(22):2800–2805, 09 2006.

##  Appendix A Data

###  A.1 Protein-Protein Interaction networks

To run the experiment presented in Section 3, we considered two different PPI
networks:

  1. 1.

HIPPIE-v2.2: Protein-protein interaction network downloaded from [ALANS16]. We
filtered the network, keeping only the edges with a score larger than 0.7.
Then, we mapped each node to its Ensembl id, removing nodes mapped on multiple
Ensembl Ids. This way, we derive a PPI network with 88,865 edges and 12,148
nodes.

  2. 2.

[GMB15]: We used the same Protein-Protein Interaction network (PPI) as in
[GMB15]. In [GMB15], the authors only considered direct physical protein
interactions with reported experimental evidence. They derive this PPI network
with the help of several data sources: TRANSFAC, [KLPK+03], IntAct [ABG+09],
MINT [CACP+09], BioGRID [CWB+10], HPRD [ea08] KEGG, BIGG [SPCP10], and CORUM
[RBW+09]. Then, we considered the main connected component of the network and
we removed self-loops (i.e., edges describing proteins’ self-interactions).
The resulting graph consists of 13396 nodes and 138405 edges.

##### Remark.

Note that the PPI of [GMB15] was used only in the internal validation on the
70 mendelian diseases, whose results are discussed in Section 3.2 of the main
article. This was done for consistency, since the network considered in
[GMB15] was explicitly curated to include all genes that are associated to at
least one of the 70 diseases under consideration, some of which might be
missing in other, even more recent, networks. Instead, all other experiments
in this paper used the more recent HIPPIE-v2.2 PPI [ALANS16].

###  A.2 Biological data sources

##### Gene annotations.

We used three different sources of gene biological information: Gene Ontology
Consortium111111http://geneontology.org/. where, for each gene, we downloaded
its biological processes, KEGG Pathways [Kan19, KSF+18, KG00] and Reactome
[JMV+19], which we used to download pathways’ annotations.

##### Pathways Network:

a network connecting proteins according to pathway interaction data, using the
R package graphite [SCCR12].

##### Tissue Network:

a co-expression network from RNA-Seq data downloaded from the Human Protein
Atlas (http://www.proteinatlas.org) (Uhlen et al., 2015). We computed Spearman
correlations of TPM (transcript per million) expression data from 32 tissues
and 45 cell lines, and selected correlation having an absolute value $\geq$
0.7 to be included in the network.

##### Gene expression data (TCGA).

Gene expression datasets were downloaded from The Cancer Genome Atlas
datasets.121212https://portal.gdc.cancer.gov/. We collected from CDC portal
the RNA-Seq data of several type of cancers. Table 1 describe the cancer types
downloaded and the number of Primary Tumor and Solid tissue normal used to
compute the set of differentially expressed genes,and co-expression network.

Table 1: The table shows all the TCGA projects downloaded to create the set of differentially expressed genes and the CO-expression network Project ID | Cancer Type | N. of Primary Tumor | N. of Solid Tissue Normal
---|---|---|---
TCGA-BRCA | Breast Cancer | 1077 | 99
TCGA-LUAD | Lung Adenocarcinoma | 519 | 58
TCGA-THCA | Papillary Thyroid Cancer | 497 | 56
TCGA-COAD | Colorectal Adenocarcinoma | 456 | 41

###  A.3 Data sources for experimental validation

##### Disease-gene associations.

We used disease-gene associations from [GMB15]. Out of a corpus of gene-
disease associations involving 69 well-characterized diseases retrieved from
OMIM131313https://www.ncbi.nlm.nih.gov/omim/ (Online Mendelian Inheritance in
Man; [ABSH08]) and PheGenI141414https://www.ncbi.nlm.nih.gov/gap/phegeni
(Phenotype-Genotype Integrator; [RHJ+14]) databases, we focused on the Breast
Cancer phenotype, which has associations with 40 genes.

##### Disease-Gene associations of Invasive Breast Cancer on several
population studies.

From cBioPortal151515https://www.cbioportal.org/ we downloaded Invasive Breast
Cancer mutated genes belonging to 5 studies [STD+12, SRG+12, BCRE+12, KDK+18,
CGB+15] (Some information about the subjects sequenced by each study is
reported in Supplementary Fig. 1) and we filtered out genes mutated in less
than 1% of the studied population and that are not present in the OncoKB
[CGP+17]. More precisely, for each study (_reference population_), a reliable
seed set was obtained, by keeping genes that i) were highly mutated in the
population study under consideration (mutation frequency $>1$%) and ii) were
known cancer genes, present in OncoKB database [CGP+17]. At the same time,
significantly mutated genes identified across the remaining four studies
provided a natural validation set. Thus, for each algorithm and for each
reference population, a candidate gene prioritized by the algorithm is
considered to be a potential disease gene if it appears mutated gene in at
least one of the other population studies. The underlying hypothesis is that
genes involved in a pathological pathway are likely to be mutated in at least
one population.

Table 2: Breast cancer drugs filtered from DrugBank to derive a set of target genes we used to validate algorithm prediction on an external dataset.  Drug Bank ID |  Drug Name |  Target Genes |  Drug Bank ID |  Drug Name |  Target Genes
---|---|---|---|---|---
DB12001 |  Abemaciclib |  CDK4, CDK6 |  DB05773 |  Trastuzumab emtansine |  ERBB2
DB01229 |  Paclitaxel |  TUBB1, BCL2, MAP4, MAP2, MAPT, NR1I2 |  DB09037 |  Pembrolizumab |  PDCD1
DB05773 |  Trastuzumab emtansine |  ERBB2 |  DB11730 |  Ribociclib |  CDK4, CDK6
DB01590 |  Everolimus |  MTOR |  DB01259 |  Lapatinib |  EGFR, ERBB2
DB12015 |  Alpelisib |  PIK3CA |  DB01006 |  Letrozole |  CYP19A1
DB01217 |  Anastrozole |  CYP19A1 |  DB09074 |  Olaparib |  PARP1, PARP2, PARP3
DB00282 |  Pamidronic acid (Aredia) |  FDPS, GGPS1, CASP3, CASP9 |  DB00351 |  Megestrol acetate |  PGR, NR3C1
DB01101 |  Capecitabine |  TYMS |  DB00563 |  Methotrexate |  DHFR, TYMS, ATIC
DB00531 |  Cyclophosphamide |  NR1I2 |  DB11828 |  Neratinib |  EGFR
DB01248 |  Docetaxel |  TUBB1, BCL2, MAP4, MAP2, MAPT, NR1I2 |  DB09074 |  Olaparib |  PARP1, PARP2, PARP3
DB00997 |  Doxorubicin |  TOP2A, NOLC1 |  DB01229 |  Paclitaxel |  TUBB1, BCL2, MAP4, MAP2, MAPT, NR1I2
DB00445 |  Epirubicin |  TOP2A |  DB09073 |  Palbociclib |  CDK4, CDK6
DB14962 |  Trastuzumab deruxtecan |  FCGR1A, TOP1 |  DB00282 |  Pamidronic acid |  FDPS, GGPS1, CASP3, CASP9
DB00445 |  Epirubicin |  TOP2A |  DB09037 |  Pembrolizumab |  PDCD1
DB08871 |  Eribulin |  BCL2, TUBB1 |  DB06366 |  Pertuzumab |  ERBB2
DB01590 |  Everolimus |  MTOR |  DB14740 |  Hyaluronidase |  TGFB1
DB00990 |  Exemestane |  CYP19A1 |  DB12015 |  Alpelisib |  PIK3CA
DB00544 |  Fluorouracil |  TYMS |  DB12893 |  Sacituzumab govitecan |  TACSTD2, TOP1, FUBP1
DB14962 |  Trastuzumab deruxtecan |  FCGR1A, TOP1 |  DB00675 |  Tamoxifen |  ESR1, ESR2, EBP, AR, KCNH2, NR1I2, ESRRG, SHBG, MAPK8, PRKCA, PRKCB, PRKCD, PRKCE, PRKCG, PRKCI, PRKCQ, PRKCZ
DB00539 |  Toremifene |  ESR1, SHBG |  DB11760 |  Talazoparib |  PARP1, PARP2
DB00947 |  Fulvestrant |  ESR1 |  DB11595 |  Atezolizumab |  CD274
DB01006 |  Letrozole |  CYP19A1 |  DB00539 |  Toremifene |  ESR1, SHBG
DB00544 |  Fluorouracil |  TYMS |  DB00563 |  Methotrexate |  DHFR, ATIC, TYMS
DB00441 |  Gemcitabine |  RRM1, TYMS, CMPK1 |  DB12893 |  Sacituzumab govitecan |  TACSTD2, TOP1, FUBP1
DB00014 |  Goserelin |  LHCGR, GNRHR |  DB11652 |  Tucatinib |  ERBB2, ERBB3
DB08871 |  Eribulin |  BCL2, TUBB1 |  DB01259 |  Lapatinib |  EGFR, ERBB2
DB00072 |  Trastuzumab |  ERBB2 |  DB12001 |  Abemaciclib |  CDK4, CDK6
DB09073 |  Palbociclib |  CDK4, CDK6 |  DB00570 |  Vinblastine |  TUBA1A, TUBB, TUBD1, TUBG1, TUBE1, JUN
DB00441 |  Gemcitabine |  RRM1, TYMS, CMPK1 |  DB01101 |  Capecitabine |  TYMS
DB04845 |  Ixabepilone |  TUBB3 |  DB00014 |  Goserelin |  LHCGR, GNRHR
Supplementary Fig. 1: Breast cancer studies selected from the cBioPortal for
Cancer Genomics. The table reports some details on the invasive Breast Cancer
populations used in the experiments. A study was selected if it included human
participants AND at least a whole exome sequencing analysis was performed
(studies performing targeted sequencing analyses were ruled out) ER: Estrogen
receptor, WES: Whole Exome Sequencing, WGS: Whole Genome Sequencing, SNVs:
Single Neucleotide Variants, CNA: Copy Number Alteration, GF: Gene Fusions,
GEA: Gene Expression Alteration

##  Appendix B Methods

###  B.1 Gene co-expression networks

At a high level, a Gene Co-expression Network (GCN) is an undirected graph,
with genes as the vertices. An undirected edge exists between two nodes $i$
and $j$ if the corresponding genes are significantly co-expressed in the
subject population.161616In the remainder we use the same letter to denote a
gene and the corresponding vertex in the PPI or GCN. Hence, differently from a
PPI, a GCN is disease-specific. Intuitively, for two genes $i$ and $v$, this
means that the transcript levels of $u$ and $j$ tend to vary in a similar
fashion across the subject population. Gene co-expression networks [CS04] are
of biological interest since co-expressed genes tend to be controlled by the
same transcriptional regulatory program, i.e., they are functionally related
or members of the same pathway or protein complex.

###  B.2 Biological Random Walks

####  B.2.1 Overview

As we decribed in the main article, we compute the personalization vector
$\mathbf{q}$ and the transition matrix $W$ by leveraging multiple biological
information sources. We remind that we refer to the biological information
associated to a gene i (e.g., the set of its annotations) as the set of its
features, denoted by features(i), such as the Gene Ontology database [Con19]
(GO in the remainder) or gene expression levels and their consequently
differential expression scores. BRW ranks genes according to three main steps:

  1. 1.

Rather than using the standard approach171717We remind that in the standard
RWR approach [KBHR08], the probability of restarting the random walk from a
given seed node (disease gene) is the same for all seeds nodes, while it is
$0$ for other nodes of the PPI., we compute a personalization vector
$\mathbf{q}$, i.e Biological Teleporting Probability (BTP), associating
individual teleporting probabilities to all nodes of the PPI. Qualitatively
speaking, the probability of teleporting to a given node increases, the higher
its feature-based similarity with the disease.

  2. 2.

In a similar fashion, we derive the transition matrix $W$ of the random walk
by also leveraging biological information, whereby the probability of a
transition between genes $i$ and $j$ tends to be higher, the more functionally
closer the genes are.

  3. 3.

Finally, we rank genes according to their Biological Random Walk (BRW) scores.

###  B.3 Biological Teleporting Probability

Both annotation data and gene expression levels allow derivation of
personalization vectors that reflect some notion of similarity between a gene
and a disease, i.e. its seed genes. In the first case, this similarity is
defined in terms of knowledge about genes’ involvement in various biological
functions and/or diseases. In the second case, similarity is defined in terms
of co-expression levels of different genes in subject cases as opposed to
expression levels in a control group, using data about a population of
patients involved in a clinical trial. It’s worth noting that this approach
allows to easily extend the algorithm on a multiomics approach and/or tissue
specific analysis, as long as it is possible to recapitulate data as a score
for each gene (i.e.chromatin accessibility) or pair of genes (i.e. co-targeted
by miRNA).

Supplementary Fig. 2: Enrichment step: starting from a seed set $S$, We filter
out $S$’ annotations that are not statistically significant using Fisher Exact
Test and FDR correction.

####  B.3.1 Computing a Personalization Vector

Building on the hypothesis that genes involved in the same disease tend to be
functionally related and thus share similar biological information, a first
approach to deriving the personalization vector $\mathbf{q}$ relies on
assigning scores to genes in the PPI using biological information.

##### Using Biological Information.

We assume we have $\ell$ sources of biological information. Assume $S$ is the
seed set. for every $j=1,\ldots,\ell$, we use $\mathcal{F}^{j}$ to denote the
subset of annotations from the $j$-th source that are associated with at least
one gene in $S$. Then, for every $j=1,\ldots,\ell$, we selected the subset of
annotation $\widehat{\mathcal{F}^{j}}$ filtering out the annotations that are
not statistically significant as shown in Supplementary Fig.2 (i.e.,
$p-value>10^{-}5$ using Fisher Exact Test and FDR correction) so that
$\mathcal{F}=\cup_{j=1}^{\ell}\widehat{\mathcal{F}^{j}}$ denotes the set of
all statistically significant annotations that are associated with genes in
$S$. Likewise, for every gene $i$ (not necessarily belonging to $S$), we
denote by $\mathcal{A}(i)$ the set of its annotations, possibly extracted from
multiple biological information sources. We next assign each gene $i$ a weight
$\theta_{i}$ as follows:

$\theta_{i}=\begin{cases}\ell&\text{if }i\in S,\\\
\sum_{j=1}^{\ell}\frac{|\mathcal{A}(i)\cap\mathcal{F}^{j}|}{|\mathcal{F}^{j}|}&\text{otherwise},\end{cases}$
(8)

if we want to weigh the seed set $S$ more than the subset of genes that are
biologically similar to $S$. Otherwise, we can assign weights as follows:

$\theta_{i}=\sum_{j=1}^{\ell}\frac{|\mathcal{A}(i)\cap\mathcal{F}^{j}|}{|\mathcal{F}^{j}|}$
(9)

If no biological information is used (i.e., we are only using knowledge of
whether a given gene belongs to the seed set or not like in [KBHR08]), each
gene $i$ is assigned a weight $\theta_{i}$ as follows:

$\theta_{i}=\begin{cases}1&\text{if }i\in S\\\ 0&otherwise\end{cases}$ (10)

At this point, components of a personalization vector $\mathbf{q}$ can be
computed as follows:

$q_{i}=\frac{\theta_{i}}{\sum_{i=0}^{N}\theta_{i}},$ (11)

with $N$ the number of genes we consider. Note that $q_{i}$ denotes the
probability that, upon teleportation, the random walkers jumps to gene $i$.

##### Using gene expression data.

RNA sequencing techniques have changed our ability to explore the molecular
mechanisms underlying complex diseases and they have been used to identify
potential disease-associated genome-wide changes in gene expression patients.
An important goal is the identification of Differentially Expressed (DE)
genes, whose expression levels systematically differ between case (Breast
Cancer tissues) and a control group (Healthy Breast tissues).

We take inspiration from the framework to construct and integrate personalized
perturbation profiles (PEEPs) from gene expression data proposed in [MGS+17],
to derive gene expression-aware personalization vectors.

PEEPs turn individual expression heterogeneity into a predictive information
by constructing personalized perturbation profiles that reflect expression
changes within a single subject. For each gene $i$ and for each subject $j$,
the expression level $l_{i}^{j}$ is compared with the distribution of the
expression level of gene $i$ within the control group. This comparison is
measured by the corresponding z-score, defined as:

$z_{ij}=\frac{l_{i}^{j}-\mu_{i}}{\sigma_{i}},$ (12)

where $\mu_{i}$ and $\sigma_{i}$ respectively denote the mean and the standard
deviation of gene $i$’s expression level in the control group. This approach
allows association of a “bar code” [MGS+17] to each subject, thus allowing to
tell differentially expressed genes ($|z_{ij}|>2.5$) from regular ones.

In general, given a case and a control group, it is possible to build a matrix
$Z\in R^{n\times m}$, where $n$ and $m$ are respectively the number of genes
and subjects:

$Z=\begin{bmatrix}z_{11}&z_{12}&\cdots&z_{1m}\\\
z_{21}&z_{22}&\cdots&z_{2m}\\\ \vdots&\vdots&\ddots&\vdots\\\
z_{n1}&z_{n2}&\cdots&z_{nm}\end{bmatrix}$ (13)

Thus, we might consider only one gene and subject at time, or we might
consider a gene $i$ as the vector $Z_{i}$ of the $z$-scores of its expression
levels in the subject group. This provides an overall description of how gene
$i$’s expression varies within the case group.

Using matrix $Z$, we identify the subset of genes that are differentially
expressed in the subject group as follows. In particular, we define a new
binary matrix $\hat{Z}$ which, for each gene $i$ and patient $j$, specifies
whether or not $i$ is differentially expressed in subject $j$. Following
[MGS+17], we assume this to happen when its $z$-score $z_{ij}$ exceeds $2.5$
standard deviations. Namely:

$\hat{Z}_{ij}=\begin{cases}1&\text{if }|z_{ij}|>2.5\\\
0&otherwise\par\end{cases}$ (14)

Next, a gene $g$ is considered to be _differentially expressed_ in the subject
group if the number of patients in which it is differentially expressed
according to the above definition exceeds the mean:

$\sum_{j}\hat{Z}_{gj}>\frac{1}{n}\sum_{i}\sum_{j}\hat{Z_{i}^{j}}.$ (15)

The subset of differentially expressed genes this identified can be used to
define a personalization vector $\mathbf{q}$, using the PPI to identify
differentially expressed genes that are also close to known disease genes
(seeds) in the PPI. In particular, we obtain a personalization vector
according to the following criteria: i) $\mathbf{q}_{i}$ is non-zero only for
differentially expressed genes; ii) $\mathbf{q}_{i}$ increases with the number
of disease genes in the vicinity of $i$. In particular, for a differentially
expressed gene $i$, we define a score $\phi_{i}$ as the proportion of its
neighbors within shortest path distance $2$ that are disease proteins [XL06],
namely:

$\phi_{i}=\frac{|N(i)\cap S|}{N(i)}+\frac{|N_{2}(i)\cap S|}{N_{2}(i)}$ (16)

where we remind that $S$ denotes the seed set in the random walk with restart.
Finally, scores are normalized to obtain restart probabilities, so that:

$\mathbf{q}_{i}=\frac{\phi_{i}}{\sum_{i=0}^{N}\phi_{i}}$ (17)

####  B.3.2 Integrating Biological Information and Gene Expression

We discussed above two orthogonal approaches to the design of personalization
vectors. The first one leverages similarities between the biological processes
associated to known disease proteins and those to be prioritized. Hence,
teleporting probabilities depend on the seed set through association with
common biological processes.

In the second case, teleporting probabilities depend on information that is
tissue-specific (the level of expression in a population of subjects affected
by a certain disease) and partly on the seed set, but this time through the
PPI’s topology. Hence, these two approaches largely rely on complementary
sources of information. In the following paragraphs, we discuss ways to
integrate these complementary sources into a unique personalization vector
that leverages both.

Assume we have computed two personalization vectors $\mathbf{q}_{1}$ and
$\mathbf{q}_{2}$, the former using biological information only, the latter
using gene expression data and the PPI.

A first approach to integrating them is considering a convex combination,
namely:

$\mathbf{q}=\alpha\mathbf{q}_{1}+(1-\alpha)\mathbf{q}_{2},$ (18)

where $\alpha\in[0,1]$. Parameter $\alpha$ allows to weigh the importance of
the different information sources we are using. In our experiments, we set
$\alpha=1/2$. Intuitively, this type of aggregation amounts to considering a
gene a potential candidate if it is statistically significant in terms of its
involvement in biological processes, gene expression levels in the subject
group, or both.

Remark. Despite its simplicity, this is a mathematically principled choice. In
particular, it is well-known [JW03] and easy to show that the stationary
distribution corresponding to the convex combination of two personalization
vectors $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$ is itself the linear combination
of the stationary distributions corresponding to $\mathbf{q}_{1}$ and
$\mathbf{q}_{2}$ respectively. Briefly put, the parameter $\alpha$ allows us
to tune the relative importance of the information provided by gene
annotations and gene expression levels respectively. Moreover, this approach
extends seamlessly to an arbitrary number of heterogeneous sources of
information (and corresponding personalization vectors).

Alternatively, we can consider a gene to be a potential candidate if it is
statistically significant both in terms of the biological processes it is
involved in (with respect to the seed of known disease genes) and of its
expression levels in the population of subjects. In this case, a natural way
to define $\mathbf{q}$ is as follows:

$\mathbf{q}_{i}=\frac{\mathbf{q}_{1i}\cdot\mathbf{q}_{2i}}{\sum_{j=1}^{n}\mathbf{q}_{1j}\cdot\mathbf{q}_{2j}}.$
(19)

Note that, again, $\mathbf{q}$ is a probability distribution over the genes.

###  B.4 Biological Transition Matrix

####  B.4.1 Computing a Transition Matrix

Similar approaches can be used to derive a transition matrix for the random
walk with restart. Both approaches rely on the PPI, differing on the way PPI’s
edges are assigned weights that reflect genes’ similarity and determine the
probabilties of edge traversals.

####  B.4.2 Using Biological Information

In this case, we consider a weighted transition matrix $W$, in which each
entry $W_{ij}$ depends on the extent to which nodes genes $i$ and $j$ of the
PPI share common annotations (i.e., they are involved in common biological
processes) that are also significant for the disease. In more detail,
considered genes $i$ and $j$, we define the following _Disease Specific
Interaction Function_ (DSI function in the remainder):

$DSI(i,j)=|\mathcal{A}(i)\cap\mathcal{A}(j)\cap\mathcal{F}|,$ (20)

where we remind that $\mathcal{A}(k)$ denotes the set of gene $k$’s
annotations, while $\mathcal{F}$ denotes the overall set of annotations that
are statistically significant for disease genes.

Intuitively, the higher $DSI(i,j)$, the more $i$ and $j$ share annotations
that are also statistically significant for the disease under consideration.

In the end, edges in the PPI will be assigned weights depending on the DIS
function as follows:

$W_{ij}=\begin{cases}c+DSI(i,j)&\text{if }(i,j)\in E,\\\
0&otherwise.\end{cases}$ (21)

Here, $c$ is a positive constant that accounts for usual sparsity of the
available datasets, so that no biological information may be available for the
end-points of a link in the PPI. In this case, the link receives a minimum
weight $c$.

####  B.4.3 Using Co-Expression Information

In this case, we use gene expression information about a population of
patients, in order to augment the PPI with tissue-specific information.
Specifically, gene expression information is used to assign weights to edges,
this time reflecting similarities between genes in terms of co-expression with
respect to the subject population.

Consider a CO-Expression network in which each pair of gene $(i,j)$ is
associated a score defined by the Pearson’s correlation coefficient
$pc_{i,j}$. We need $pc_{i,j}$ to be a probability distribution on PPI network
for each gene $i\in V$. Thus, we define the probability of the random walker
to go from node i to one of its neighbors j in the PPI network as:

$W_{ij}=\frac{pc_{ij}}{\sum_{k\in N(i)}pc_{ik}}$ (22)

where N(i) is the set of i’s neighbors in the PPI network.

####  B.4.4 Computing an Aggregate Transition Matrix

In the previous sections, we have seen how we can derive random walk’s
transition matrices using only information about biological processes (i.e.,
annotations) or gene expression data from a population of subjects. While
similar approaches apply to other prospective information sources, we show
below how transition matrices derived from the aforementioned sources can be
integrated into an _aggregate transition matrix_ , leveraging both information
sources.

Assume we have computed transition matrices $W_{1}$ and $W_{2}$ using
biological and gene expression information respectively. We proceed in a way
that is similar to the approach used to derive an aggregate personalization
vector.

A first approach is considering a convex combination. Namely, for
$\beta\in[0,1]$, we shall consider the transition matrix

$W=\beta W_{1}+(1-\beta)W_{2}.$

Note that choosing $\beta=1$ corresponds to only considering biological
annotations, whereas $\beta=0$ corresponds to only leveraging gene expression
data. So, $\beta$ is a parameter, whose tuning allows us to weigh the
importance of one source of information with respect to the other. This
approach can be extended to the case of more than two information sources in
the obvious way.

###  B.5 Network Randomization

For the sake of completeness, in this section, we succintly describe the
degree-preserving, network randomization algorithm proposed in [MKI+03]. Given
an initial network $G(V,E)$, the algorithm involves carrying out a series of
Monte Carlo switching steps whereby a pair of edges $(A\rightarrow
B,C\rightarrow D)$ is selected at random and the ends are exchanged to give
$(A\rightarrow D,C\rightarrow B)$.

Supplementary Fig. 3: Randomization process: a) shows the original edges in
the network. b) shows the edges computed after the switching phase

Supplementary Fig. 3 shows the original edges and the switched edges after one
step of the Markov Chain Monte-Carlo Switching Algorithm. However, the
switching phase is not performed in two different cases:

  * •

If it generates multiple edges. For instance if it chooses the pair of edges
$(A\rightarrow D,C\rightarrow D)$ a duplicate edges is created as shown in
figure 4.

Supplementary Fig. 4: Multiple edge generation: a) shows the original edges in
the network. b) shows the duplicated edge computed after the switching phase

  * •

If it generates self-edges. If it randomly select the pair of edges
$(B\rightarrow A,A\rightarrow C)$, a self-loop is created as show in figure

Supplementary Fig. 5: Multiple edge generation: a) shows the original edges in
the network. b) shows the self-loop computed after the switching phase

The entire process is repeated exactly $Q\cdot|E|$ times, where $E$ is the
number of edges in the original network and $Q$ is chosen large enough to have
a good mixing.

##  Appendix C Results

###  C.1 Performance indices

To assess accuracy in disease gene prioritization, for each disease analyzed,
we used a 100-fold random validation approach, with each random split
consisting of a training set of seed genes (corresponding to 70% of the total)
that were provided as input to the algorithms and a test set, accounting for
the remaining 30% of the known disease genes . To evaluate the power of each
heuristic, we considered two standard prediction indices in Data Mining:

  * •

Recall@$k$: this is the fraction of the test set that is successfully
retrieved in the top $k$ positions of the ranking computed by the algorithm.

  * •

nDCG: The normalized Discounted Cumulative Gain (nDCG) is proven to be able to
select the better ranking between any two, substantially different rankings.
For binary classification the nDCG is given by:

$\frac{\sum_{i\in
P}\frac{1}{log_{2}(i+1)}}{\sum_{i=0}^{|P|}\frac{1}{log_{2}(i+1)}}$

Where summation in the numerator runs over all positive instances, while
summation in the denominator quantifies the ideal case, in which positive
instances appear in the top positions of the algorithm’s ranked list.

###  C.2 Multi-Omics Integration

Supplementary Fig. 6: Benchmark on each cancer disease analyzed to pick the
best combination of hyper parameters $(\alpha,\beta,r)$. Tab a) shows the
ccomparison with respect to Recall@k and Tab b) comparison with respect to
nDCG, for different values of the teleporting probability. Regardless of the
choice of $r$, results (in terms of Recall@k and nDCG) were best when
$\alpha=1$ and $\beta=0.5$. These values correspond to i) considering only
ontologies to the purpose of computing personalization vectors and ii) giving
equal importance to ontologies and gene expression when computing the
transition matrix of the random walk.

To identify the best combination of hyper-parameters, we tested how different
integration approaches affect BRW’s performances on the same validation
benchmark discussed in the main article, in terms of Recall@k and nDCG. In
more detail, we considered 3 different validation approaches: i) Monte Carlo
Random Sampling, ii) Drug Target Discovery and iii) Drug Enrichment.

##### Monte Carlo Sampling.

We considered all cancer phenotypes analyzed in the main article, and for each
phenotype we considered different values of $\alpha$ and $\beta$, fixing the
restart probability $r$ to 0.5 in order to reduce grid search space. We found
out that manually curated disease genes have a bias on the ontologies (GO
biological processes, KEGG and Reactome pathways). Indeed, the best
performance are obtained for large values of $\alpha$ as shown in
Supplementary Fig. 7. I.e., it is better to teleport to biologically similar
proteins than differentially expressed genes. As regards the aggregated
transition matrix, performance in terms of Recall@k/nDCG is maximized when
$\beta=0.5$. In other words, performance is maximized when co-expression
information is combined with biological information (and the PPI). Once we
discovered the best combination of $(\alpha,\beta)$, we tried different values
of the restart probability, with results summarized in Supplementary Fig. 6.
As Supplementary Fig. 6 shows, results are not greatly affected by the value
of the restart probability, as long this as this exceeds the value $0.25$.

##### Drug target discovery and drug enrichment.

In this case, we considered Breast Cancer phenotype and we analyzed drug
enrichment and performance in drug target discovery of each combination of
hyper-parameters $(\alpha,\beta,r)$, in terms of Recall@k and nDCG. Table 3
shows the result of the drug enrichment, Recall@k, nDCGs and number of FDA
approved, Breast Cancer drugs that were annotated, for each combination of
parameters. As shown in the table, different parameter combinations perform
best for different validation tasks. For example, parameter combinations
corresponding to a bias toward ontologies perform best in the task of
identifying Breast Cancer drug targets of approved drugs, highly scoring drug
targets in the first 200 positions, whereas biasing the random walk towards
differentially expressed genes and the co-expression network performs best in
prioritizing genes that are enriched in many relevant drugs. Indeed,
Paclitaxel, Fulvestrant, Doxorubicin,Tamoxifen, Methotrexate, and letrozole,
well-known approved Breast Cancer drugs, obtain are significantly enriched
when gene expression is weighted more than ontologies. Possible explanations
for this phenomenon are given in Sections 3.4.1 and 4.

Supplementary Fig. 7: Benchmark on each cancer disease analyzed to identify
the best combination of pair of hyper parameters $(\alpha,\beta)$. Each row
represents the following tumor type analyzed in this manuscript: BRCA, LUAD,
THCA, and COAD. Each row plots the comparison of different runs of the BRW for
same value of $\alpha$ and different values of $\beta$, while restart
probability is set to 0.5 in all experiments.

Hyper Parameter Combination ($\alpha$, $\beta$, $r$) |  (0.25, 0.25, 0.25) |  (0.25, 0.25, 0.5) |  (0.25, 0.25, 0.75) |  (0.25, 0.5, 0.25) |  (0.25, 0.5, 0.5) |  (0.25, 0.5, 0.75) |  (0.25, 0.75, 0.25) |  (0.25, 0.75, 0.5) |  (0.25, 0.75, 0.75) |  (0.5, 0.25, 0.25) |  (0.5, 0.25, 0.5) |  (0.5, 0.25, 0.75) |  (0.5, 0.5, 0.25) |  (0.5, 0.5, 0.5) |  (0.5, 0.5, 0.75) |  (0.5, 0.75, 0.25) |  (0.5, 0.75, 0.5) |  (0.5, 0.75, 0.75) |  (0.75, 0.25, 0.25) |  (0.75, 0.25, 0.5) |  (0.75, 0.25, 0.75) |  (0.75, 0.5, 0.25) |  (0.75, 0.5, 0.5) |  (0.75, 0.5, 0.75) |  (0.75, 0.75, 0.25) |  (0.75, 0.75, 0.5) |  (0.75, 0.75, 0.75)
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
estradiol |  18.7 |  10.8 |  10.1 |  11.9 |  7.4 |  8.8 |  10.1 |  7.4 |  8.1 |  11.8 |  7.9 |  5.9 |  12.2 |  7.0 |  6.2 |  10.9 |  7.1 |  6.2 |  8.5 |  4.2 |  3.3 |  11.5 |  4.2 |  3.7 |  11.6 |  5.2 |  3.7
imatinib |  15.5 |  6.9 |  5.9 |  11.0 |  6.4 |  5.8 |  8.3 |  7.1 |  5.5 |  13.5 |  7.9 |  5.9 |  10.6 |  7.4 |  6.9 |  8.6 |  7.4 |  7.1 |  10.5 |  7.8 |  7.1 |  10.4 |  8.0 |  6.1 |  10.3 |  7.3 |  6.4
arsenic trioxide |  9.5 |  6.3 |  5.1 |  11.0 |  5.8 |  5.0 |  10.3 |  7.1 |  5.3 |  15.6 |  9.8 |  5.5 |  15.0 |  10.8 |  6.9 |  15.3 |  12.3 |  6.5 |  15.6 |  10.8 |  11.8 |  17.7 |  13.9 |  10.8 |  17.1 |  13.9 |  11.3
bexarotene |  7.2 |  4.6 |  4.8 |  9.5 |  5.2 |  4.1 |  9.8 |  5.7 |  3.7 |  5.7 |  6.7 |  4.3 |  8.9 |  7.0 |  4.3 |  9.3 |  7.4 |  4.3 |  5.2 |  3.7 |  3.6 |  8.0 |  5.2 |  4.6 |  11.0 |  6.5 |  4.0
paclitaxel |  4.8 |  6.3 |  4.6 |  4.7 |  5.0 |  4.5 |  4.2 |  5.2 |  4.2 |  3.4 |  4.2 |  4.7 |  3.8 |  3.9 |  4.7 |  3.5 |  4.2 |  4.5 |  3.0 |  2.9 |  2.7 |  4.0 |  2.5 |  2.4 |  4.6 |  3.1 |  2.5
sirolimus |  10.1 |  4.5 |  4.3 |  11.0 |  4.6 |  4.2 |  10.1 |  4.7 |  4.6 |  8.9 |  7.0 |  4.3 |  13.4 |  6.6 |  3.8 |  15.3 |  7.1 |  3.9 |  12.2 |  7.8 |  5.0 |  17.7 |  11.3 |  4.9 |  19.1 |  13.4 |  5.3
fulvestrant |  9.1 |  4.4 |  4.2 |  5.7 |  2.6 |  4.1 |  4.0 |  2.9 |  3.8 |  5.7 |  3.1 |  2.0 |  4.7 |  2.2 |  2.1 |  2.9 |  2.0 |  2.3 |  3.0 |  0.9 |  1.1 |  2.5 |  1.1 |  0.9 |  2.4 |  1.1 |  1.0
valproic acid |  8.5 |  3.7 |  3.5 |  7.3 |  3.2 |  3.5 |  4.4 |  3.3 |  3.4 |  5.7 |  4.2 |  3.1 |  5.1 |  3.6 |  3.1 |  5.8 |  3.1 |  3.1 |  3.8 |  3.0 |  2.4 |  5.1 |  3.8 |  2.4 |  5.9 |  4.7 |  2.8
vitamin c |  9.1 |  4.4 |  3.5 |  5.5 |  3.4 |  3.2 |  6.6 |  3.3 |  2.9 |  7.0 |  4.5 |  2.0 |  6.6 |  3.6 |  2.0 |  6.1 |  3.3 |  2.1 |  6.9 |  3.3 |  2.4 |  6.6 |  4.7 |  2.5 |  8.8 |  5.3 |  2.5
doxorubicin |  7.6 |  4.1 |  3.5 |  5.7 |  3.6 |  3.6 |  6.6 |  4.4 |  3.8 |  5.7 |  5.8 |  4.5 |  6.6 |  4.9 |  4.1 |  6.5 |  5.3 |  3.9 |  5.1 |  5.6 |  4.8 |  7.6 |  6.6 |  4.6 |  7.6 |  6.7 |  4.6
tamoxifen |  12.4 |  5.3 |  3.2 |  5.6 |  3.3 |  2.6 |  6.3 |  4.0 |  2.5 |  6.2 |  4.9 |  2.4 |  5.7 |  4.3 |  2.9 |  7.4 |  4.0 |  3.2 |  5.5 |  3.5 |  3.7 |  5.7 |  3.5 |  3.8 |  6.4 |  4.0 |  4.3
methotrexate |  14.2 |  6.5 |  3.2 |  10.2 |  5.6 |  2.6 |  6.8 |  6.0 |  2.6 |  9.7 |  6.9 |  2.9 |  8.3 |  5.4 |  3.2 |  6.9 |  4.7 |  3.2 |  6.8 |  3.0 |  2.2 |  7.0 |  3.0 |  1.8 |  6.3 |  3.6 |  1.8
letrozole |  8.7 |  4.9 |  3.2 |  4.0 |  2.4 |  3.1 |  2.6 |  2.9 |  2.6 |  3.1 |  3.9 |  1.5 |  2.7 |  3.0 |  1.9 |  1.8 |  3.0 |  2.3 |  3.0 |  1.6 |  2.0 |  3.1 |  1.8 |  2.0 |  2.3 |  2.3 |  2.0
etoposide |  8.9 |  2.9 |  2.8 |  5.1 |  2.6 |  2.8 |  5.7 |  3.0 |  3.1 |  6.7 |  3.6 |  2.0 |  6.2 |  2.7 |  2.0 |  6.3 |  2.7 |  2.1 |  6.1 |  3.9 |  2.2 |  7.5 |  4.4 |  1.9 |  6.4 |  4.0 |  1.9
deferoxamine |  4.5 |  2.8 |  2.7 |  3.7 |  4.1 |  2.2 |  4.4 |  4.2 |  2.2 |  5.9 |  3.6 |  2.1 |  3.7 |  4.2 |  2.1 |  4.5 |  4.2 |  2.7 |  4.3 |  4.2 |  2.2 |  3.8 |  4.2 |  2.8 |  5.4 |  5.1 |  2.8
vorinostat |  7.7 |  4.0 |  2.6 |  7.1 |  3.6 |  2.6 |  8.4 |  4.7 |  2.5 |  7.5 |  5.2 |  4.3 |  7.6 |  4.7 |  4.1 |  9.0 |  5.2 |  4.0 |  5.2 |  6.0 |  4.2 |  7.6 |  6.9 |  4.6 |  9.1 |  7.3 |  5.0
cisplatin |  7.4 |  3.9 |  2.4 |  4.6 |  3.4 |  2.5 |  4.7 |  3.5 |  2.4 |  3.9 |  3.1 |  1.7 |  3.9 |  3.0 |  1.7 |  4.1 |  3.1 |  1.7 |  4.3 |  3.7 |  2.6 |  4.0 |  3.1 |  2.6 |  6.0 |  4.0 |  2.9
cetuximab |  2.8 |  2.8 |  2.3 |  4.5 |  3.5 |  2.8 |  4.3 |  3.8 |  2.8 |  5.1 |  4.6 |  4.0 |  5.1 |  4.9 |  3.7 |  4.9 |  5.5 |  3.6 |  5.3 |  7.2 |  6.2 |  5.8 |  6.6 |  6.2 |  5.6 |  6.3 |  6.4
zinc acetate |  4.7 |  1.7 |  2.0 |  4.7 |  2.1 |  2.0 |  7.1 |  2.2 |  2.0 |  4.1 |  3.6 |  1.6 |  7.8 |  3.5 |  1.6 |  7.2 |  4.0 |  1.6 |  3.6 |  2.1 |  2.1 |  7.9 |  2.6 |  2.1 |  6.1 |  4.1 |  2.1
mifepristone |  4.6 |  3.3 |  1.7 |  5.0 |  3.7 |  1.7 |  5.0 |  3.4 |  1.7 |  3.9 |  3.8 |  2.0 |  5.0 |  4.2 |  1.7 |  4.5 |  4.3 |  1.7 |  3.9 |  3.0 |  1.7 |  4.6 |  2.9 |  2.1 |  6.9 |  2.6 |  2.1
trovafloxacin |  5.1 |  2.3 |  1.6 |  3.4 |  2.5 |  1.4 |  4.3 |  3.3 |  1.1 |  5.7 |  4.0 |  2.0 |  5.1 |  3.6 |  2.0 |  4.9 |  4.1 |  2.0 |  5.0 |  4.2 |  2.5 |  5.5 |  4.3 |  2.5 |  6.2 |  3.9 |  2.5
plicamycin |  5.8 |  4.1 |  1.6 |  5.7 |  3.9 |  1.8 |  5.5 |  4.4 |  1.8 |  5.4 |  4.7 |  2.9 |  4.9 |  6.6 |  3.2 |  4.2 |  6.4 |  3.4 |  4.5 |  5.5 |  5.5 |  4.2 |  5.3 |  5.7 |  4.8 |  5.3 |  5.7
hydrocortisone |  2.7 |  1.6 |  1.3 |  4.8 |  2.1 |  1.8 |  5.7 |  2.5 |  2.1 |  4.0 |  3.3 |  2.2 |  5.7 |  4.2 |  2.6 |  6.2 |  4.2 |  2.6 |  5.5 |  5.5 |  3.6 |  5.8 |  5.1 |  4.0 |  7.3 |  4.7 |  4.1
vincristine |  2.0 |  1.1 |  0.7 |  1.6 |  1.5 |  0.6 |  2.4 |  1.5 |  0.6 |  2.3 |  1.5 |  0.5 |  3.2 |  2.3 |  0.4 |  4.4 |  2.3 |  0.6 |  2.3 |  1.5 |  0.6 |  4.4 |  1.8 |  0.6 |  6.3 |  2.3 |  0.8
hydroquinone |  3.7 |  2.3 |  0.6 |  3.3 |  2.7 |  0.8 |  4.5 |  3.1 |  0.8 |  4.0 |  3.2 |  2.2 |  3.7 |  3.9 |  2.2 |  4.5 |  3.5 |  2.6 |  3.6 |  3.5 |  2.3 |  4.1 |  3.5 |  2.4 |  5.6 |  4.0 |  2.4
isotretinoin |  1.5 |  0.9 |  0.4 |  1.5 |  1.4 |  0.4 |  2.1 |  1.4 |  0.5 |  2.0 |  1.4 |  2.3 |  2.7 |  1.5 |  2.3 |  2.8 |  2.0 |  2.3 |  4.1 |  4.7 |  5.3 |  3.5 |  4.0 |  5.3 |  2.8 |  4.1 |  5.3
Recall@200 |  0.16 |  0.09 |  0.07 |  0.19 |  0.11 |  0.09 |  0.19 |  0.12 |  0.09 |  0.19 |  0.16 |  0.12 |  0.19 |  0.19 |  0.14 |  0.23 |  0.19 |  0.14 |  0.21 |  0.25 |  0.21 |  0.23 |  0.26 |  0.23 |  0.25 |  0.26 |  0.23
nDCG |  0.35 |  0.33 |  0.32 |  0.37 |  0.34 |  0.32 |  0.40 |  0.34 |  0.32 |  0.39 |  0.36 |  0.34 |  0.40 |  0.37 |  0.34 |  0.41 |  0.37 |  0.34 |  0.39 |  0.40 |  0.37 |  0.40 |  0.40 |  0.37 |  0.41 |  0.41 |  0.38
N. of Breast Cancer Annotated Drugs |  15 |  9 |  6 |  17 |  9 |  7 |  17 |  12 |  7 |  17 |  15 |  6 |  17 |  17 |  9 |  17 |  17 |  9 |  17 |  17 |  13 |  17 |  17 |  15 |  17 |  17 |  15

Table 3: This table describes drug enrichment and drug target discovery
results for each combination of hyper-parameters $(\alpha,\beta,r)$. The table
is divided into four sections: The first one (the header) represents all
hyper-parameter combinations tested. Each combination is described by the
corresponding triple $(\alpha,\beta,r)$. The second section shows the
$-log_{10}P_{value}$ of drugs enriched in at least one of the combinations. We
considered only drugs with $P_{value}<10^{-5}$ in at least one of the
combinations. Drugs in bold are FDA-approved to treat breast cancer. The third
section shows the portion of Breast cancer’s drug targets predicted by each
combination in the top 200 positions (Recall@200) and the nDCG. Finally, the
last section shows the number of Breast cancer’s FDA-approved drugs with at
least a drug target predicted in the first 200 positions by each combination
of Hyper-Parameters.

###  C.3 Robustness and Bias

Supplementary Fig. 8: BRW and RWR-M comparison. Tab a) shows the comparison
between BRW and RWR-M on ranking test known disease genes when at least one
network is randomized. As we can see, enriched ontologies (GO biological
Processes, Pathways) correlate positively between seed and validation set.
Consequently, the noise induced by the randomized PPI is mitigated.
Unfortunately, even if RWR-M relies on several biological information (PPI-
Network, CO-Expression Network and Pathways Network) it is not able to
mitigate the noise as much is mitigated by Biological Random Walks. Tab b)
illustrates how multi-omics integration affect the bias induced by curated
ontologies (GO, Reactome) and the PPI network

.

In a first experiment, we performed a mean 100-fold validation for each cancer
phenotype by randomly sampling 70% of known disease genes and using the rest
to compute the average Recall@k of each tested algorithm. To begin, we tried
different values of the restart probability to understand how the biologically
weighted teleporting probability is correlated with the test genes. We found
out that this correlation is strong, and thus when a very high restart
probability is used, BRW improves its prediction performances, as shown in
Supplementary Fig. 8 a).

In a second round of experiments, we tried to understand what biological
information used by RWR-M was critical for test gene retrieval. To this
purpose, we started by only randomizing only the PPI network and keeping the
original pathways network used by RWR-M. In a second step, we did the
opposite. We discovered that the bias induced by Graphite is more critical in
predicting test genes than using the original PPI network. Indeed,
Supplementary Fig. 8 a) shows that the average Recall of the first 200
positions computed on the randomized PPI network is always higher than the
average Recall@k compute when we randomize the (pathways) Graphite network.
Finally, we executed RWR-M using randomized versions of both networks as
inputs to the algorithm, and its prediction performance dropped drastically.

Furthermore, to understand if candidate genes predicted using randomized PPI
networks are more biologically meaningful than those returned using the
original PPI, we considered the test suite proposed by [LBLB21]. We evaluated
RWR-M and BRW on the Lung Cancer phenotypic condition using the BIOGRID PPI
network, both implemented in the test suite. Since RWR based algorithms depend
on restart probability (i.e., greater values of $r$ make the algorithms depend
less on the topology of the network and, in the extreme case when $r=1$, the
Markov chain converges to the personalization vector making the topology of
the network completely useless), we set $r=0.25$. Moreover, we chose
$\alpha=1$ and $\beta=1$ when running BRW; in this way, we isolated the bias
induced by manually curated ontologies and not by gene expression. We then
validated each predicted gene set using gseapy KEGG Pathways enrichment
analysis implemented in the test suite and we compared the outcomes returned
using the original PPI with those obtained using its randomized counterparts:
uniform, scale-free, epected degree, and shuffled.

Supplementary Fig. 8 b) shows the comparison between candidate genes predicted
using a randomized network and the original one. The x axis represents all
randomized methods that are available in the test suite, while the y axis
represents the negative log scale of p-values calculated by the test suite
using KEGG enrichment analysis provided by gseapy on predefined Lung Cancer
pathways.

###  C.4 Comparison with [GMP+19]

BRW is the evolution of a framework we proposed in [GMP+19] in its embryonic
form. Both rely on the hypothesis that genes involved in the same disease
should be more likely to share biological similarities. [GMP+19] is different
from the current work in significant ways. Our initial approach integrated
only Gene Ontology to define the personalization vector and the transition
matrix of the random walk, something that was done in an ad-hoc way. In
contrast, the framework proposed in this paper proposes a principled approach,
described in Section 2 of the main paper and Sections B.3.2 and B.4.4 of this
supplementary material, whereby multiple, heterogeneous biological datasets
such as KEGG, Reactome, and GO, can be seamlessly integrated to define a
biased random walk with restart. In particular, this new version effectively
integrates differential expression and co-expression data to improve the
biological meaning of candidate genes prioritized in the top K positions. To
demonstrate the power of incorporating multi-omics data (i.e., Proteome and
Transcriptome) and several manually curated ontologies (i.e., KEGG, Reactome,
and GO), in this section, we compare our framework with the algorithm proposed
in [GMP+19]. In the paragraphs that follow, the latter is referred to as BRW
CIBCB 2019.

Supplementary Fig. 9: Old vs New Approach comparison: a) and b) compare
respectively the Recall@K and nDCG of the previous version of BRW with the one
presented in this manuscript. c) compares BRW CIBCB 2019 with BRW on the
percentage of drug target prioritized in the top 200 positions. Finally, d)
compares the drug enrichment computed using the first 200 candidate genes
prioritized by the frameworks.

We compared our approach with [GMP+19] along two axes, similarly to what we
did for other baselines: i) we performed a 100-fold Montecarlo random sample
(internal) validation, dividing the set of known disease genes into a seed set
and a test set, which we used to compute Recall@K and nDCG; ii) we further
performed an external validation, to this purpose focusing on Breast Cancer
and comparing the two competing approaches in terms of the set of
statistically significant drugs they retrieve if we restrict to candidate
genes ranked the first 200 positions.

Supplementary Figures 9 a) and b) compare our approach with BRW CIBCB 2019 in
terms of average Recall@K and nDCG computed over the aforemention 100-fold
Montecarlo runs. BRW’s ability to leverage manually curated information from
KEGG and Reactome and gene expression data (i.e., differential expression and
co-expression network) drastically improves the algorithm’s performance in
prioritizing known candidate genes from the test set.

In the case of drug targets, the fraction of drug targets discovered by BRW
consistently increases with the number $k$ of top ranked genes we consider,
showing a two-fold improvement over [GMP+19] (Figure 9 c)). This behaviour is
further confirmed if one considers Breast Cancer drugs that are enriched by
genes in top 200 positions. In particular, Figure 9 d) shows that BRW returns
a list of candidates enriched in doxorubicin, methotrexate and tamoxifen, all
well-known drugs to treat Breast Cancer (https://www.cancer.gov/about-
cancer/treatment/drugs/breast).

Table 4: Drugs approved to treat Breast Cancer which have at least one drug
target predicted by at least one framework. If an algorithm prioritize a
target of a specific drug, its cell is highlighted and it contains 1. The
combination of Hyper-parameters used by BRW is $\alpha=0.75$, $\beta=0.75$,
$r=0.25$

Drug Bank ID | Drug Name | Associated Condition | BRW | DADA | DIAMOnD | RWR-M | RWR
---|---|---|---|---|---|---|---
DB12001 | Abemaciclib | Metastatic HR + HER2 - breast cancer | 1 | 0 | 0 | 0 | 0
DB00675 | Tamoxifen | Breast Cancer | 1 | 1 | 1 | 1 | 1
DB11828 | Neratinib | Breast Cancer | 1 | 1 | 0 | 1 | 1
DB01259 | Lapatinib | Metastatic Breast Cancer | 1 | 1 | 0 | 1 | 1
DB09073 | Palbociclib | Metastatic Breast Cancer | 1 | 0 | 0 | 0 | 0
DB11730 | Ribociclib | Metastatic Breast Cancer | 1 | 0 | 0 | 0 | 0
DB00072 | Trastuzumab | Breast Cancer | 1 | 1 | 0 | 1 | 1
DB00570 | Vinblastine | Refractory Breast cancer | 1 | 1 | 1 | 1 | 0
DB11760 | Talazoparib | Metastatic Breast Cancer | 1 | 0 | 1 | 0 | 0
DB09074 | Olaparib | Metastatic Breast Cancer | 1 | 0 | 1 | 0 | 0
DB00624 | Testosterone | Inoperable, metastatic Breast cancer | 1 | 1 | 1 | 1 | 1
DB06710 | Methyltestosterone | Metastatic Breast Cancer | 1 | 1 | 1 | 1 | 1
DB01248 | Docetaxel | Metastatic Breast Cancer | 1 | 1 | 0 | 0 | 1
DB01248 | Docetaxel | Locally Advanced Breast Cancer (LABC) | 1 | 1 | 0 | 0 | 1
DB08871 | Eribulin | Refractory, metastatic Breast cancer | 1 | 1 | 0 | 0 | 1
DB00286 | Conjugated estrogens | Metastatic Breast Cancer | 1 | 1 | 1 | 0 | 1
DB00481 | Raloxifene | Invasive Breast Cancer | 1 | 1 | 1 | 0 | 1
DB00445 | Epirubicin | Breast Cancer | 0 | 0 | 1 | 0 | 0
DB01204 | Mitoxantrone | Metastatic Breast Cancer | 0 | 0 | 1 | 0 | 0
DB00997 | Doxorubicin | Metastatic Breast Cancer | 0 | 0 | 1 | 0 | 0