Network-based coverage of mutational profiles reveals cancer genes

TL;DR

This paper presents nCOP, a network-based method that identifies cancer-related genes by analyzing individual mutational profiles within protein interaction networks, effectively uncovering both common and rare cancer genes across diverse tumor samples.

Contribution

The study introduces a novel optimization approach using integer linear programming to detect small gene subnetworks covering many individuals, improving cancer gene discovery beyond existing methods.

Findings

nCOP outperforms existing methods in identifying cancer genes.

The approach uncovers both frequently and infrequently mutated cancer genes.

Large-scale evaluation across 24 cancer types validates the method's effectiveness.

Abstract

A central goal in cancer genomics is to identify the somatic alterations that underpin tumor initiation and progression. This task is challenging as the mutational profiles of cancer genomes exhibit vast heterogeneity, with many alterations observed within each individual, few shared somatically mutated genes across individuals, and important roles in cancer for both frequently and infrequently mutated genes. While commonly mutated cancer genes are readily identifiable, those that are rarely mutated across samples are difficult to distinguish from the large numbers of other infrequently mutated genes. Here, we introduce a method that considers per-individual mutational profiles within the context of protein-protein interaction networks in order to identify small connected subnetworks of genes that, while not individually frequently mutated, comprise pathways that are perturbed across (i.e., "cover") a large fraction of the individuals. We devise a simple yet intuitive objective function that balances identifying a small subset of genes with covering a large fraction of individuals. We show how to solve this problem optimally using integer linear programming and also give a fast heuristic algorithm that works well in practice. We perform a large-scale evaluation of our resulting method, nCOP, on 6,038 TCGA tumor samples across 24 different cancer types. We demonstrate that our approach nCOP is more effective in identifying cancer genes than both methods that do not utilize any network information as well as state-of-the-art network-based methods that aggregate mutational information across individuals. Overall, our work demonstrates the power of combining per-individual mutational information with interaction networks in order to uncover genes functionally relevant in cancers, and in particular those genes that are less frequently mutated.

Figures

b). We observe that the coverage on the validation sets is much poorer. Though coverage of course increases as more nodes are added, it never exceeds $50\%$ even when $\alpha$ is increased to $1$ or when we have nearly perfect coverage on the training set, despite adding many more nodes. This poor performance is consistent with the expectation that synonymous mutations do not result in altered protein sequences and do not disturb cellular pathways. Hence, given the differences observed between using missense versus silent mutation data when varying settings for $\alpha$ and comparing training

b). We observe that the coverage on the validation sets is much poorer. Though coverage of course increases as more nodes are added, it never exceeds $50\%$ even when $\alpha$ is increased to $1$ or when we have nearly perfect coverage on the training set, despite adding many more nodes. This poor performance is consistent with the expectation that synonymous mutations do not result in altered protein sequences and do not disturb cellular pathways. Hence, given the differences observed between using missense versus silent mutation data when varying settings for $\alpha$ and comparing training

b). We observe that the coverage on the validation sets is much poorer. Though coverage of course increases as more nodes are added, it never exceeds $50\%$ even when $\alpha$ is increased to $1$ or when we have nearly perfect coverage on the training set, despite adding many more nodes. This poor performance is consistent with the expectation that synonymous mutations do not result in altered protein sequences and do not disturb cellular pathways. Hence, given the differences observed between using missense versus silent mutation data when varying settings for $\alpha$ and comparing training

b). We observe that the coverage on the validation sets is much poorer. Though coverage of course increases as more nodes are added, it never exceeds $50\%$ even when $\alpha$ is increased to $1$ or when we have nearly perfect coverage on the training set, despite adding many more nodes. This poor performance is consistent with the expectation that synonymous mutations do not result in altered protein sequences and do not disturb cellular pathways. Hence, given the differences observed between using missense versus silent mutation data when varying settings for $\alpha$ and comparing training

Equations

$$\alpha X+(1-\alpha )Size(G^{\prime }),$$

$$p_i$$

$$\sum _{i:(i,j)\in E}{y}_{i,j}$$

$$\sum _{j:(i,j)\in E}{y}_{i,j}$$

$$\sum _i{x}_i$$

