extfieldsofstudy:
Computer Science
Biology
Medicine

s2fieldsofstudy:
Computer Science
Biology

text:
A new correlation clustering method for cancer mutation analysis

Cancer genomes exhibit a large number of different alterations that affect many genes in a diverse manner. It is widely believed that these alterations follow combinatorial patterns that have a strong connection with the underlying molecular interaction networks and functional pathways. A better understanding of the generative mechanisms behind the mutation rules and their influence on gene communities is of great importance for the process of driver mutations discovery and for identification of network modules related to cancer development and progression. We developed a new method for cancer mutation pattern analysis based on a constrained form of correlation clustering. Correlation clustering is an agnostic learning method that can be used for general community detection problems in which the number of communities or their structure is not known beforehand. The resulting algorithm, named $C^3$, leverages mutual exclusivity of mutations, patient coverage, and driver network concentration principles; it accepts as its input a user determined combination of heterogeneous patient data, such as that available from TCGA (including mutation, copy number, and gene expression information), and creates a large number of clusters containing mutually exclusive mutated genes in a particular type of cancer. The cluster sizes may be required to obey some useful soft size constraints, without impacting the computational complexity of the algorithm. To test $C^3$, we performed a detailed analysis on TCGA breast cancer and glioblastoma data and showed that our algorithm outperforms the state-of-the-art CoMEt method in terms of discovering mutually exclusive gene modules and identifying driver genes. Our $C^3$ method represents a unique tool for efficient and reliable identification of mutation patterns and driver pathways in large-scale cancer genomics studies.

