The use of gene interaction networks to improve the identification of cancer driver genes

DriverNet

DriverNet uses a greedy algorithm to identify driver genes from a bipartite graph combining mutation frequency and differential expression. It utilizes the protein functional interaction network constructed by Wu, Feng & Stein (2010), which was constructed from various sources of information such as curated pathways with non-curated data including protein–protein interactions, gene co-expression, protein domain interaction, gene ontology (GO) annotations, and text-mined protein interactions. These provide various small molecules, proteins, complexes, post-translationally modified proteins, and nucleic acid sequences which are mapped onto the genome using online repositories such as UniProt (2010) and Entrez Genes (Maglott et al., 2011). It determines which interactions form part of the network by using a Bayes classifier, and eliminates those that do not fit. DriverNet predicts driver genes by considering the effect of mutated genes on the gene expression levels of interacting partners.

Using threshold cut-off values genes are categorized as expressed or not. In Fig. 2, blue nodes partition of the bipartite graph represent mutated genes whilst nodes in red represent their expression status for different patients. Genes in red are significantly expressed. The gene interaction network connects nodes between the two sets. In the identification of candidate driver genes, the guiding principle is to select as many red nodes as possible using the fewest number of blue nodes. At each stage of the greedy algorithm, mutated genes with the highest number of significant connections (such as G4 in Fig. 2) are selected as candidate driver genes.

VarWalker

VarWalker uses a Random Walk with Restart (RWR) algorithm. The network it uses was constructed using the Human Protein Reference Database (Keshava Prasad et al., 2009), a manually curated resource. It includes protein–protein interactions, catalytic reactions, and protein translocation events that have been evaluated against other repositories of human protein–protein interaction data in the public domain (Mathivanan et al., 2006). This network shows the Cancer Genome Census (CGC) genes (Futreal et al., 2004) tend to be located more closely to each other than other genes. Specifically, 71% of CGC genes are directly connected and 26% have a shortest path of two. VarWalker uses this trait to nominate candidate driver genes by ascertaining consensus across multiple samples for mutated genes that converge. An initial gene filtering process removes long genes that are more frequently mutated due to size.

The RWR method calculates a vector V that represents the proximity between a given node and all other nodes in the network by solving Eq. (1) in Fig. 3. The gene network is represented by a matrix A (see Fig. 3), if a gene i links to a gene j, i.e. i–j, then Aij<>0. In this case, its value would be the probability of moving to node j from i. If gene i does not link to gene j, then Aij = 0. In Fig. 1, from G4, it is possible to move directly to one of the three other nodes. The probability of moving to a directly connected node is proportional to the number of outgoing nodes from G4, in this case 1/3. The RWR is applicable as a proximity metric because after a sufficiently long time interval, the probability of being at G4 at a random time provides a measure of the proximity between G4 and all the other nodes. Figure 3 is the matrix representation of this equation for the 12-gene network in Fig. 1. (1)${\displaystyle \mathrm{V}=\left(1-r\right)\mathrm{AV}+r\mathrm{P}.}$

In this example a restart probability value of 0.1 is used for r. The 12 by 12 matrix A presented in Fig. 3 is derived from the 12-node network in Fig. 1. P is the vector in which the ith element holds the probability that the walker is at node i at time 0. In this case we start at G4, so the fourth element of P is 1, and all others are zero (see Fig. 3). The value of V is then calculated, to satisfy Eq. (1). The solution of this equation is vector V given in Fig. 3: ${\displaystyle \mathrm{V}=\left(0.13,0.10,0.13,0.22,0.13,0.05,0.05,0.08,0.04,0.03,0.04,0.02\right).}$

This solution indicates nearby nodes (G1, G3, and G5) with higher scores of 0.13. We can also determine G9 and G11 are equally distant from G4. With a large network, this can be computationally intensive, thus, this matrix solution can be replaced by an iterative solution.

DawnRank

DawnRank selects potential driver genes based on their impact on the overall differential expression of its downstream genes in the interaction network. In this network, all redundant edges are collapsed to single edges when aggregating networks from different databases. With this method an individual patient sample is used rather than a large cohort, so drivers are identified on a personalized level. This single patient approach is totally independent of the mutation frequency, and can therefore be considered focused on finding more infrequent or rare drivers. It classifies rare and even patient-specific mutations. This is the use of the ‘long tail phenomenon’ when selecting driver genes, which considers cancer mutations as being characterized by a small number of frequently mutated genes and a large number of infrequently mutated genes. Selected genes are compared to CGC and Pan Cancer standard driver gene list (Cancer Genome Atlas Research et al., 2013; Tamborero et al., 2013) for validation.