$$p_i,x_i,y_{i,j}$$

Full text

Software download: http://www.cs.princeton.edu/$\sim$mona/software/ncop.html

Introduction

Large-scale cancer genome sequencing consortia, such as The Cancer Genome Atlas (TCGA) [35], the International Cancer Genome Consortium (ICGC) [37] and other smaller, cancer-specific studies have sequenced the protein-coding regions of thousands of tumor samples across tens of different cancer types. Initial analyses of these data have revealed that while there may be numerous somatic mutations in a tumor that result in altered protein sequences, very few are likely to play a role in cancer development  [2, 41, 14]. Therefore, a major challenge in cancer genomics is to develop methods that can distinguish the so-called “driver” mutations important for cancer initiation and progression from numerous other “passenger” mutations.

Early statistical approaches have identified cancer-driving genes by highlighting those genes that are mutated more frequently in a cohort of patients than expected by chance according to some background model [44, 10, 25]. However, the genetic underpinnings of cancer are highly heterogeneous: even when considering a single cancer type, very few genes are found to be somatically mutated across large numbers of individuals [18]. Further, genes altered only in a few individuals may also be important for tumorigenesis and cancer progression [34]. Clearly, these rarely mutated but cancer-relevant genes cannot be detected by purely frequency-based approaches.

A promising alternative viewpoint is to consider somatic mutations in the context of pathways instead of genes. In particular, it has been proposed that alterations within any of several genes comprising the same pathway can have similar consequences with respect to cancer development, and that this contributes to the mutational heterogeneity evident across cancers. Consistent with this, numerous analyses of TCGA data have shown that certain known pathways are frequently altered across tumor samples of a particular cancer via mutations in different genes [36, 28]. Early studies have leveraged this observation by analyzing known pathways for enrichment of somatic mutations [21, 4] and pinpointing those that are significantly mutated across patients [42, 40]. The power of these studies is somewhat limited, however, as our knowledge of pathways is incomplete and new pathways cannot be identified by these approaches.

De novo discovery of cancer-relevant pathways using large-scale protein interaction networks has thus been the focus of several newer methods (e.g., [39, 4, 7, 29, 5]). In particular, since protein-protein interaction networks have a modular organization [15, 32], proteins taking part in the same pathways and processes tend to be close to each other in the network. One prominent class of techniques leverages this modular structure by propagating mutational information through protein interaction networks and deriving pathways from the induced subnetworks [39, 26, 20, 1]. For instance, Vandin et al. [39] diffuse a “heat” signal arising from the frequency with which proteins are somatically mutated across a cohort of samples to uncover cancer-relevant modules while Hofree et al. [17] approach the problem from a different angle, using biological network information to stratify cancer subtypes. A recent pan-cancer network analysis [26] affirms the power of diffusing mutational data across protein interaction networks, especially for uncovering rarely mutated cancer genes. However, such diffusion approaches can be highly influenced by frequently mutated genes [26], and further, these methods do not consider whether most patients have mutations in any of the identified pathways.

Here we present a novel network-based approach to tackle cancer mutational heterogeneity by utilizing per-individual mutational profiles. Our method is based on the expectation that if a pathway is relevant for cancer, then (1) many individuals will have a somatic mutation within one of the genes comprising the pathway and (2) the genes comprising the pathway will interact with each other and together form a small connected subcomponent within the larger network. Therefore, given a biological network as well as patient sample data consisting of somatic point mutations, the goal of our approach is to find a set of candidate genes that both “cover” the most patients (i.e., individuals have mutations in one or more of these genes) and are connected in the network (i.e., these genes are likely to participate in the same cellular pathway or process). In contrast to network diffusion approaches, our framework focuses on per-individual mutational profiles and as a result, the “influence” of frequently mutated genes is not spread through the network. We note that network-based coverage approaches have been previously introduced to uncover pathways that are dysregulated [38, 6, 23] or mutated [9, 24] across cohort of samples. However, either patients were required to be covered by these approaches [38, 6, 23, 24], in some cases multiple times (which is especially relevant for dysregulated genes, since there many of them), or these approaches were designed for data sets with significantly fewer mutations [9]; both cases lead to very different optimizations and algorithms that are not effective for the task at hand.