Introduction
Rapid advances in high-throughput sequencing technologies have provided unique opportunities for analyzing large numbers of cancer genomes. However, the complexity of genomic alterations in cancer remains a challenge that has to be overcome in order to fully characterize the functional roles of various mutations. Cancer genomes often exhibit a large number of different mutations that affect genes in diverse manners. But the vast majority of these mutations do not have significant impact on tumorigenesis [18]. A central question in cancer genomics is how to distinguish "driver" mutations, which contribute to tumorigenesis, from functionally neutral "passenger" mutations. Such driver mutations (e.g., point mutations or copy number changes) are of critical importance in elucidating key biological pathways perturbed in cells that eventually lead to cancer.
Many computational methods have been developed to facilitate the discovery of driver mutations [6,11,16,22,26]. Due to the high level of inter-tumor heterogeneity, two patients with the same cancer may have vastly different drivers and most cancer mutations occur in very low frequency in the patient population. Therefore, approaches relying on simple recurrence or frequency of mutations do not work well in practice. Recently, pathway-based and network-based models have been shown to be effective not only in determining common driver mutations and mutation patterns, but also in pinpointing the key biological pathways and subnetworks affected by driver mutations [4,20,[29][30][31]. Such methods have a unique advantage in so far that in addition to mutation analysis, they take into account gene interactions as an added source of prior knowledge.
In parallel, a number of methods have been proposed to identify driver pathways, i.e., groups of genes that may interact together in combinatorial patterns to promote tumorigenesis. The authors of (9) described a method called MEMo, and subsequently used it to show that mutually exclusive modules based on known networks can aid in determining groups of genes that contribute to tumorigenesis. These gene groups, or modules, are jointly highly recurrent, have similar pathway impact in terms of biological processes, and their corresponding mutations tend to be mutually exclusive, meaning that very often only one gene in each gene group is mutated at a given time in any given patient. This mutual exclusivity rule in cancer pathways is supported by the observations that, in general, one mutated gene suffices to perturb the function of its corresponding pathway. Multiple mutations would require significantly higher energy investments on the part of cancer cells, and are hence selected against. Dendrix [41] (and later Multi-Dendrix [23]) was developed to identify driver pathways de novo using mutual exclusivity and coverage (recurrence) principles, without relying on known network information that has the potential to improve the discovery process of new modules. More recently, CoMEt [24] was proposed to address an inherent bias in Dendrix and Multi-Dendrix that resulted in high frequency mutations being significantly more likely to be included in mutually exclusive modules.
However, while Dendrix, Multi-Dendrix and CoMEt have the ability to identify mutually exclusive modules de novo, there still have significant limitations. These methods are typically very inefficient when it comes to applying them on large-scale datasets with large parameter setting. Also, some of these methods are randomized in nature and no guarantees exist that multiple runs of the methods will produce compatible results. Furthermore, almost all methods are able to identify only a small number of modules with limited number of genes, as cluster sizes are critical algorithmic parameters from the perspective of computational tractability.
To overcome these and other shortcomings of existing methods, we introduce a novel method called Cancer Correlation Clustering (C 3 ) to directly tackle the problems of integrating diverse sources of evidence regarding driver pattern behavior and eliminating computational bottlenecks associated with large cluster sizes or cluster numbers. C 3 uses an optimization framework specifically developed for the driver discovery task, where data is converted to a simple set of optimization weights that do not require the algorithm to change upon incorporation of new data sources. In addition to this flexibility, C 3 has low computational cost, and it allows for adding relevant problem constraints while retaining good theoretical performance guarantees.
The paper is organized as follows. A basic introduction of the principles of correlation clustering is provided in Section 2. Section 3 contains a description of the weight computation methods, the algorithmic clustering approach based on the computed weights, and the evaluation criteria used to compare C 3 and CoMEt. Section 4 contains the main results of our analysis, contrasting the performance of C 3 and CoMEt on breast cancer and glioblastoma data. A discussion of our findings and concluding remarks are given in Section 5. A rigorous mathematical performance analysis of C 3 may be found in the Supplementary Materials.

Approach
The basic idea behind the C 3 approach is correlation clustering, an agnostic learning technique first proposed in [3]. In the most basic form of the clustering model, one is given a set of objects and, for all or some pairs of objects, one is also given an assessment as to whether the objects are "similar" or "dissimilar". This information is described using a complete graph with labeled edges: each object is represented by a vertex of the graph, and the assessments are represented by edges labeled with either a "+" symbol, for similar objects, or a "-" symbol, for dissimilar objects. The goal is to partition the objects into clusters so that the edges within clusters are mostly positive and the edges between clusters are mostly negative. Unlike in many other clustering models, such as k-means [19], the number of clusters is not fixed ahead of time and finding the optimal number of clusters is part of the problem.
The similarity assessments need not be mutually consistent: for example, if the graph contains a triangle with two positive edges and one negative edge, then we must either group the endpoints of the negative edge together, erroneously putting a negative edge inside a cluster, resulting in a "negative error" or else we must group them separately, forcing one of the positive edges to erroneously go between clusters, resulting in a "positive error". When a perfect clustering is not possible, we seek an optimal clustering: one that minimizes the total number of "errors." This form of correlation clustering is known to be NP-hard, but depending on the graph topology, various constant or logarithmic approximation guarantees exist.
The authors of [3] also proposed a weighted version of the correlation-clustering problem, where the edges of the graph receive weights between −1 and 1 rather than simply receiving + or − labels: an edge with weight x incurs cost 1+x 2 if it is placed between clusters and cost 1−x 2 if it is placed within a cluster. A more general weighted formulation was introduced in [7,8], and this is the formulation we subsequently consider. In this model, each edge e is assigned two nonnegative weights, w + e and w − e . A clustering incurs cost w + e if e is placed between clusters, and incurs cost w − e if e is placed within a cluster. If no restrictions are placed on the weights w + e and w − e , then it is possible to have edges with w + e = w − e = 0; these edges are effectively absent from the graph, so there is no loss of generality in assuming that the graph is a complete graph. In order to arrive at problems that have efficient constant approximation algorithms, one needs to place certain restrictions on w + e and w − e . The probability constraints give a natural restriction on the edge weights w + e + w − e = 1 for every edge e. Another widely studied restriction is the triangle inequality restriction, where one requires w − uw ≤ w − uv + w − vw for all distinct vertices u, v and w. The analytic approach pursued in this work operates on the following model: genes which show sufficiently large mutation prevalence in cancer patients represent vertices of the complete graph to be clustered. Note that in this work we only use the top 5% of genes ordered by mutation frequency. This equates to 170 genes in glioblastoma (GBM) and 130 genes in breast cancer (BRCA). The weights w + e and w − e of an edge e connecting two genes g 1 and g 2 are weighted sums of the mutual exclusivity and coverage strength, as well as an adequately chosen measure of network distance and expression similarity.
More precisely, the negative weights w − e are chosen to be relatively small if the endpoint genes describing the edge are deemed to be mutually exclusive in cancer patients. A small negative weight encourages placing mutually exclusive genes within the same cluster. The positive weights jointly depend on the coverage, network distance and expression correlation of the endpoint genes: the larger the joint coverage, co-expression and inverse of the network distance of the endpoint genes, the more likely they will end up in the same cluster. Precise mathematical formulations of the weight functions will be provided in the next section.
To control the size of the resulting clusters so as discourage uninformative singleton and giant clusters, we developed two new correlation clustering algorithms that use cluster sizes as problem parameters that may be chosen by the users. These cluster size bounds also allow for more accurate comparison with other methods which operate with inherent cluster size constraints. Notice that unlike in the aforementioned known methods, the cluster sizes have no bearing on the complexity of our algorithm nor on their overall approximation quality.
The driver discovery approaches closest to C 3 are Multi-Dendrix [23] and CoMEt [24]. Multi-Dendrix is an integer linear programming clustering algorithm that ensures that the genes within a cluster have mutation patterns that satisfy mutual exclusivity and coverage: for any two genes in a cluster, the number of patients in which these genes are mutated at the same time is relatively small; in addition, a large portion of the patients should have at least one mutation in each cluster. CoMEt uses a statistical score of mutation exclusive that conditions upon the frequency of each alteration, alleviating the inherent bias caused by frequently mutated genes. Compared to Multi-Dendrix, C 3 uses a weighted linear programming relaxation instead of an integer linear program which significantly improves the versatility and running time of the algorithm. Furthermore, the weights allow for straightforward incorporation of heterogeneous sources of evidence into the clustering method and the algorithm itself remains unchanged with the addition of new data. This flexibility comes at the cost of C 3 providing only an approximate solution. Nevertheless, given the inherently approximate nature of weight selection and parametrization of both algorithms, this does not appear to be a significant shortcoming. Also, empirical evaluations suggest that the approximation algorithms produce results very close to the optimal solution. Another advantage is that if one needs to change the combinatorial conditions that the clusters satisfy, Multi-Dendrix cannot be easily adapted, and a whole new algorithm needs to be developed. On the other hand, our algorithms are very flexible and most combinatorial patterns can be easily incorporated in the same algorithm, by simply changing the weights assigned to the edges.

Methods
Before rigorously describing our algorithmic methods, we introduce some relevant notation and explain how to estimate appropriate clustering weights based on available data.

Clustering Weights
Let G(V, E) be a complete graph, where V (G) denotes the set of vertices and E(G) denotes the set of edges of the graph G, respectively. The symbol e ∈ E(G) or e = uv with u, v ∈ V (G) is used to denote a generic edge. Each edge is assigned a positive weight w + e and a negative weight w − e . Recall the interpretation of these weights: for two distinct vertices u, v ∈ V (G), w + uv is the cost of placing u and v in different clusters; consequently, by making the positive weight of an edge large, one can discourage placing the corresponding two genes into different clusters. Similarly, w − uv is the cost of placing u and v in the same cluster, and hence making this weight large discourages placing the corresponding two genes into the same cluster. In the rest of this section, we will explore different ways of defining the weights w − uv and w + uv ; in order to avoid confusion between the different definitions, each weight we define will include a parenthetical abbreviation, so that, for example, w + (c) uv will refer to the positive weight of uv defined according to the coverage criteria, while w + (c,n) uv will refer to the positive weight of uv according to the coverage and network criteria.
The weights are computed using four types of datasets: gene alteration data, gene copy number variation (CNV), network information (NI), and gene expression (GE) data. Let n p denote the number of samples (i.e., patient genomes available) and let n g denote the number of genes. Also, let A ∈ {0, 1} ng×np denote the matrix containing alteration data: if gene i is altered in sample j, we set A(i, j) = 1; otherwise, we set A(i, j) = 0. Also, let C be an n g ×n p matrix representing the CNV data: we set C(i, j) = 0 if there is no change in the copy number of gene i in sample j; otherwise, we choose an integer value reflecting the deviation of the CNV number from its baseline value. Hence, the CNV matrix contains both positive and negative values corresponding to the copy number changes of the corresponding gene in each sample.
To combine the CNV with alteration, we referred to the following method. Using the matrices A and C, we formed a new binary matrix M ∈ {0, 1} ng×np such that and M(i, j) = 1, otherwise. In this formulation, l cnv and h cnv are lower and upper bounds that may be chosen by the user. These bounds determine what is deemed to be a significant CNV change. In our tests, we set l cnv = −1 and h cnv = 3, although other options are clearly possible. It is worth pointing out that more conservative CNV thresholds tend to decrease coverage, while more relaxed CNV assumptions tend to decrease mutual exclusivity. Based on the procedure above, we arrive at one "mutation" matrix M which we use instead of the matrices A and C. Finally, let Z ∈ R ng×np be the matrix corresponding to z-scores of gene expression data: here, Z(i, j) denotes the z-score of the expression of gene i in sample j. More precisely, if the raw expression of gene i in sample j equals x ij , then Z(i, j) = x ij−µ i σ i ; µ i denotes the average expression of gene i and σ i denotes its standard deviation.

Clustering Weights Determined Based on Mutual Exclusivity and Coverage (ME-CO)
The idea behind our approach is to impose the mutual exclusivity constraint through the weights w − e and coverage constraint through the weights w + e s. For each gene (i.e., vertex) u, let S(u) denote the set of patients in which u is mutated. Note that we use the matrix M to determine if a mutation in the gene exists, either due to sequence alteration or CNV. Then, for any u, v ∈ V (G), the negative weights are chosen according to where a is a scaling parameter to be chosen by the user. The intuition behind the use of the factor a is that if a is large (e.g., empirically, a value of a > 3 is deemed large), the mutual exclusivity is enforced strictly, while if a is small, (e.g. a < 3), the genes in each cluster may not be highly mutually exclusive. Also, note that To capture the coverage property through the positive weights, observe that if two genes increase the coverage significantly, their positive weight should be large so that they are encouraged to be placed in the same cluster. To determine the positive weights, we first form the set and ∆ denotes the symmetric difference of two sets. A large value for the symmetric difference D(u, v) suggests that the vertices u and v should be placed in the same cluster, since they increase the coverage of the cluster.
Given the set D, we define T (J) to be the Jth percentile of the values in D. In all our runs, we used the default value of J = 95, although this choice may be governed by the user as well. The positive weights are chosen according to: Note that by this definition, 0 ≤ w + (c) e ≤ 1 for any e ∈ E(G).
In order to ensure that the positive and negative weights meet the constraints imposed by our algorithm needed to ensure a constant approximation guarantee, we require that for all u, v ∈ V (G), w − (e) uv +w + (c) uv ≥ 1. This leads to the additional constraints:

Clustering Weights Determined Based on Mutual Exclusivity, Coverage, and Network Information (NI-ME-CO)
The comprehensive results of pan-cancer studies reported in a number of recent papers (23,32,24,15) have revealed the important connection between network topology and cancer driver distribution patterns. More precisely, the effect of deleterious mutations on the phenotype may be suppressed through a particular configuration of the corresponding protein complexes, and at the same time, the strength of the effect of a mutation may be emphasized through another configuration. As an example, most of the variants observed in healthy individuals seem to appear at the periphery of the interactome, and they do not seem to influence network connectivity. In contrast, cancer driver somatic mutations tend to occur in central, internal regions of the interactome and within highly co-integrated components. It appears that no previous attempts were made to more precisely quantify the network distances between driver variants, which prompted us to perform the following analysis. We first computed the pairwise (shortest) network distances between genes in a large pathway comprising 8, 726 genes from Ciriello et al. [9] via an implementation of the standard Dijkstra algorithm [38]. In this test, we randomly selected 1, 000 pairs in order to reduce the computational burden of running Dijkstra's algorithm O(8726 2 ) times. By using the most complete known driver list from the Cancer Gene Census (CGC) [13], we computed the same distances for driver genes, this time for all pairs of genes. The resulting distribution of shortest paths is depicted in Fig. 1. One can clearly observe that the average shortest distance between drivers is significantly smaller than the average shortest distance between two randomly selected genes. A permutation test confirms this observation, and we calculated a p-value of less than 0.001. These findings suggest that when determining potential driver mutations, one should make use of network distance and connectivity information. This may be accomplished in our approach by adjusting the positive weight of edges connecting two genes: if both endpoint genes were to be drivers, they should be sufficiently central to a given pathway, close to other drivers or to each other.
For this purpose, we consider an undirected graph corresponding to the gene network, denoted by G ; in this graph, which is assumed to be known a priori and which in this work was retrieved from the KEGG Database, each vertex corresponds to a gene. The graph is not complete, but rather relatively sparse, and each edge represents an interaction between genes. As before, we let n p and n g = |V (G)| = |V (G )| denote the total number of patients and the total number of genes in our dataset, respectively. For each vertex u ∈ V (G ), we let N (u) denote the set of neighbors of u and let N (u) = N (u) ∪ {u}. The first step in assigning the positive weights is to determine the set Note suggests that the genes u and v are well connected and likely to be involved in the same pathway [9], and that the corresponding genes should be clustered together.
Given the set F , we define T (J ) to be the J th percentile of the values in F . For any u, v ∈ V (G), the positive weights are then chosen according to the following formula: where w 1 , w 2 ≥ 0, w 1 + w 2 = 1, and as before, the coverage weight equals and the network weight equals Again, in order to ensure that for all The weights w 1 , w 2 may be chosen in such a way as to emphasize the importance of either coverage or network information. We suggest using w 1 = w 2 = 1/2 in a coverage/network only test, although our analysis reveals that emphasizing one criterion over the other offers improved algorithm performance on some datasets.

Clustering Weights Determined Based on Mutual Exclusivity, Coverage, and Gene Expression Data (EX-ME-CO)
Similar to the case of network information, expression data may be incorporated through the positive weights, using the assumption that coexpressed genes may be involved in the same function or cancer pathway. Hence, highly coexpressed genes should be encouraged to cluster together.
To explain how to incorporate gene expression data into the clustering procedure, assume that z(u) and z(v) denote the vectors of time-evolving expression values corresponding to genes u and v, respectively. The first step in assigning the positive weights is to determine the set G = {g(u, v)}, where for every pair of genes u, v, Here, a, b denotes the classical inner product of the vectors a and b, while ||a|| stands for the 2 norm. A large value for g(u, v) indicates that the expression vectors of u and v are highly correlated and hence should be clustered together. Also, note that 0 ≤ g(u, v) ≤ 1 for all u and v.
Given the set G, we let T (J ) denote the J th percentile of the values in G. For any u, v ∈ V (G), the positive weights are chosen according to the following formula: where w 1 , w 2 ≥ 0, w 1 + w 2 = 1, and and Hence, all the conditions are satisfied for the weights, except possibly the third one. In order to make sure that for all u, v ∈ V (G), w − (e) uv + w + (c,x) uv ≥ 1, we include an additional condition that Note that other combinations of datasets may be used, with appropriate changes in the weights. For example, incorporating coverage, network information as well as expression information into a positive weight may be accomplished by setting where w 1 , w 2 , w 3 ≥ 0, w 1 + w 2 + w 3 = 1.

Clustering Algorithms
The bounded cluster size correlation clustering problem for driver gene inference may be formulated as follows. Let K be a "hard" bound on the size of the driver clusters, and let the positive w + and negative weights w − be chosen according to a desired combination of datasets, as described in the previous section. The optimum clustering may be found by solving the integer linear program (ILP) below.
In this formulation, and for a fixed edge e = uv, x uv = 1 implies that u and v should belong to different clusters and x uv = 0 implies that the two vertices should belong to the same cluster. Note that the triangle inequality (17) ensures that if u and z are in the same cluster and z and v are in the same cluster, then u and v are also in the same cluster. Any clustering of the vertices can be described using the variables x e . For a fixed clustering, the objective function is the cost associated with that clustering.
Solving the ILP is NP-hard. We hence relax the problem by changing the integer constraint x e ∈ {0, 1} to an interval constraint x e ∈ [0, 1]. This relaxation leads to a classical LP, the solution of which may be fractional. To obtain a valid clustering, the fractional solutions have to be subsequently rounded to produce integer solutions. Unfortunately, known rounding algorithms (33) for this problem tend to produce small clusters, often as small as single-vertex clusters. For our study, we hence slightly modify the LP by removing the cluster size constraint (19), which we move directly to the rounding procedure as described in Algorithm 1.
Algorithm 1 is closely based on the rounding algorithm described in Charikar et al. [7,8]. The idea behind the rounding algorithm is to pivot on one vertex, examine its closest neighbors, where closeness is governed by the value of the output variables x e of the LP, and partition large neighborhoods if needed to get clusters of size at most K + 1. In the Appendix of the Supplementary Materials, we prove that the LP and Rounding Algorithm 1 provides a 9-approximation for the ILP problem, given that the parameter α is set to 2/7 and given that the weights obey the following constraints: • w + e ≤ 1 for every edge e, and • w + e + w − e ≥ 1 for every edge e. The above inequalities were addressed in the weight selection process through normalization, as described in the previous section.
Let the "pivot vertex" u be an arbitrary element of S.

Evaluation methods
We evaluated the performance of both C 3 and CoMEt in terms of their ability to detect mutually exclusive, high-coverage, and biologically relevant gene clusters. We ran both methods using mutation and CNV data collected from TCGA, pertaining to breast cancer (BRCA) [28] and glioblastoma (GBM) [5]. In addition to GBM and BRCA, we also considered kidney cancer (KIRC) and ovarian cancer (OV), but the available patient data appeared limited at this stage to allow for statistically significant and comprehensive results. We accessed the TCGA provisional data using the cBioPortal platform [14] on August 14, 2015. We ran both methods using the same alteration dataset. We evaluated both point mutations and indels, and for CNVs, we used the GISTIC thresholds [27] of -1 and 3 as our cut-offs. To focus on mutations with high frequency, we only selected genes in the top 95 percentile of alteration frequencies, thereby obtaining 130 genes spanning 959 patient samples in BRCA and 170 genes spanning 291 patient samples in GBM.
To test the effects of cluster sizes and the quality of our results, we ran both C 3 and CoMEt to find clusters of sizes upper bounded by 5, 6, 7, 10, and 15. CoMEt is naturally designed to discover the "most" mutually exclusive gene sets. Due to the fact that correlation clustering will cluster all genes in a dataset, we only compared the top ten most mutually exclusive gene sets generated by C 3 with those of CoMEt.
We ran CoMEt with 1, 000 iterations each and 3 initialization points to ensure both timely and consistent runs. For C 3 , we ran the C 3 clustering method for all combinations of weights w 1 , w 2 , w 3 ∈ {0, 0.25, 0.5, 0.75, 1} that satisfy w 1 + w 2 + w 3 = 1, but selected to report only results for the weight parameters w 1 = 0.167 (coverage), w 2 = 0.333 (network information) and w 3 = 0.333 (expression data). We observe that the choice of the weights may be completely governed by the user, and that the increase in one weight may produce better results in one performance category while reducing the performance in another category. Our choice of high weights for expression and network information was governed by the need to increase the ability of the C 3 algorithm to detect biologically significant clusters. Furthermore, the patient coverage criteria appears to be less relevant than pathway coverage and some other coverage properties that have not been explicitly investigated in the literature. We used four statistical methods to assess the performance of the algorithms which reflect both the statistical and biological significance of the clusters found.
Mutual Exclusivity. To evaluate the degree of mutual exclusivity in a cluster, we performed a Fisher Exact Test [12] for each pair of genes in the cluster. The Fisher Exact Test uses the hypergeometric distribution to calculate the probability of observing a 2 × 2 contingency table of n total samples, with a samples that has an alteration in two genes (say, g i and g j ), b samples with an alteration in gene g i only, and c samples with an alteration in gene g j only. If d is the number of samples with no alteration in either gene, then the probability of co-mutation is evaluated according to a+b a c+d c n a+c (20) We evaluated the overall exclusivity of a cluster as the median value of each pairwise test for each pair of genes g i , g j in the network. The pairwise Fisher's method has also been used by Mutex suite to establish mutual exclusivity [2]. However, because in our context the Fisher Exact Test is used as an evaluation rather than discovery tool, we used the median pairwise p-value rather than the maximum p-value to get a better sense of the overall exclusivity of genes in a cluster. It is also important to note that while CoMEt has an built-in method that generalizes the exclusivity test to 2 k contingency table for a cluster size k, the exponential size of their test set makes evaluation for large cluster sizes computationally impractical. An alternative test for overall mutual exclusivity is a permutation test, as implemented by MEMo, which compares the exclusivity of a gene set by sampling random gene sets and comparing the patients with multiple alterations.
Coverage. To compare and evaluate the overall coverage of a cluster found by C 3 or CoMEt, we calculated the proportion of patients with at least one alteration in a gene belonging to the given cluster.
Network Clustering. We also performed a pathway analysis for the potential drivers. As pointed out in the previous section, driver genes tend to be, on average, closer to each other in a pathway compared to randomly selected genes. Using this dogma, we calculated the average pairwise distance between each each pair of genes g i , g j within a discovered cluster using Dijkstra's Algorithm. As before, our tests were performed on 8726 genes from [9].
Biological Significance. In addition to testing the performance of the algorithm with respect to mutual exclusivity and coverage, we also investigated the biological significance of the C 3 and CoMEt methods from the perspective of gene discovery and pathway analysis. Although there is no overarching gold standard to determine biological significance, a commonly accepted metric employed by MEMo, Dendrix, Mutex, CoMEt and other similar tools is to count the number of known driver genes found within the best clusters according to the mutual exclusivity principle. These clusters usually contain several known driver genes. To determine the driver gene-based biological significance, we calculated the proportion of drivers found in the ten most mutually-exclusive C 3 and CoMEt clusters using a comprehensive, curated list of known drivers from the Cancer Gene Census [13].
It is important to point out that while the four benchmarks we define are a reliable way to test the performance of CoMEt and C 3 , no perfect benchmark exists for detecting mutually exclusive and biologically significant genes clusters. As with many previous methods regarding this topic, the criteria are in a sense circular in that some of the same parameters that we maximize are the parameters that we use to evaluate the method. This is the reason why we use multiple benchmarks to evaluate the method.

Results
In what follows, we demonstrate that C 3 outperforms CoMEt in almost all of the aforementioned benchmarking criteria. As a rule of thumb, C 3 can be made to outperform CoMEt in any chosen single criteria or pairs of criteria by adjusting the weights appropriately. This observation follows from that fact that the weights trade off the strengths of different modeling assumptions. We supplement our statistical analysis with a discussion of the biological relevance of our findings, and explore the role of the new potential drivers found by C 3 within their driver gene communities. In particular, we discuss the significance of large mutually exclusive clusters that cannot be recovered by other methods. Recall that we restrict our attention to the ten best performing clusters according to mutual exclusivity, as this approach was used in the original evaluation process of the CoMEt algorithm.

Performance evaluation
The results of our extensive comparison between C 3 and CoMEt, regarding mutual exclusivity, coverage, driver identification, and pathway-level evaluation, are shown in Fig. 2. Both algorithms were tested on the same server with a 256GB RAM memory. Both methods ran uninterruptedly when the cluster sizes were constrained to k = 5, 6, 7, and 10. CoMet reported segfault memory errors for k = 15, and for this case, only C 3 was benchmarked.
To assess the biological significance of the two methods in terms of their ability to cluster high-impact drivers from the CGC repository together, we compared the results of C 3 and CoMEt both to each other and to a "baseline" value equal to the average proportion of drivers in the ten most mutually-exclusive clusters found, in this case 0.067, using uniform random sampling of genes (see Fig. 2A). In BRCA, we found that C 3 detected a median driver proportion of 0.160 and CoMEt detected a median driver proportion of 0.117 in the top ten clusters. C 3 outperforms CoMEt for each cluster size. We also used a Mann-Whitney Rank Sum test (34) to compare the overall performance of the algorithms with respect to mutual exclusivity, for all cluster sizes. We chose a rank-sum test because it is unclear that the drivers are following a normal distribution due to the small amount of data available. The results show that C 3 outperforms CoMEt (p-value of 0.0079) in terms of amount of drivers in clusters. C 3 also outperforms CoMEt on GBM, with a median proportion of drivers per cluster equal to 0.170, compared to a 0.12 proportion of drivers per cluster found by CoMEt. This finding holds for every cluster size, with a rank-sum test p-value of 0.0361. Both methods succeed in finding biologically significant drivers within clusters exhibiting high mutual exclusivity, and both methods significantly outperform the expected number of drivers per cluster in the random setting (p-value 1.594×10 −5 and p-value 1.312×10 −3 for C 3 and CoMEt, respectively).
We next tested the clusters found by each method based on their mutual exclusivity (see Fig. 2B). To do so, we used the previously described pairwise Fisher's exact test to obtain a p-value for each of the top ten clusters of the two methods. For better visualization, we performed a negative log transform on the p-values, and plotted the transformed p-value distribution. Hence, in this system, larger values indicate more mutually exclusivity. We again used a Mann-Whitney rank-sum test to evaluate the performance of C 3 and CoMEt. For BRCA, one can see that while both methods have significant median exclusivity values (p = 7.541 × 10 −6 for C 3 and p = 3.337 × 10 −4 for CoMEt), C 3 has an overall more significant p-values for each cluster size. The median p-value of C 3 for each cluster size is lower than its CoMEt counterpart except for the case k = 10. However, C 3 does have superior performance overall with a rank-sum p-value of p = 4.0202 × 10 −4 . For GBM, the median exclusivity results are not as strong as for the BRCA set, for both the C 3 and CoMEt method. C 3 has a median p-value of 0.3795 as opposed to CoMEt's 0.5022. The general drop in significance may be attributed to a lower confidence of the Fisher test due to a small number of samples available; recall that the GBM set involved 291 samples, compared to 959 BRCA samples. This indicates that one should look at individual significant clusters to evaluate mutual exclusivity. Even for the reduced median p-value regime, C 3 outperforms CoMEt in significance, having lower median p-values for each cluster size. Overall, the C 3 p-values are consistently and significantly lower than those produced by CoMEt for mutual exclusivity (the rank-sum test p-value equals 0.04401).
The results of the coverage tests are depicted in Fig. 2C. In the coverage benchmark, CoMEt outperforms C 3 for GBM, but neither method outperforms the other BRCA. In BRCA, both methods show comparable performance, with a median result for the fraction of samples covered equal to 0.5505 for C 3 , and 0.5662 for CoMEt. This rather poor performance of both methods is observed for all values of k, with no p-value based on a Student T-test 43 being less than 0.05. The largest difference in coverage recorded for the two methods is present for k = 6. In conclusion, there appears to be no statistical difference between C 3 and CoMEt in terms of BRCA coverage percentage (p-value of 0.5127). In GBM, the median p-value for coverage difference is more pronounced. The median coverage of C 3 is 0.632 and the median coverage of CoMEt is 0.696. CoMEt finds significantly higher-coverage clusters according to the Student T-test, with p-value 0.0345, and the most pronounced coverage percentage differences exist for small values of k (0.3745 vs. 0.6495 for k = 5 C 3 and CoMEt, respectively). It is also important to note the wide distribution of coverage score values produced by C 3 for small k: the IQR (Interquartile range) value is roughly 0.35 for k = 5. The most likely reason behind this result is that our test weights were chosen to boost the relevance of mutual-exclusivity and biological significance rather than coverage. Mutual exclusivity accounts for 100% of the negative weights of edges, while coverage accounts for only 16.7% of the positive weights. We justify this weight choice by the fact that it leads to multiple significant cluster discovery and by our belief that coverage is a less significant driver property compared to mutual exclusivity. We also emphasize that a potentially biologically more relevant coverage constraint should pertain to important pathway, rather than patient sample coverage.
As already mentioned in the previous sections, one advantage of C 3 is that the user can adjust the weights according to her/his own belief about the significance of patient coverage. For example, by changing the averaging weights in our GBM run to w 1 = 0.60 (coverage), w 2 = 0.20 (network), and w 3 = 0.20 (expression), we obtain a coverage percentage of 0.7903 for k = 5. However, this excellent coverage comes at a cost of a less significant mutual exclusivity score (fractional value 0.4288) and a lower proportion of detected drivers (fractional value 0.1267). As may be seen from the above example, C 3 is highly customizable and can be adapted to the user's specification to best reflect the scope and preferences of the analysis.
The last setting in which we analyzed C 3 and CoMEt uses the distances of drivers in the network as performance criteria (see Fig. 2D). Here, we calculated the average pairwise distance between all pairs of genes clustered together. We used the Student T-test to determine the statistical significance of this value. We also compared the values for both algorithms based on 1000 randomly selected genes by using a permutation test. For BRCA, we found no significant performance difference between the two methods in terms of the average pairwise distance: 3.110 for C 3 and 3.070 for CoMEt, with a p-value of 0.9330. In GBM, C 3 showed a smaller average pairwise distance of 2.908 compared to CoMEt's 3.097. This difference is statistically significant, with a p-value of 0.0379. The small average network distance results of C 3 for GBM, coupled with the low coverage, leads to the conclusion that C 3 favors niche, exclusive clusters in biologically relevant cancer pathways. Hence, the method may be useful for discovering specific molecular cancer subtypes. Both methods had an average pairwise distance well below the permutation benchmark of 3.903: the p-values of both C 3 and CoMEt were less than 2 × 10 −16 for each cancer type.
In conclusion, from our detailed evaluation we conclude that although C 3 does not outperform CoMEt with respect to all four evaluation criteria, but only three of them, the C 3 performance indicates a strong overall propensity to select biologically more relevant and mutually exclusive clusters, and with a higher degree of flexibility compared to CoMEt.

Discovering potential driver pathways
We examine next the potential of the C 3 algorithm to detect clusters whose genes may be novel cancer driver candidates. We focus our search on clusters that contain biologically significant driver genes and known biological network interactions, and exhibit high mutual exclusivity and coverage. We also focused on the large cluster size regime, as results in this domain have not been previously reported in the literature. Two examples are shown in Fig. 3 and Fig. 4.
In BRCA, one candidate cluster with several potential novel driver genes is the cluster containing PTEN, HUWE1, CNTNAP2, GRID2, CACNA1B, CYSLTR2, MYH1 depicted in Fig. 3. The genes in the candidate cluster are mutually exclusive (p-value = 0.0084). The genome landscape of this cluster is dominated primarily by mutations in PTEN and HUWE1, and secondarily by homozygous deletions in PTEN and CYSLTR2. The most frequently altered gene in this set is a common driver gene PTEN, a tumor suppressor gene that negatively regulates the AKT/PKB apoptosis pathway [39]. The remaining six genes in the cluster are potential driver candidates. HUWE1 is a part of the Mule multidomain complex of the HECT domain family of E3 ubiquitin ligases responsible for apoptosis suppression, DNA damage repair, and transcriptional regulation [21]. CNT-NAP2 is a neurexin protein with functions in cell-to-cell adhesion and epidermal growth factor and was found to be hypomethylated in breast cancer cell lines [37]. Hypomethylation and the association with epidermal growth factors, coupled with a large number of amplifications in the alteration landscape of CNTNAP2 suggest potential oncogenic functions of the gene. GRID2 is an ionotropic glutamate receptor that is frequently deleted in lymphomas [35]. CACNA1B codes for a N-type calcium channel which is responsible for calcium influx. Defects in the calcium influx channel can lead to alteration in the apoptosis, proliferation, migration and invasion pathways of breast cancer [1]. CYSLTR2 is a proinflammatory cysteinyl leukotriene receptor that plays a role in cancer cell differentiation and is associated with breast cancer survival rates [25]. MYH1 is a myosin heavy chain protein that plays a role in cell signaling and pro-apotosis pathways.
Perhaps more important than the propensity of each individual gene to be a driver is the collective interaction pattern of the seven genes in the cluster in a cancer pathway. From Fig. 3, it is clear that the each gene in the cluster interacts with each other in a tightly-connected community with no gene more than three nodes away when plotted in the network, using the cBioPortal visualization tool. The seven genes in the cluster PTEN, HUWE1, CNTNAP2, GRID2, CACNA1B, CYSLTR2, MYH1 are strong candidates to define a novel driver pathway. This conclusion is reinforced by the presence of high impact common drivers (TP53, MYC, AKT, and PIK3R1) which define several important cancer pathways such as apoptosis, DNA repair, and cell cycle arrest [40,42].
We also examined a cluster containing potential cancer drivers relevant for GBM. In GBM, we found a cluster of size 10 with four known drivers and many potential drivers. The cluster includes GLI1, WNT2, BRAF, PLCG1, FAS, CREBBP, BRCA2, GLI2, PIK3R5, VAMP3 (see Fig. 4). This large cluster has a p-value of 0.0901 in terms of mutual exclusivity, which is actually low as compared to other GBM clusters. The cluster also contains several important driver genes such as WNT2, BRAF, BRCA2 and CREBBP which encompass pathways such as sonic hedgehog signaling, cell fate determination, cell growth and apoptosis, checkpoint activation, and DNA repair. Additionally, six out of the ten members are within the same compact network community (GLI1, PLCG1, FAS, CREBBP, BRCA2, PIK3R5). Of these six genes, GLI1 and GLI2 are hedgehog signaling genes that are common and first isolated in glioblastoma. These genes are responsible for cell differentiation and stem cell self-renewal [10]. PLCG1 is involved in intracellular transduction of receptor-mediated tyrosine kinase activators, and it has been classified as a biomarker in GBM [36]. FAS is a cell surface receptor that mediates apoptosis. FAS is known as a histological hallmark of GBM, affecting both apoptosis and necrosis factors [17]. Finally, PIK3R5 is a subunit of phosphatidylinositol 3-kinases who together have important effects on cell growth, proliferation, differentiation, motility, survival and intracellular trafficking.

Discussion and Conclusion
We described a novel method, termed C 3 , which has the potential to precisely and efficiently identify clusters of gene modules with mutually exclusive mutation patterns. The C 3 algorithm uses large-scale cancer genomics datasets which are pre-processed to yield parameters governing novel constrained correlation clustering techniques. The optimization criteria used in clustering include patterns of mutual exclusivity of mutations, patient sample coverage, and network driver concentration. There are several major advancements of our method when compared to previously known approaches. Unlike methods that use randomized approaches without the guarantee that multiple runs of the methods on the same data will produce compatible results (such as CoMEt), C 3 is consistent. Also, C 3 has a complexity that does not depend on the chosen cluster sizes, and is hence much more appropriate for large cluster problems than other methods. Furthermore, it partitions the gene set and hence creates clusters covering all genes used in the analysis. This is to be contrasted with the results produced by other methods that tend to identify only a small number of modules with limited number of genes.
None of the previous methods were able to identify clusters utilizing different sources of information via a weighting mechanism. This is important because it gives us flexibility to focus more on certain aspects based on the analysis. For example, we can focus more on mutual exclusivity instead of coverage to identify clusters specific to a group of samples which may facilitate the discovery of subtype-specific modules.
By addressing the above challenges, we believe our new method C 3 represents a unique tool to efficiently and reliably identify mutation patterns and driver pathways in large-scale cancer genomics studies.

Supplement -Theoretical Performance Guarantees and Synthetic Data Evaluations
Suppose that we have a weighted correlation-clustering instance, where each edge e incurs a cost w − e if it is placed within a cluster and incurs a cost w + e if it is placed between clusters. We also assume that each cluster may contain at most K + 1 vertices, where K is a fixed constant. If no constraint on the size of the clusters is desired, one can simply set K = |V |.
We wish to find an approximation of the minimum-cost clustering. As already mentioned in the Approach Section, the C 3 method may not work for fully general weights, and, due to the integrality gap in the general case, no algorithm based on LP rounding can work in the fully general weight case. We therefore require the following weight constraints: • w + e ≤ 1 for every edge e, and • w + e + w − e ≥ 1 for every edge e.
Note that these constraints are satisfied if the weights obey the probability constraints w + e + w − e = 1 for all e, but they also apply to much more general choices of weights, which is relevant for the clustering problem at hand.
Recall that we adopted the convention that x uu = 0 for all u.We also observe that one can think of the quantity x uv as a "distance" between u and v. Before starting the analysis of our algorithm, we record an observation based on the triangle inequality: Observation 1. Let x be a feasible LP solution, and let wz be an edge. For any vertex u, we have Given any feasible solution x to the linear program, we use the size-bounded clustering C generated via Algorithm 1, which we express in terms of a parameter α whose optimal value will be determined later in the proof. The idea behind the rounding is to use classical rounding to obtain a solution to the non-size-bounded problem, and if overlarge clusters are produced, to partition them into small-sized clusters.
The standard rounding algorithm comes with a guarantee that cost(C) ≤ 6 cost(x). Here, we cannot ask for such a guarantee: it is possible, for example, that cost(x) = 0 and cost(C) > 0, if all input edges are positive but there are too many vertices to fit into a single cluster. Instead, we prove our approximation ratio by bounding cost(C) in terms of a different lower bound on the cost of an optimal size-bounded clustering.
Let C 1 = cost(x). Choose Y ⊆ E(G) to minimize e∈Y w + e subject to the constraint that each vertex v has at most K incident edges which do not lie in Y . Let C 2 = e∈Y w + e . We think of C 2 as a lower bound on positive weight that must be cut in order to produce a clustering with all clusters of size at most K + 1, since the edges contained within clusters form a subgraph of maximum degree at most K. Thus, the cost of an optimal size-bounded clustering is at least max{C 1 , C 2 }. In particular, an optimal clustering costs at least 7 9 C 1 + 2 9 C 2 . We will use the following lower bound on C 2 : For v ∈ V (G), the excess weight at v, written e(v), is defined by Proof. Choose Y as described above. We have where the second equality follows from the fact that each edge of Y is counted once at each of its endpoints.
We use the "charging" idea in rounding, where one can make "charges" both against the individual contribution of each edge, as well as a global "bank" that will be paid for using C 2 .
Case 1: A singleton cluster {u} is output. The total cluster cost when outputting a singleton cluster {u} is v∈S−u w + uv , while the total LP cost accrued by edges incident to u is v∈S−u w + uv x uv . If the singleton {u} is output, then we have since either the average in step 3 was "too high", or else T was empty. Now for each v ∈ T we have x uv ≤ α.
For such x uv we have 1 − x uv ≥ x uv , since α < 1/2. This yields the following lower bound on the LP cost of uv: where the last inequality uses the bound w + uv + w − uv ≥ 1. Thus, each edge uv has LP cost at least x uv , and so the edges joining u and T have total LP cost at least α |T | /2. Each such edge uv incurs cluster-cost w + uv , which is at most 1. Thus, charging (2/α)x uv to each edge uv for v ∈ T gives enough charge to pay for the cluster-cost of edges with v ∈ T , while each edge is charged at most 2/α times its LP cost.
For v ∈ S − (T ∪ {u}), we have x uv > α, so each edge uv incurs LP cost at least αw + uv and incurs cluster cost at most w + uv . Thus, charging such each edge (1/α)w + uv x uv pays for the cluster cost of these edges. Case 2: A nonsingleton cluster {u} ∪ T is output. We first consider edges inside {u} ∪ T , then we consider edges joining {u} ∪ T with S − ({u} ∪ T ). Finally, we consider edges joining T i with T j for i = j.
Edges within {u} ∪ T . Suppose vz is an edge contained in {u} ∪ T , so that vz incurs cluster-cost w − vz and LP-cost at least w − vz (1 − x vz ). By the definition of T , we have x uv , v uz ≤ α. Hence, by Observation 1, we have to these edges pays for their cluster-cost. Edges joining {u}∪T with S −({u}∪T ). Let z be a vertex outside {u}∪T . A cross-edge for z is an edge from z to {u} ∪ T . We show that the cross-edges for z have total cluster-cost that is at most max{ 1 1−2α , 2 α } times their total LP-cost. Note that whenever vz is a cross-edge, we have x uv ≤ α, by the definition of T . Each cross-edge incurs vz cluster-cost w + vz and LP-cost w + vz x vz + w − vz (1 − x vz ). If in fact x uz ≥ 1−α, then Observation 1 yields x vz ≥ x uz −x uv ≥ 1−2α for every cross-edge v, yielding LP-cost at least 1 1−2α w + vz ; thus, such edges have cluster-cost at most 1 1−2α times their LP cost. It remains to handle the case x uz ∈ (α, 1 − α). Here, we seek a lower bound on the total LP-cost of the cross-edges for z. Note that the total cluster cost of these edges is v∈{u}∪T w + vz , which is at most |T | since each w + vz ≤ 1. By Observation 1, we have x vz ≥ x uz − x uv and 1 − x vz ≥ 1 − x uz − x uv for each edge vz. It follows that the total LP-weight of the cross-edges for z is at least Using the bounds v∈{u}∪T x uv ≤ α(|T | + 1)/2 and x uv ≤ α for v ∈ T , we bound the subtracted sum as follows: Note that in the last inequality we rely on the fact that w + vz + w − vz ≥ 1. It follows that the total LP cost is at least a linear function in x uz . We consider the behavior of this function on the interval (α, 1 − α). When x uz = α, the lower bound simplifies to which is at least α(|T | + 1)/2 since α < 1/2 implies (1 − α)w − vz ≥ αw − vz . When x uz = 1 − α, the lower bound simplifies to which is again at least α(|T | + 1)/2 since (1 − α)w + vz ≥ αw + vz . Thus, when x uz ∈ (α, 1 − α), we conclude that the total cluster-cost of the cross-edges for z is at most 2/α times the total LP-cost of those edges. In all cases, charging max{ 1 1−2α , 2 α } times the LP-cost of the edges pays for their cluster-cost.
Edges joining T i with T j for i = j. These are the trickiest edges to deal with. The problem is that these edges may have low LP-cost and high cluster-cost, so we cannot just charge these edges to pay for their cluster cost. Instead, we will charge edges inside the T i , and we will charge against the bank.
Let v ∈ T i . A cross-edge for v is an edge vw with w ∈ T j for j = i. Note that every edge joining T i , T j with i = j has at least one endpoint that does not lie in T p ; as such, we will assume that v / ∈ T p . Define quantities A v and B v as follows: {w + e : e is a crossing edge for v}, Since |T i | = K + 1, we have where the second inequality follows from w + vz + w − vz ≥ 1 and the final equality follows from the fact that |T i | = K + 1 for i < p. On the other hand, since each w + e ≤ 1, we have We charge the quantity e(v) to the bank, and we charge w − vz to each edge vz with z ∈ T i . This yields total charge at least A v , which is the total cluster-cost of the cross-edges for v. Note that for each edge vz charged this way, Observation 1 yields 1 − x vz ≥ 1 − x uv − x vz ≥ 1 − 2α. Hence, v charges the edge vz at most 1 1−2α times its LP-cost. Observe that each edge in T is only charged this way at its endpoints. (Furthermore, edges whose endpoints lie in different T i , T j with both i, j < p are actually paid for twice.) In total, we have paid for all the cluster-costs by making the following charges: • Edges vz within each cluster T i were charged at most 1 1−2α w − vz (1 − x vz ) to pay for themselves plus at most 2 1−2α to pay for edges joining T i to T j for i = j. Thus, their total charge is at most 3 1−2α times their total LP cost.
• The bank was charged v∈V e(v). Thus, by Lemma 1, the total charge to the bank is at most 2C 2 .
It follows that the total cost is minimized when 3 1−2α = 2 α , which yields α = 2/7. For this choice of α, the total charge is at most 7C 1 + 2C 2 , which is at most 9 times the lower bound of 7 9 C 1 + 2 9 C 2 ." Given that the theoretical guarantees only establish approximation results, it is appropriate to compare the solution of the proposed relaxation with the solution of the original ILP on synthetic data. For this purpose, we created a simple graph comprising 35 vertices grouped into five clusters of size 6 and one cluster of size 5: in order to ensure this cluster separation, we selected the weights as follows. For u, v ∈ V (G), we set w + uv = γ, if u and v are in the same cluster (21) w + uv = 1 − γ, otherwise, and w − uv = 1 − w + uv. We ran the ILP, the C 3 algorithm, and the clu for γ ∈ {0.6, 0.7, 0.8, 0.9, 0.99}; all three algorithms recovered the correct clusters without making errors. We subsequently modified the clusters by randomly reversing the positive and negative weights of up to 20 edges. All three algorithms recovered the correct clusters even in this scenario. Since large weight perturbations may change the clusters in a way that the ground truth becomes undetectable, for larger scale perturbation we decided to compare the solution of the C 3 algorithm to that of the ILP method. We again ran the ILP and C 3 algorithms on a graph with 35 vertices, where the positive weights, w + e , were chosen randomly according to a multinomial distribution; we also considered a number of choices for the parameters of the distributions. The negative weights were chosen according to the formula w − e = 1 − w + e . For each set of weights, we compared the value of the objective function for the results obtained using the two algorithms. We observed that in all cases the value of objective functions, denoted by f , satisfied which shows that in practice, the C 3 algorithm performed much better than suggested by the theoretical analysis. Furthermore, the clusters generated by the ILP and the C 3 algorithm showed more than 90% overlap in terms of participating nodes.