DawnRank uses the PageRank family of algorithms to rank genes based on incoming links. PageRank (Page et al., 1999) was developed to measure the human interest in web pages. It is used by Google’s search engine to rank web pages. To illustrate this, consider a sub-network from our 12-node network in Fig. 1 as shown in Fig. 4A. Each of the 4 genes G1, G3, G4, and G5 is initially given the same rank. The initial rank for each of the 4 genes is therefore 0.25. Fig. 4A illustrates the topology of this simplified network at initial state t = 0. The next iteration involves updating the rank of each gene by adding up the rank of each incoming gene divided by the number of outgoing links from it. This is illustrated in Fig. 4B. Thus, new ranks shown in Fig 5 are calculated as follows: ${\displaystyle \begin{array}{cccc}\mathrm{G1}=\mathrm{G3}∕2=0.25∕2=0.125;\mathrm{G3}=\mathrm{G1}∕2=0.25∕2=0.125;\hfill & \hfill & \hfill & \hfill \\ \mathrm{G4}=\mathrm{G1}∕2+\mathrm{G3}∕2+\mathrm{G3}∕1=0.25∕2+0.25∕2+0.25∕1=0.375;\mathrm{G5}=0.\hfill & \hfill & \hfill & \hfill \end{array}}$

The ranks of genes may be weighted so that a gene is given a higher rank, even though there are fewer links to it, if more important genes link to it. The output of the PageRank algorithm is a list of genes and their rankings based on the gene network configuration. A high PageRank score for a mutated gene in cancer indicates that the gene is more likely to be a driver.

For DawnRank’s implementation of the PageRank algorithm the initial rank value for each gene would be 1/11,648, as the network of genes consists of 11,648 members. A gene linked to many other genes with high ranks receives a high rank. This process is modeled using states, the transitions from one state to another depending only on the current state rather than a preceding state. This is the Markov property, where each iteration is equally probable. The difference in the ranks between time t = 0 and t = 1 is computed recursively as r_t+1 and r_t, until it converges to an insignificant value (epsilon). It can also stop after a set number of iterations, which is 100 for DawnRank.

Interaction network construction

Network data files from the DriverNet, DawnRank and VarWalker packages were used. These were transformed into individual graphs and combined using the igraph library (Csardi & Nepusz, 2006) as shown in the dataflow diagram in Fig. 6.

For each of the three individual networks and the combined network we computed an interaction score assuming independence. To infer these interaction scores, we combined for each interaction between gene G_i and gene G_j, two scores:

• Q_ij is the number of graphs the interaction between gene G_i and gene G_j occurs in, represented on a common scale [0, 1]. This is where q_ij = 0 represents no information about the interaction, and q_ij = 1 represents strong evidence for the interaction as it occurs in all the graphs;

• R_ij is a count of the number of v-structures the edge G_i->G_j is part of in the network projected onto a common scale [0, 1] (see Fig. 7).

Q_ij and R_ij were represented as matrices, with genes identifying both rows and columns, q_ij and r_ij are the scores for the interaction between gene G_i and gene G_j. These two were combined as S_ij as in Eq. (2). This schema for combining scores allows us to adjust w between [0, 1] depending on the confidence in each of these contributors to the final weighting. For our implementation we used equal weighting by setting w to 0.5 (2)${\displaystyle {\mathrm{S}}_{\mathrm{ij}}={w\mathrm{Q}}_{\mathrm{ij}}+\left(1-w\right){\mathrm{R}}_{\mathrm{ij}}.}$

Assessing the linear bias correction

In order to quantify the impact of the dependency on the interaction scores, we compared the sum of the interactions from the individual graphs to those produced from the combined network for those interactions common in all three networks. The individual scores were summed as in Eq. (3) across the 3 graphs, and projected onto a common scale [0, 1]. (3)${\displaystyle \mathrm{Sum}=\sum {\mathrm{S}}_{\mathrm{i}}.}$

Eq. (4) was used to calculate the bias. We applied a linear regression between the summed values and the calculated values using Eq. (2) for the combined network. (4)${\displaystyle \text{Combined network}=\alpha \ast \mathrm{sum}+\beta .}$

The parameters α and β represent the bias for the interaction weights of the combined network.

Testing the weighted network

Interactions with a weight less than or equal to the cut-off value of 0.17 were discarded. The resultant network was used to analyse the prediction of driver genes by DriverNet and DawnRank.

Data consisting of 504 samples of breast cancer (BRCA) initially from TCGA were derived from DawnRank. These included somatic mutation and differential gene expression data between the cancer and normal transcriptome. Driver genes were predicted by DawnRank using its standard network and the weighted combined network. Mutation and expression datasets consisting of 178 cervical cancer samples were downloaded from the Catalogue of Somatic Mutation in Cancer (COSMIC) (Forbes et al., 2011). These were transformed into two binary matrices where the rows were patients and the columns were genes. For the purpose of our analysis, expression values between the range −2 and 2 were considered to be normal. Thus in the expression matrix, if a z-score value was >2.0 or <−2 the binary matrix element was set to 1 (TRUE), otherwise it was set to 0 (FALSE). Glioblastoma Multiforme (GBM) samples from The Cancer Genome Atlas (TCGA) (Cancer Genome Atlas Research, 2008) were used from DriverNet. These were represented by 2 matrices with 200 rows and 1,255 columns. Driver genes were predicted using DriverNet for GBM and cervical cancer.

The analysis included sensitivity, specificity, accuracy and receiver operating characteristic (ROC) with Area under the curve (AUC) measures (Zhu, Zeng & Wang, 2010).

Genes predicted as candidate driver genes were classified as true positives if they were present in the 125 driver genes from Vogelstein et al. (2013). Details on this analysis can be found in Supplementary Files.