We devise a simple yet intuitive objective function that balances identifying a small subset of genes with covering a large fraction of individuals. Our objective has just a single parameter that is automatically set using a series of cross-validation tests, eliminating the need of many previous approaches to manually select values for various thresholds and parameters. We develop an integer linear programming formulation to solve this problem and also give a fast heuristic algorithm. We apply our method—network-based coverage of patients (nCOP)—to 24 cancer types from TCGA and uncover both well-known cancer driver genes as well as new potential cancer-related genes. We compare nCOP to previous methods that do not use network information, including a state-of-the-art frequency-based method [25] and a “set cover” version of our approach that attempts to find a set of genes that covers cancer samples without considering network connectivity, and demonstrate nCOP’s superior power in detecting known cancer genes and in zooming in on rarely mutated ones. Finally, we compare nCOP to a recent network-based method that aggregates mutational information [5] and show that our per-patient approach readily outperforms it.

Methods

Overview. In this section, we give an overview of our methodology (see also Supplementary Figure 1); each part is described in more detail in the subsequent sections.

The biological network is modeled as an undirected graph $G$ where each vertex represents a gene, and there is an edge between two vertices if an interaction has been found between the corresponding proteins. We annotate each node in the network with the IDs of the individuals having one or more mutations in the corresponding gene. Our goal is to find a relatively small connected component $G^{\prime}$ such that most patients have mutations in one of the genes within it. A small subgraph is more likely to consist of functionally related genes and is less likely to be the result of overfitting to the set of individuals whose diseases we are analyzing. However, we would also like our model to have the greatest possible explanatory power—that is, to account for, or cover, as many patients as possible by including genes that are mutated within their cancers. We formulate our problem to balance these two competing objectives with a parameter $\alpha$ that controls the trade-off between keeping the subgraph small and covering more patients.

For a fixed value of $\alpha$, we have developed two approaches to solve the underlying optimization problem. One is based on linear programming and the other is a fast greedy heuristic (see below). We use the greedy heuristic in the context of a carefully designed cross-validation procedure to automatically select a value for $\alpha$ that results in good coverage of patients but avoids overfitting to them. Once $\alpha$ is selected, this value is used within our objective function and we next analyze the entire patient cohort. In particular, multiple independent trials using $\alpha$ are run on randomly chosen subsets of the patient data, as we have found that introducing a little bit of randomness helps increase performance as compared to a single run on the full data set. Each trial outputs a subgraph, and our final aggregated output is an ordered list of candidate genes ranked by how frequently each has been selected over the trials.

General formulation. Each vertex $v_{j}$ is associated with a set $C_{j}$ containing the IDs of the individuals who have somatic mutations in the corresponding gene. A patient with ID $i$ is covered if $i\in\bigcup\limits_{v_{j}\in G^{\prime}}C_{j}$, and uncovered otherwise. We formulate our problem as that of finding a connected subgraph $G^{\prime}$ of $G$ so as to minimize

[TABLE]

where $X$ is the fraction of patients that do not have an alteration in a gene included in $G^{\prime}$ (i.e., they are uncovered), $Size(G^{\prime})$ is the size of the subgraph, and $0\leq\alpha\leq 1$ is a fixed parameter controlling the trade-off between keeping the subgraph $G^{\prime}$ small and covering more patients. We note that our problem is similar, though not identical, to the Minimum Connected Set Cover Problem [31], a NP-hard problem.

A simple and natural measure for the size of a subnetwork is its number of nodes (i.e., $Size(G^{\prime})=|G^{\prime}|$). However, longer genes may tend to acquire more mutations simply by chance. We correct for that by associating with each node $v_{j}$ a weight $w_{j}$ that is equal to the ratio of the length of the gene to the total number of mutations it has. The size of the subcomponent is then defined as $Size(G^{\prime})=\sum\limits_{v_{j}\in G^{\prime}}w_{j}$. This way, genes having longer length will be weighted more, correcting for a possible bias towards selecting longer genes. We note that since our objective function balances the fraction of uncovered patients with the size of the graph, we would like the size of the graph to be between 0 and 1; thus, we normalize each node weight by dividing by the unnormalized size of what we call a fully covering subgraph $G^{f}$—a connected subgraph of $G$ that covers all patients. (In practice, we compute $G^{f}$ using the greedy heuristic described below, with $\alpha=1$).

Integer linear programming formulation. The problem of finding a minimum connected subgraph that covers as many patients as possible can be solved using constraint optimization. Let $n$ be the number of patients in our sample. For each patient $i$, we define a binary variable $p_{i}$ that is set to 1 if patient $i$ is covered by the chosen subgraph $G^{\prime}$, and 0 otherwise. For each vertex (or gene) $v_{j}$, we define a binary variable $x_{j}$ that is set to 1 if the vertex is included in the chosen subgraph $G^{\prime}$, and 0 otherwise. It is straightforward to set up constraints to ensure that a patient is considered uncovered if none of its mutated genes are part of $G^{\prime}$, and covered if at least one of its mutated genes is selected as part of $G^{\prime}$ (see Equations (1) and (2) below).

The challenging part of the ILP is setting up constraints to ensure that the chosen nodes form a connected subgraph $G^{\prime}$. For this task, we employ a flow of commodity technique [11], which we now briefly describe. We inject $|G^{\prime}|$ units of flow into $G^{\prime}$ (i.e., we inject $\sum x_{i}$ units of “flow” into a vertex that is included in the chosen subnetwork). Flow can move from one vertex to any of its neighbors in the network, and each vertex removes exactly one unit of flow as the flow passes through it. All flow must be removed from the subnetwork, and we set the constraints so that this is possible only if the subnetwork $G^{\prime}$ is connected. For the source of the flow we use an artificial external node $v_{extr}$. The main issue is that we do not know which node $v_{extr}$ should be connected to, as we do not know the nodes of $G^{\prime}$ in advance. To resolve this, we decide that $v_{extr}$ connects to the node that covers the largest number of patients $v_{max}$; this is equivalent to determining in advance that $v_{max}\in G^{\prime}$, though as an alternate approach we could also decide to choose this node probabilistically and run the ILP several times. Finally, to handle the flow constraints, for each edge $(i,j)\in E$, we introduce integer variables $y_{i,j}$ and $y_{j,i}$ to represent the amount of flow from node $i$ to node $j$ and from node $j$ to node $i$, respectively. The full integer linear program is:

*minimize *$\alpha(n-\sum\limits_{i}p_{i})/n+(1-\alpha)\sum\limits_{j}x_{j}w_{j}$

subject to

[TABLE]

Equation (1) ensures that a patient is considered covered if one of his or her or somatically mutated genes is included in $G^{\prime}$. Equation (2) ensures that a patient is not considered covered if none of his or her somatically mutated genes is chosen to be part of the subgraph. Equations (3), (4) and (5) enforce the connectivity requirement. Equation (3) requires that the flow going out of each vertex in the chosen subnetwork is 1 less than the flow coming in. Equation (4) requires that if a vertex is not part of the chosen subgraph, the flow going through it is 0. Equation (5) sets the amount of flow injected into the subgraph to be equal to the number of chosen nodes.

Greedy heuristic. Solving the ILP yields an exact solution but is computationally difficult. Thus, we have also developed an efficient greedy heuristic. Our heuristic procedure initializes $G^{\prime}$ by randomly choosing the first gene from among the five most mutated genes, with probability proportional to the number of patients it is found mutated in. It then expands the subgraph $G^{\prime}$ iteratively as follows. At each iteration, all vertices that are at most distance 2 from a vertex in $G^{\prime}$ are examined and the one that improves the objective function the most is chosen; any ties are broken uniformly at random. If this vertex is not directly adjacent to the nodes in the subnetwork, the intermediary node is also added. The heuristic terminates when no improvement to the objective is possible. We repeat this heuristic multiple times, as it is probabilistic.

In practice, the greedy heuristic finds a solution that is on average $\sim$90% of the best value for the objective function as determined by the ILP formulation using CPLEX [19]. For example, on the glioblastoma dataset of 277 individuals, the ILP finds 61 genes covering 90% of the patients when using $\alpha$=0.5. In comparison, for this value of $\alpha$, the greedy heuristic finds on average 66 genes covering 88% of the patients. In the rest of the paper, we use the greedy optimization as it has comparable performance to the ILP, while being much faster.

Parameter selection and solution aggregation. Our approach to uncover a subnetwork of mutated genes that covers many patients has a single parameter, $\alpha$. Large values of $\alpha$ result in a larger number of selected genes that cover more patients, yet may contain more irrelevant genes; this may especially be a factor in the current analysis if there are many samples where missense mutations are not the driving event. We devised a simple but effective data-driven cross-validation technique to choose an appropriate $\alpha$ for a set of cancer samples. In particular, we split our samples into training, validation and test sets [16]. A test set of (10%) of the patients is completely withheld. While varying $\alpha$ in small increments in the interval $(0;1)$, the remaining data is repeatedly split (100 times for each value of $\alpha$) into training (80%) and validation (20%) sets. For each split, the greedy heuristic algorithm is run on the training set to find $G^{\prime}$. The fractions of patients covered (by the selected $G^{\prime}$) in the training and validation sets are compared. The parameter $\alpha$ is selected where performance on the validation sets deviates as compared to the training sets. While this can be done visually, for all results reported here we do this automatically using a simple two-rule procedure that selects the smallest $\alpha$ for which the difference in average coverage between the training and validation set exceeds $5\%$ and for which average performance on the validation set is within $10\%$ from the maximum observed one for any $\alpha$. Finally, the coverage of patients on the (completely withheld) test set is computed to ensure it is similar to the one on the validation set.

Once $\alpha$ is chosen for a set of cancer samples, we repeatedly (1000 times) run the algorithm on this set, each time withholding a fraction ($15\%$) of the patients in order to introduce some randomness in the process. Genes are then ranked by the number of times they appear in $G^{\prime}$. In practice, we have found that this improves performance as compared to running the algorithm once on the full data set.

Data sources and pre-processing. We downloaded all level 3 cancer somatic mutation data from The Cancer Genome Atlas (TCGA) [35] that was available as of October 1, 2014. This data consists of a total of 19,460 genes with somatic point mutations across 24 cancer types. For each cancer, samples that are obvious outliers with respect to their total number of mutated genes are excluded. See Supplementary Table 1 for a list of the cancer types, the cancer-specific thresholds to determine outlier samples, the number of patient samples considered for each cancer type, and other statistics about the TCGA somatic mutation dataset.

We use two different biological networks in our analysis: HPRD [30] (Release 9_041310) and BioGrid (Release 3.2.99, physical interactions only) [33]. Biological networks can exhibit several nodes with very high connectivity, often due to study bias. As such high connectivity destroys the usefulness of the network information, we remove all nodes whose degrees are clear outliers. For BioGrid, this removes UBC, APP, ELAVL1, SUMO2, CUL3. For HPRD, this removes no nodes. For both networks, we exclude the nine longest genes (TTN, MUC16, SYNE1, NEB, MUC19, CCDC168, FSIP2, OBSCN, GPR98) as they tend to acquire numerous mutations by chance while covering many patients.

To further handle the connectivity arising within the networks due to high-degree nodes, we filter edges using the diffusion state distance (DSD) metric introduced in [3]; the DSD metric captures the intuition that edges between nodes that also share interactions with low degree nodes are more likely to be functionally meaningful than edges that do not (and thus are assigned closer distances). For each edge, the DSD scores (as computed by the software of [3]) between the corresponding nodes are Z-score normalized, and edges with Z-scores $>0.3$ are removed. We note that the overall performance of our approach improves when performing this filtering (data not shown), supporting the claim of [3] that preprocessing a biological network in this manner is an important step. The final number of nodes and edges, respectively, for the filtered networks are 9,379 and 36,638 for HPRD; and 14,326 and 102,552 for BioGrid.

Other approaches. To ascertain the contribution of network information, we compare nCOP to two approaches that do not use network information: (1) MutSigCV 2.0 [25], a state-of-the-art method that identifies genes that are mutated more frequently than expected according to a background model, and (2) a set cover approach that tries to find mutated genes that simply cover as many patients as possible. We formulate the set cover approach as an ILP that tries to find a good cover consisting of $k$ vertices. Using the same notation as for nCOP, the set cover objective is to *maximize *$\sum\limits_{i}p_{i}$, subject to Equations (1) and (2) of nCOP, and with the additional constraints that $\sum\limits_{j}x_{j}\leq k$ and $\sum\limits_{j}x_{j}\geq k$. We also compare nCOP to HOTNET2 [26] and Muffinn [5], two recent network-based approaches that aggregate mutational information. All comparisons are made when running on exactly the same cancer and (if used) network data. MutSigCV 2.0, Hotnet2 and Muffinn are run with default parameters.

Performance evaluation. To evaluate the gene rankings of all the tested methods, we use the curated list of 517 cancer census genes (CCGs) available from COSMIC [13]. All genes in this list are considered as positives, and all other genes are considered as negatives. Though we expect that there are genes other than those already on the CCG list that play a role in cancer, this is a standard approach to judge performance (e.g.,see [20]) and gives us an idea of how methods are performing as cancer genes should be highly ranked by methods that perform well. Since only the top predictions by any method are relevant for cancer gene discovery, we judge performance by computing the area under the precision-recall curve (AUPRC) using the top 100 genes predicted by each method. If a method returns less than 100 genes, we extend the precision-recall curve to 100 genes assuming that it performs as a random classifier. We note that reasonable changes to the number of predictions considered does not change our overall conclusions (data not shown).

Results

We run nCOP, using the greedy heuristic algorithm, on somatic point mutation data from 24 different TCGA cancer types. Results in the main paper use the HPRD network [30] for all analysis and show some results using kidney renal clear cell carcinoma (KIRC) with 416 samples as an examplar.

Increasing patient coverage while preventing overfitting. We first demonstrate that, across the 24 cancer types, our training-validation-test set framework is a highly effective approach for choosing an $\alpha$ that balances patient coverage with subnetwork size. For all cancers, when using somatic missense mutation data, as $\alpha$ increases, the total number of genes in the chosen subnetwork $G^{\prime}$ increases (as expected), as does the fraction of patients in the training set that are covered by these genes (Figure 1a and Supplementary Figure S2). For smaller values of $\alpha$, coverage on the validation sets closely matches that obtained on the training sets; that is, the sets of genes chosen using patients in the training sets are also effective in covering patients in the corresponding validation sets. For KIRC, when $\alpha=0.5$, genes chosen using the training sets cover on average nearly 70% of patients in the corresponding validation sets. The fact that a small subnetwork can be found that covers a large fraction of previously unseen patients is consistent with the hypothesis that a shared pathway or process plays a role in most (but not all) of these patients’ cancers.

For larger values of $\alpha$ ($>0.6$ for KIRC), however, coverage on the validation sets lags behind that observed on the training sets. For even larger values of $\alpha$ ($>0.85$ for KIRC), the algorithm selects many genes, and eventually increases the coverage for most cancers on the training sets to nearly 100%. However, larger values of $\alpha$ do not substantially increase coverage of the withheld patients. This difference between the training and validation curves captures the overfitting of the model and also illustrates the trade-off between covering more patients and keeping the solution parsimonious. We note that the eventual plateau of the validation curve is consistent across cancer types (Supplementary Figure S2) and on different networks (data not shown). For each cancer type, values of $\alpha$ are selected by our automated procedure (see Methods); this value is $\alpha=0.5$ for the KIRC dataset shown in Figure 2a.

As a control, we repeat the same procedure using only synonymous mutations (Figure 1b). We observe that the coverage on the validation sets is much poorer. Though coverage of course increases as more nodes are added, it never exceeds $50\%$ even when $\alpha$ is increased to $1$ or when we have nearly perfect coverage on the training set, despite adding many more nodes. This poor performance is consistent with the expectation that synonymous mutations do not result in altered protein sequences and do not disturb cellular pathways. Hence, given the differences observed between using missense versus silent mutation data when varying settings for $\alpha$ and comparing training and validation sets, our formulation appears to be well-suited for extracting cancer-relevant information from mutational profiles and interaction networks.

nCOP effectively uses network information to uncover known cancer genes. Having shown in the previous section how to select a value for the only parameter in the model and having demonstrated that our formulation is effective in choosing genes that are mutated in previously unobserved patient samples, we next evaluate nCOP’s performance in predicting cancer genes.

We first consider how well nCOP recapitulates known cancer genes (CCGs) [12]. We find that on KIRC, our top predictions include a high fraction of CCG genes (Figure 2). For example, genes that are always output are VHL, PTEN, and BAP1, three well-known cancer genes. Other known cancer genes such as JUN, BLM and ARID1A, are also highly ranked.

We next illustrate the power of our network-based method by comparing its performance on the KIRC data set to approaches that do not consider network information (Figure 2). First, we consider a set cover version of our approach that does not use network information at all. We find that for the same number of predicted genes, our approach consistently has a larger fraction of CCGs, demonstrating the advantage of using network information. We also considerably outperform a state-of-the-art frequency-based approach, MutSigCV 2.0 [25].

We next compare nCOP to these two non-network approaches across all 24 cancer types. In particular, we compute the $\log$ ratio of the area under the precision recall curve (AUPRC) of our approach versus each of the other approaches on each cancer type (Figure 3). We outperform MutSigCV 2.0 in 21 of the 24 cancers and the set cover approach in all cancers, demonstrating the clear advantage of using network information; the performance improvement of nCOP over the set cover approach is particularly notable as the main difference of these approaches is the additional use of network information by nCOP. In several cancers, the performance improvements of nCOP are substantial. For example, nCOP shows a four-fold improvement over MutSigCV 2.0 in predicting cancer genes for liver hepatocellular carcinoma (LIHC) and a nearly eight-fold improvement over MutSigCV 2.0 on pheochromocytoma and paraganglioma (PCPG). Importantly, nCOP uncovers rarely mutated CCGs genes that network-agnostic methods fail to. For instance, KLF6 and TCF7L2, mutated respectively in only one and two individuals in LIHC, and WT1, mutated in two individuals in PCPG, are all uncovered by nCOP because of their proximity to other mutated genes in the network. In contrast, frequency-driven MutSigCV 2.0 and network-agnostic set cover both fail to uncover these CCG genes. We note that the suboptimal performance of nCOP in the case of uterine carcinosarcoma (UCS) is due to the fact that for this particular data set, two genes TP53 and FBXW7 cover $85\%$ of the patients, with only four genes being sufficient to cover $95\%$ of the patients. This renders nCOP incapable to effectively leverage network information as it returns only four genes. Nevertheless, the overall superior performance of nCOP as compared to these two non-network based approaches on the vast majority of cancers demonstrates its considerable power.

Having shown that nCOP better identifies cancer-relevant genes than two approaches that do not use network information, we next consider whether the specific way in which nCOP uses network information is beneficial. Towards this end, we compare the effectiveness of nCOP in uncovering cancer genes to Muffinn [5], a method published within the last year that considers mutations found in interacting genes. We find that in 20 out of the 24 cancer types, nCOP outperforms Muffinn (Figure 4). We also compare nCOP to Hotnet2  [26], a cutting-edge network diffusion method. As Hotnet2 does not output a ranked list of genes, we could not compute the area under the precision recall curve. Instead, examining the complete list of genes highlighted by both methods, we observe that nCOP exhibits significantly better precision while trailing slightly in recall (Supplementary Figure S3).

nCOP reveals novel genes, including those that are rarely mutated. In addition to ranking known cancer genes highly, nCOP also gives high ranks to several non-CCG genes that may or may not be implicated in cancer, as our knowledge of cancer-related genes is incomplete. We observe that non-CCG that are highly ranked by nCOP tend to be less frequently mutated. For example, among these novel predictions for KIRC are SALL1, NR5A2, and UBE2A which have all recently been implicated to play a role in cancers [43, 27, 45].

When we consider the proteins output by our procedure in at least half the runs for KIRC, we find that they are enriched in many KEGG and GO pathways relevant for cancer, including MicroRNAs in cancer, Choline metabolism in cancer, Central carbon metabolism in cancer and PI3K-Akt signaling pathway (Bonferroni-corrected $p<0.001$, hypergeometric test). When we remove all known CCG genes and consider only new predictions, thus performing a much harder test, only one pathway, Thyroid hormone signaling pathway, is enriched. Interestingly, it has been shown that thyroid hormones play a role in kidney growth and development [22] and four of our non-CCG predictions are part of that pathway, together with four known cancer genes. Interestingly, two of the non-CCG genes we find in this pathway are rarely mutated, and are not highlighted by MutSigCV 2.0 or set cover. This illustrates the power of nCOP to zoom in on potentially cancer-relevant modules consisting of rarely mutated genes.

Potential biases, performance on randomized data, and robustness. Due to space and time constraints, we briefly describe some additional tests we performed to show that nCOP is robust and well-behaved. First, we have tested nCOP on another network and have shown that it remains effective in identifying CCGs (Supplementary Figures S4). Second, to confirm the importance of network structure to nCOP, we have run nCOP on two types of randomized networks, degree-preserving and label shuffling, and have shown that (as expected) overall performance deteriorates across the cancer types (Supplementary Figure S5); we note that though network structure is destroyed, per-gene mutational information is perserved, and so highly mutated genes are still output. Third, to make sure that the novel genes we uncover are not driven by patients with large numbers of passenger mutations (i.e., that the novel genes are not likely to be passenger genes), we have compared the overall number of mutations for patients having missense mutations only in CCG genes but not in any predicted non-CCG (or novel) genes to the total number of such mutations for patients having missense mutations only in novel genes but not in any CCG genes (Supplementary Figure S6), and have found that patients with only mutations in novel genes do not harbor more mutations. Finally, to make sure that genes are not more likely to be picked because they have higher degree, we have confirmed that newly predicted do not tend to exhibit higher degree than known cancer genes; indeed, among all novel genes found across all cancer types, most have degree between $<15$, and there are only a couple with high degree ($\geq 50$).

Discussion

In this paper, we have shown that nCOP, a method that incorporates individual mutational profiles with protein–protein interaction networks, is a powerful approach for uncovering cancer genes. Our method is based on an intuitive mathematical formulation and demonstrates higher precision than other state-of-the-art methods in detecting known cancer genes. Further, our approach is particularly beneficial in highlighting infrequently mutated genes that are nevertheless relevant for cancer. Our approach therefore complements existing frequency-based methods (e.g., [25]) that generally rely on comparisons to background mutational models and lack the statistical power to detect genes mutated in fewer individuals.

In the future, nCOP can be extended in a number of natural ways. First, while nCOP currently analyzes only mutations within genes that affect protein coding, other alterations are also commonly observed in cancers. For example, copy number variants (CNVs) are found frequently in cancers and can play critical functional roles  [46]. Although nCOP does not currently use CNV information, our framework can be extended to incorporate this data. Indeed, as the numbers of CNVs and point mutations found within each cancer genome appear to be inversely related [8], considering both types of alterations will increase the power of our approach. Second, nCOP may also benefit from incorporating gene weights that reflect its likelihood to play a role in cancer; in our current work, we consider a gene’s length but no other gene-specific attributes are considered. Such gene weights may be derived from existing approaches to detect frequency of mutation or to assess the functional impact of mutations. Finally, while nCOP can output groups of genes that are not part of a single connected component due to our randomized aggregation procedure, extending nCOP’s core algorithms to explicitly consider multiple subnetworks corresponding to distinct pathways may be a particularly promising avenue for future work.

We conclude by noting that researchers can use our framework to rapidly and easily prioritize cancer genes, as nCOP requires only straightforward inputs and runs on a desktop machine. Indeed, nCOP’s efficiency, robustness, and ease of use make it an excellent choice to investigate cancer as well as possibly other complex diseases. As sequencing costs plummet and cancer and other disease sequencing mutational data become more abundant, the predictive power of our method should only increase. In sum, we expect that our method nCOP will be of broad utility, and will represent a valuable resource for the cancer community.

Supplementary Tables and Figures

The following pages contain a table summarizing the TCGA data we use along with 6 supplementary figures that support the findings of the main paper.