Efficient algorithms to discover alterations with complementary functional association in cancer
Rebecca Sarto Basso Dorit S. Hochbaum Fabio Vandin







Efficient algorithms to discover alterations with complementary functional association in cancer 





Rebecca Sarto Basso1,
Dorit S. Hochbaum1,
Fabio Vandin2,3,4*


1 Department of Industrial Engineering and Operations Research, University of California at Berkeley, CA, USA.


2 Department of Information Engineering, University of Padova, Padova, Italy.


3 Department of Computer Science, Brown University, Providence, RI, USA.


4 Department of Mathematics and Computer Science, University of Southern Denmark, Odense, Denmark.




























* Corresponding author: fabio.vandin@unipd.it






§ ABSTRACT


Recent large cancer studies have measured somatic alterations in an unprecedented number of tumours. These large datasets allow the identification of cancer-related sets of genetic alterations by identifying relevant combinatorial patterns. Among such patterns, mutual exclusivity has been employed by several recent methods that have shown its effectivenes in characterizing  gene sets associated to cancer.  Mutual exclusivity arises because of the complementarity, at the functional level, of alterations in genes which are part of a group (e.g., a pathway) performing a given function. The availability of quantitative target profiles, from genetic perturbations or from clinical phenotypes, provides additional information that can be leveraged to improve the identification of cancer related gene sets by discovering groups with complementary functional associations with such targets.

In this work we study the problem of finding groups of mutually exclusive alterations associated with a quantitative (functional) target. We propose a combinatorial formulation for the problem, and prove that the associated computation problem is computationally hard. We design two algorithms to solve the problem and implement them in our tool . We provide analytic evidence of the effectiveness of  in finding high-quality solutions and show experimentally that  finds sets of alterations significantly associated with functional targets in a variety of scenarios. In particular, we show that our algorithms find sets which are better than the ones obtained by the state-of-the-art method, even when sets are evaluated using the statistical score employed by the latter. In addition, our algorithms are much faster than the state-of-the-art, allowing the analysis of large datasets of thousands of target profiles from cancer cell lines. We show that on one such dataset from project Achilles our methods identify several significant gene sets with complementary functional associations with targets. Software available at: <https://github.com/VandinLab/UNCOVER>.




§ INTRODUCTION

Recent advances in sequencing technologies now allow to collect genome-wide measurements in large cohorts of cancer patients (e.g., <cit.>). In particular, they allow the measurement of the entire complement of somatic (i.e., appearing during the lifetime of an individual) alterations in all samples from large tumour cohorts. The study of such alterations has lead to an unprecedented  improvement in our understanding of how tumours arise and progress <cit.>. One of the main remaining challenges is the interpretation of such alterations, in particular identifying alterations with functional impact or with relevance to therapy <cit.>.

Several computational and statistical methods have been recently designed to identify driver alterations, associated to the disease, and to distinguish them from random, passenger alterations not related with the disease <cit.>. The identification of genes associated with cancer is complicated by the extensive intertumour heterogeneity <cit.>, with large (100-1000's) and different collections of alterations being present in tumours from different patients and no two tumours having the same collection of alterations <cit.>. Two main reasons for such heterogeneity are that [i)]

  * most mutations are passenger, random mutations, and,  more importantly, 

  * driver alterations target cancer pathways, 

groups of interacting genes that perform given functions in the cell and whose alteration is required to develop the disease. Several methods have been designed to identify cancer genes using a-priori defined pathways <cit.> or interaction information in the form of large interaction networks <cit.>.

Recently several methods (see Section <ref>) for the de novo discovery of mutated cancer pathways have leveraged the mutual exclusivity of alterations in cancer pathways. Mutual exclusivity of alterations, with sets of genes displaying at most one alteration for each patient, has been observed in various cancer types <cit.>.  The mutual exclusivity property is due to the complementarity of genes in the same pathway, with alterations in different members of a pathway resulting in a similar impact at the functional level. Mutual exclusivity has been successfully used to identify cancer pathways in large cancer cohorts <cit.>.

An additional source of information that can be used to identify genes with complementary functions are quantitative measures for each samples such as: functional profiles, obtained for example by genomic or chemical perturbations <cit.>;  clinical data describing, obtained for example by (quantitative) indicators of response to therapy; activation measurements for genes or sets of genes, as obtained for example by single sample scores of Gene Set Enrichment Analysis <cit.>. The employment of such quantitative measurements is crucial to identify meaningful complementary alterations since one can expect mutual exclusivity to reflect in functional properties (of altered samples) that are specific to the altered samples.



 §.§ Related work



Several recent methods have used mutual exclusivity signals to identify sets of genes important for cancer <cit.>. RME <cit.> identifies mutually exclusive sets using a score derived from information theory.  Dendrix <cit.> defines a combinatorial gene set score and uses a Markov Chain Monte Carlo (MCMC) approach for identifying mutually exclusive gene sets altered in a large fraction of the patients. Multi-Dendrix <cit.> extends the score of Dendrix to multiple sets and uses an integer linear program (ILP) based algorithm to simultaneously find multiple sets with mutually exclusive alterations. CoMET <cit.> uses a generalization of Fisher exact test to higher dimensional contingency tables to define a score to characterize mutually exclusive gene sets altered in relatively low fractions of the samples. WExT <cit.> generalizes the test from CoMET to incorporate individual gene  weights (probabilities) for each alteration in each sample. WeSME <cit.> introduces a test that incorporates the alteration rates of patients and genes and uses a fast permutation approach to assess the statistical significance of the sets. TiMEx <cit.> assumes a generative model for alterations and defines a test to assess the null hypothesis that mutual exclusivity of a gene set is due to the interplay between waiting times to alterations and the time at which the tumor is sequenced. MEMo <cit.> and the method from <cit.> employ mutual exclusivity to find gene sets, but use an interaction network to limit the candidate gene sets. The method by <cit.> and PathTIMEx <cit.> introduce an additional dimension to the characterization of inter-tumor heterogeneity, by reconstructing the order in which mutually exclusive gene sets are mutated. None of these methods take quantitative targets into account in the discovery of significant gene sets and sets showing high mutual exclusivity may not be associated with target profiles (Fig. <ref>).




<cit.> recently developed the repeated evaluation of variables conditional entropy and redundancy (REVEALER) method, to identify mutually exclusive sets of alterations associated with functional phenotypes. REVEALER uses as objective function (to score a set of alterations) a re-scaled mutual information metric called information coefficient (IC).  REVEALER employs a greedy strategy, computing at each iteration the conditional mutual information (CIC) of the target profile and each feature, conditioned on the current solution. REVEALER can be used to find sets of mutually exclusive alterations starting either from a user-defined seed for the solution or from scratch, and <cit.> shows that REVEALER finds sets of meaningful cancer-related alterations.



 §.§ Our contribution


In this paper we study the problem of finding sets of alterations with complementary functional associations using alteration data and a quantitative (functional) target measure from a collection of cancer samples. Our contributions in this regard are fourfold. First,
we provide a rigorous combinatorial formulation for the problem of finding groups of mutually exclusive alterations associated with a quantitative target and prove that the associated computational problem is NP-hard. Second, we develop two efficient algorithms, a greedy algorithm and an ILP-based algorithm to identify the set of k genes with the highest association with a target; our algorithms are implemented in our method fUNctional Complementarity of alteratiOns discoVERy (). Third, we show that our algorithms identify highly significant sets of genes in various scenarios; in particular, we compare  with REVEALER on the same datasets
used in <cit.>, showing that  identifies solutions of  higher quality than REVEALER while being on average two order of magnitudes faster than REVEALER. Interestingly, the solutions obtained by  are better than the ones obtained by REVEALER even when evaluated using the objective function (IC score) optimized by REVEALER.
Fourth, we show that the efficieny of  enables the analysis of a large dataset from Project Achilles with thousands of target measurements and tens of thousands of alterations. On such dataset  identifies identifies several statistically significant associations between target values and mutually exclusive alterations in genes sets
sets,  with an enrichment in well-known cancer genes and in known cancer pathways.



§ MATERIALS AND METHODS


This section describes the problem we study and the algorithms we designed to solve it, that are implemented in our tool . We also describe the data and computational environment for our experimental evaluation.



 §.§ : Functional complementarity of alterations discovery


The workflow of our algorithm  is presented in Fig. <ref>. takes in input information regarding 

  * the alterations measured in a number of samples (e.g., patients or cell lines), and

  * the value of the target measure for each patient.

 then identifies the set of mutually exclusive alterations with the highest association to the target, and employs a permutation test to assess the significance of the association. Details regarding the computational problem and the algorithms used by  are described in the following sections. The implementation of  is available at <https://github.com/VandinLab/UNCOVER>.


 


  §.§.§ Computational problem



Let J={j_1,…,j_m} be the set of samples and let G = {g_1,…, g_n} be the set of genes for which we have measured alterations in J. We are also given a target profile, that is for each sample j ∈ J we have a target value w_j ∈ℝ which quantitatively measures a functional phenotype (e.g., pathway activation, drug response, etc.). For each sample j∈ J we also have information on whether each g ∈ G is altered or not in j. 
Let A_g be the set of patients in which gene g ∈ G is mutated. We say that a patient j ∈ J is covered by gene g∈ G if j ∈ A_g i.e. if gene g is mutated in sample j. Given a set of genes S ⊂ G, we say that sample j ∈ J is covered by S if j ∈∪_g ∈ S A_g.

The goal is to identify a set S of at most k genes, corresponding to k subsets S_1, S_2,… S_k where for each subset S_i we have that S_i ⊆ J, such that the sum of the weights of the elements covered by S is maximized. We also penalize overlaps between sets when an element is covered more than once by S by assigning a penalty p_j for each of the additional times j is covered by S. As penalty we use the positive average of the normalized target values if the original weight of the element was positive. If the original weight of the element was negative we assign a penalty equal to its weight.


Let c_S(j) be the number of sets in S_1, …, S_k that cover element j ∈ J.

Therefore for a set S of genes, we define its weight W(S) as:


    W(S) = ∑_j ∈∪_s ∈ SA_s w_j - ∑_j ∈∪_s ∈ SA_s(c_S(j)-1) p_j


The  problem
Given a set J of samples, sets A_g_1,…,A_g_n describing alterations of genes G={g_1,…,g_n} in the set J, weights w_j and penalties p_j>0 for each sample j ∈ J find the S of ≤ k elements maximizing W(S).


The following results defines the computational hardness of the problem above.



The   problem is NP-hard.


The proof is by reduction from the Maximum Weight Submatrix Problem (MWSP) defined and proved to be NP-hard in <cit.>. The MSWP takes as input an m×  n binary matrix A and an integer k > 0 and requires to find the m × k column sub-matrix M̂ of A that maximizes the objective function |Γ(M)| - ω(M), where Γ(M) is the set of rows with at least one 1 in columns of M and ω(M)=∑_g ∈ M |Γ({g})| - |Γ(M)|.

Given an instance of Maximum Weight Submatrix Problem, we define an instance of the  as follow: the set of samples J corresponds to the rows of A, the set of genes G corresponds to the columns of A, and the set S_g of samples covered by gene g∈ G is the subset of the rows in which g has a 1. By setting w_j = 1 and p_j = 1 for all j ∈ J, we have that the objective function of MWSP corresponds to the weight W(S) for the   therefore the optimal solution of the  corresponds to the optimal solution of MWSP.





  §.§.§ ILP formulation



In this subsection we provide an ILP formulation for the  problem. Let x_i be a binary variable equal to 1 if set i ∈ G is selected and x_i=0 otherwise. Let z_j be a binary variable equal to 1 if element j is covered and z_j=0 otherwise. Let  y_j denote the number of sets in the  solution covering element j. Finally, let w_j be the weight of element j and  p_j be the penalty for element j

In our ILP formulation, the following constraints need to be satisfied by a valid solution:

    
  * the total number of sets in the solution is at most k: ∑_i x_i≤ k
    
  * for each element j ∈ J we have:  y_j=∑_i: j ∈ S_i x_i
    
  * for each element j ∈ J, if j is covered by the current solution then the number of times j is covered in the solution is at least 1:  y_j≥ z_j
    
  * for each element j ∈ J, if j is covered by at least one element in the current solution then j is covered: z_j≥ y_j/k.


With the variables defined above, the score for a given solution is 

    z(q) = max∑_j=1^m (w_j+p_j) z_j-∑_j=1^m p_j y_j.


z(q) constitutes the objective function of our ILP formulation.



  §.§.§ Greedy algorithm



Since solving ILPs can be impractical for very large datasets, we also design a k-stage greedy algorithm to solve the  problem. During each stage the algorithm picks 1 set A_i to be part of the solution; this is done by first computing the total weight of each subset which is defined as the sum of the weights of its elements W(A_i)=∑_j ∈ A_i w_j. Then the algorithm finds the subset A_i of maximum positive weight and adds it to the solution. It may be that at some stage ℓ no additional set of positive weight can be selected, in this case, the solution obtained after stage ℓ-1 will be our output.  
At the end of the iteration the weight of every element j that belonged to the chosen set A_i is set to the negative of penalty p_j, in order to penalize future selections of the same elements. The greedy algorithm is described in Algorithm <ref>.



We note that our greedy algorithm is analogous to the greedy algorithm for the Maximum k-Coverage problem <cit.> with the difference that rather than eliminating the elements already selected we change their weight to a penalty. Also, assuming it is acceptable to return less than k sets, we only pick a set if it has a positive weight. The running time of the algorithm is O(kmn) where m = number of samples and n = number of alterations.

While the greedy algorithm may not return the optimal solution,
we prove that it provides guarantees on the weight of the solution it provides.


Let S^* the optimal solution of the  and Ŝ be the solution returned by the greedy algorithm. Then W(Ŝ) ≥ W(S^*)/k.



Note that the weight of subsets in the optimal solution W(S^*) can only be lower compared to the original weight of the subsets, since the only weight update operation performed is to substitute positive weights of elements selected with a negative penalty. 

The first subset Ŝ_̂1̂  selected by our algorithm is the set of maximum weight out of all subsets and therefore W(Ŝ_̂1̂) ≥ W(S_ℓ^*) for ℓ=1..k. By the pigeonhole principle, one of these subsets in the optimal solution must cover at least W(S^*)/k worth of elements. Thus W(Ŝ_̂1̂) ≥ W(S^*)/k.
Therefore the first subset selected by the algorithm already gives a 1/k approximation of the optimal solution. In subsequent iterations of the algorithm we only pick additional sets if they have a positive weight so our approximation ratio can only improve. 


We also prove that the bound above is tight


There are instances of the  such that W(Ŝ) = W(S^*)/k.


The proof is in the Supplementary Material.

While the proposition above is based on an extreme example, our experimental analysis shows that in practice the greedy algorithm works well and often identifies the optimal solution. We therefore analyze the greedy algorithm under a generative model in which there is a set H of k genes with mutually exclusive alterations associated with the target, while each genes g  ∈ G ∖ H is mutated in sample j with probability p_g independently of all other events. We also assume that the weights w_j are such that ∑_j∈ J w_j = 0 and for each j: |w_j ≤ 1|. (In practice this is achieved by normalizing the target values before running the algorithm, by subtracting to each w_j the average value ∑_j∈ J w_j/m and then dividing the result by the maximum of the absolute values of the transformed w_j's.) Note that this last condition implies that |p_j| ≤ 1 for all j. We also assume that for genes in H the following assumptions hold:

    
  * the set H has an association with the target, i.e.: [W(H)] ≥m/c' for a constant c'≥ 1.
    
  * each gene of H contributes to the weight of H, i.e. for each S ⊂ H and each g ∈ H ∖ S we have [W(S ∪{g})] - [W(S)]≥W(H)/k c^” for a constant c^”≥ 1.


The following shows that, if enough samples from the generative model are considered, the greedy algorithm finds the set H associated with the target with high probability.


If m ∈k^2 ln (n/δ) samples from the generative model above are provided to the greedy algorithm, then the solution of the greedy algorithm is H with probability ≥δ.


The proof is in the Supplementary Material.



  §.§.§ Statistical significance



To assess the significance of the solution reported by our algorithms we use a permutation test in which the dependencies among alterations in various genes are maintained, while the association of alterations and the target is removed. In particular, a permuted dataset under the null distribution is obtained as follows: the sets A_g, g∈𝒢 are the same as observed in the data; the values of the target are randomly permuted across the samples. 

To  estimate the p-value for the solutions obtained by our methods we used the following standard procedure: [1)] 
  * we run an algorithm (ILP or greedy) on the real data 𝒟, obtaining a solution with objective function o_𝒟; 
  * we generate N permuted datasets as described above; 
  * we run the same algorithm on each permuted dataset; 
  * the p-value is the given by (e+1)/(N+1), where e is the number of permuted datasets in which our algorithm found a solution with objective function ≥ o_𝒟.



 §.§ Data and computational environment





  
Alteration Data. We downloaded data for the Cancer Cell Line Encyclopedia on 25^th September, 2017 from 
<http://www.broadinstitute.org/ccle>. In particular we used the mutation (single nucleotide variants) and copy number aberrations (CNAs) which are derived from the original Cancer Cell Line Encyclopedia (CCLE) mutations and CNA datasets. The file we used is . It consists of a binary (0/1) matrix across 1,030 samples and 48,270 gene alterations in the form of mutations, amplifications and deletions, with a 1 meaning that the alteration is present in a sample, and a 0 otherwise.



  
Target Data. In terms of target values we use the same datasets used by <cit.> to compare the performance of  with REVEALER.
In particular we used the following files available through the Supplementary Material of <cit.>: , which consists of measurements of a β-catenin reporter in 81 cell lines; , which includes NFE2L2  enrichment profiles for 182 lung cell lines; , which contains MEK-inhibitor PD-0325901 sensitivity profile in 493 cancer cell lines from the Broad Novartis CCLE14l; and , which  corresponds to the feature KRAS from a subset of 100 cell lines from the Achilles project dataset. 
In all these cases we considered the same direction of association (positive or negative) between alterations and the target as in <cit.>. Since our algorithm is very efficient we then decided to run it on a large dataset from Project Achilles  (<https://portals.broadinstitute.org/achilles>), that uses genome-scale RNAi and CRISPR-Cas9 perturbations to silence or knockout individual genes. In particular, we use the whole 2.4.2 Achilles dataset () available from the project website. This dataset provides phenotype values for 5711 targets, corresponding to genes silenced by shRNA. The phenotype values correspond to ATARiS <cit.> gene (target) level scores, quantifying the cell viability when the target gene is silenced by shRNA. These scores are provided for 216 cell lines <cit.>, with 205 of them appearing in CCLE.



  
Data Preprocessing. To be consistent with REVEALER we discarded features with high or low frequency, in particular features present in less than 3 samples or more than 50 samples were excluded from our analyses. The only exception was the MEK-inhibitor example, where the high frequency threshold was changed to be 100 since the number of original samples was substantially higher (i.e., 493) for this case.  From the Achilles dataset we excluded targets that have at least one missing value, in particular in this case we exclude 21 of the 5711 sets of target scores.
In all our experiments we normalized the target values before running the algorithm, by subtracting to each weight w_j the average value ∑_j∈ J w_j/m and dividing the result by the standard deviation of the (original) w_j's, in order to have both positive and negative target values.



  
Simulated Data. 
We investigated how effective   is at finding selected alterations in a controlled setting, where the ground truth is known. We generated target values according to a normal distribution with mean 0 and standard deviation 1. We tested dataset with 200, 600, 1000 and 10000 samples. For each dataset we considered the 38002 gene alterations present in CCLE and for each of them we placed alterations in the samples independently of all other events with the same frequencies as they appear in CCLE. To be consistent with the preprocessing done on rel data we filtered alterations to only have alterations with frequencies between 0.1 and 0.25. We also generated a set T of 5 features to have an association with the target values. This association was varied throughout the experiments to cover different percentages of positive and negative targets. In particular we generated the selected features to cover 100%, 80%, 60%, 40% of the positive target values and 5%, 10%, 15%, 20% of the negative target values respectively, chosing random subsets of samples with positive or negative target values. We will refer to the parameter indicating the percentage of samples with positive target values selected as P and to the parameter for the percentage of samples with negative target values selected as N. We divided the number of targets covered by each of the 5 mutations equally.



  
Computing Environment and Solver Configuration. To describe and solve an ILP we used AMPL 20150516 and CPLEX 12.6.3.  All parameters in CPLEX were left at their default values. We implemented our greedy algorithm in Python 3.6.1. We run our experiments (with the exception of experiments conducted on simulated data) on a MacBook Air with 1.7 GHz Intel Core i7 processor, 8 GB RAM and 500 GB of local storage. In order to make a time comparison with REVEALER we also run the R code provided in <cit.> on the same machine, using R 3.3.3. All the parameters were left at their given values except for the number of permutations used to calculate their p-value, which we changed in order to compare the running time of the methods excluding the time needed to compute p-values.
Experiments on simulated data were conducted on local nodes of a computing cluster. Each node had the following configuration: four 2.27 GHz CPUs, 5.71 GB RAM and 241 GB of storage. 



§ RESULTS AND DISCUSSION


We tested  on a number of cancer datasets in order 
to compare its results with state-of-the-art algorithms and to test whether  allows the 
analysis of large datasets. In particular we used four datasets described in <cit.> to compare
 with REVEALER, and then performed a large scale analysis using targets from the Achilles project dataset and alterations from the Cancer Cell Line Encyclopedia (CCLE). 



 §.§ Comparison with REVEALER



We run the greedy algorithm and the ILP from  on the same four datasets considered by the REVEALER publication <cit.>. We used the same values of k used in <cit.>, that is k=3 for all the datasets, except from the KRAS dataset where k=4 was used. For each dataset we recorded the solution reported by the greedy algorithm, the solution reported by the ILP, the value of the objective functions for such solutions and the running time to obtain such solutions. For ILP solutions, we also performed the permutation test (see Materials and methods) to compute a p-value using 1000 permutations. The results are reported in Table <ref>, in which we also show the results from REVEALER (without initial seeds). Fig. <ref> shows alteration matrices and the association with the target for the solutions identified by .



We can see that the greedy algorithm identifies the same solution of the ILP based algorithm in three out of four cases, and that the runtime of the ILP and the runtime of greedy algorithm are comparable and very low (<40 seconds) in all cases. In contrast, the running time of REVEALER is much higher (>1000 seconds in most cases).  (We included all preprocessing in the reported runtimes in table 1 to ensure a fair comparison with REVEALER; not including preprocessing our running times are all under 10 seconds.)  Comparing the alteration matrices of the solutions by  and the ones of solutions by REVEALER (Supplementary Fig. <ref>) we note that alterations in solutions by  tend to have higher mutual exclusivity and to be more concentrated in high weight samples than alterations in solutions by REVEALER. As expected, the value of the objective function we use is much lower for solutions from REVEALER than for solutions from our algorithm.



We then compared the solutions obtained by our algorihtms with the solutions from REVEALER in terms of the information coefficient (IC), that is the target association score used in <cit.> as a quality of the solution. Surprisingly, in two out of four datasets , which does not consider the IC score, identifies solutions with IC score higher (by at least 5%) than the solutions reported by REVEALER. For the other two cases, in one dataset the IC score is very similar (0.50 vs 0.49) while in the other case the IC score by REVEALER is higher (0.7 vs 0.67) but the solution reported by REVEALER differs from the solution reported by  by 1 gene only.
Interestingly, the latter is the only case where the solution from the ILP has a p-value >0.1 (p<0.03 in all other cases), and therefore the solutions (by our methods and by REVEALER) for such dataset may be, at least in part, due to random fluctuations of the data.

In most cases the solutions by  and by REVEALER are very similar, with cancer relevant genes identified by both methods. For NFE2L2 activation, both methods identify KEAP1, a repressor of NFE2L2 activation <cit.>. For MEK-inhibitor, both methods find BRAF, KRAS, and NRAS, three well knwon oncogenic activators of the MAPK signaling pathway, which contains MEK as well. For KRAS essentiality, both methods report mutations in KRAS in the solution. For β-catenin activation, both methods identify CTNNB1 mutations and APC mutations, that is known to be associated to β-catenin activation <cit.>.
These results show that  identifies relevant biological solutions that are better than the ones identified by REVEALER when evaluated using our objective function and also when evaluated according to the objective function of REVEALER with a running time that is on average two orders of magnitude smaller than required by REVEALER.



 §.§ Results on simulated data


    For each combination we generated 10 simulated dateset as described in Materials and methods. Each dataset contains a planted set of 5 alterations associated with the target. We used both the greedy algorithm and the ILP from   with k=5 to attempt to find the 5 correct alteration, and evaluated our algorithms both in terms of fraction of the correct (i.e., planted) solution reported and running time.

    
 As shown in Fig. <ref>, the greedy algorithm is faster than the ILP for all datasets, and the difference in running time increases as the number m of samples increases, with the runtime of the greedy algorithm being almost two order of magnitude smalle than the runtime of the ILP for m=1000 samples. In addition, for a fixed number of samples and alterations, the running time of the greedy algorithm is constant, that is it does not depend on the properties of the planted solution, while the running time of the ILP varies greatly depending on these parameters. For m=10,000 samples the running time of the ILP becomes extremely high, so we restricted to consider only two sets of paramters (p-n=0.95 and p-n=0.2). In this case the ILP took between 44 minutes and 7 hours to complete, while the greedy algorithm terminates in 5 minutes.
 
 In terms of the quality of the solutions found, as expected the ILP outperforms the greedy (Fig. <ref>) but the difference among the two tends to disappear when the number of samples is higher. In addition, since the ILP finds the optimal solution, we can see that for a limited number of samples we may not reliably identify the planted solution with 200 samples unless the planted solution appears almost only in positive targets and in almost all of them (p-n=0.95), while for m=1000 we can reliably identify the planted solution using both the ILP and the greedy algorithm even when the association with the target is weaker (p-n=0.6). When m=10,000, both the ILP and the greedy algorithm perform well in terms of the quality of the solution: the ILP finds the correct alterations on every experiment and the greedy identifies the whole planted solution in all experiments but one for p-n=0.2, for which it still reports a solution containing 4 genes in the planted solution.
 
 These results show that for a large number of samples the greedy algorithm reliably identifies sets of alterations associated with the target, as predicted by our theoretical analysis, and is much faster than the ILP. For smaller sample size the ILP identifies better solutions than the greedy and has a reasonable running time.
 
 
 
 
 

 

 


 §.§ Analysis of Achilles project data




The efficiency of  renders the analysis of a large number of targets, such as the ones available through the Achilles project, possible. After preprocessing the dataset included 5690 functional phenotypes as targets, and for each of these the CCLE provides alteration information for 205 samples and 31137 alterations. In total we have therefore run 10380 instances (i.e., 5690 targets screened for positive and for negative associations) looking for both positive and negative association with target values. Since the number of samples (205) is relatively small, we have run only the ILP from  on the whole Achilles dataset and looked for solutions with k=3 genes.
The runtime of  to find both positive and negative associations, including preprocessing, is 24 hours. Based on the runtime required on the instances reported in <cit.> (see the Comparison with REVEALER section), running REVEALER on this dataset would have required about 5 months of compute time.

To identify statistically significant associations with targets in the Achilles project dataset we used a nested permutation test. We first run the permutation test with 10 permutations on all instances (i.e., on all targets for both positive association and negative association). We then considered all the instances with the lowest p-value (1/11) and performed a permutation test with 100 permutations only for such instances. We the iterated such procedure once more, selecting all the instances with lowest p-value (1/101) and performing a permutation test with 1000 permutations only for such instances.
For positive association we found 60 solutions with p-value < 0.001, and for negative association we found 102 solutions with p-value < 0.001. The solutions with p-value < 0.001 (with 1000 permutations) are reported in Supplementary Table <ref> and  <ref>. See Supplemental Fig. <ref> for some corresponding alteration matrices.

The genes in the solutions by  with p-value 1/1001 are enriched (p=2× 10^-12 by Fisher exact test) for well-known cancer genes, as reported in <cit.>. We also tested whether genes in solutions by  (with p-value 1/1001) are enriched for interactions, by comparing the number of interactions in  <cit.> among genes in such solution with the number of interactions in random sets of genes of the same cardinality. Genes in solutions by  are significantly enriched in interactions (p=7× 10^-3 by permutation test). In addition, the genes in solutions by  are also enriched in genes in well-known pathways: 12 KEGG pathways <cit.> have a significant (corrected p ≤ 0.05) overlap with genes in solutions by  and four of these (endometrial cancer, glioma, hepatocellular carcinoma, EGFR tyrosine kinase inhibitor resistance) are cancer related pathways. In addition, the targets (i.e., genes) with solutions of p-value 1/1001 are enriched (p=10^-3 by permutation test) for interactions in  and for well-known cancer genes as reported in <cit.>. These results show that  enables the identification of groups of well known cancer genes with significant associations to important targets in large datasets of functional target profiles. For example, for target (i.e., silenced gene) TSG101, related to cell growth,  identifies the gene set shown in Figure <ref> as associated to reduced cell viability. ERBB2 is a well known cancer gene and CDH4 is frequently mutated in several cancer types, and both are associated to cell growth.






§ CONCLUSION



In this work we study the problem of identifying sets of mutually exclusive alterations associated with a quantitative target profile. 
We provide a combinatorial formulation for the problem, proving that the corresponding computational problem is NP-hard. We design two efficient algorithms, a greedy algorithm and an ILP-based algorithm, for the identification of sets of mutually exclusive alterations associated with a target profile. We provide a formal analysis for our greedy algorithm, proving that it returns solutions with rigorous guarantees in the worst-case as well under a reasonable genarative model for the data. We implemented our algorithms in our method  , and showed that it finds sets of alterations with a significant association with target profiles in a variety of scenarios. By comparing the results of UNCOVER with the results of REVEALER on four target profiles used in the REVEALER publication <cit.>, we show that  identifies better solutions than REVEALER, even when evaluated using REVEALER objective function. Moreover,  is much faster than REVEALER, allowing the analysis of large datasets such as the dataset from project Achilles, in which  identifies a number of associations between functional target profiles and gene set alterations.
Our tool  (as well as REVELEAR) relies on the assumption that the mutual exclusivity among alterations is due to functional complementarity. Another explanation for mutual exclusivity is the fact that each cancer may comprise different subtypes, with different subtypes being characterized by different alterations <cit.>.  can be used to identify sets of mutually exclusive alterations associated with a specific subtype whenever the subtype information is available, by assigning high weight to samples of the subtype of interest and low weight to samples of the other subtypes.




§ ACKNOWLEDGEMENTS

  This work is supported, in part, by NSF grant IIS-124758 and by the University of Padova grants SID2017 and Algorithms for Inferential Data Mining, funded by the STARS program. This work was done in part while the authors were visiting the Simons Institute for the Theory of Computing, supported by the Simons Foundation. A short version of this work was accepted to RECOMB 2018.
 

















10

Brennan:2013aa
Brennan CW, Verhaak RGW, McKenna A, Campos B, Noushmehr H, Salama SR, et al.
The somatic genomic landscape of glioblastoma.
Cell. 2013;155(2):462–77.
doi:10.1016/j.cell.2013.09.034.

Cancer-Genome-Atlas-Network:2015aa
Cancer Genome Atlas Network.
Comprehensive genomic characterization of head and neck squamous cell
  carcinomas.
Nature. 2015;517(7536):576–82.
doi:10.1038/nature14129.

Cancer-Genome-Atlas-Research-Network:2013aa
Cancer Genome Atlas Research Network.
Comprehensive molecular characterization of clear cell renal cell
  carcinoma.
Nature. 2013;499(7456):43–9.
doi:10.1038/nature12222.

Cancer-Genome-Atlas-Research-Network:2014aa
Cancer Genome Atlas Research Network.
Integrated genomic characterization of papillary thyroid carcinoma.
Cell. 2014;159(3):676–90.
doi:10.1016/j.cell.2014.09.050.

cancer2017integrated
Network CGAR, et al.
Integrated genomic characterization of oesophageal carcinoma.
Nature. 2017;541(7636):169–175.

cancer2017integratedB
Network CGAR, et al.
Integrated Genomic Characterization of Pancreatic Ductal
  Adenocarcinoma.
Cancer cell. 2017;32(2):185.

garraway2013lessons
Garraway LA, Lander ES.
Lessons from the cancer genome.
Cell. 2013;153(1):17–37.

mcgranahan2017clonal
McGranahan N, Swanton C.
Clonal heterogeneity and tumor evolution: past, present, and the
  future.
Cell. 2017;168(4):613–628.

pmid24479672
Raphael BJ, Dobson JR, Oesper L, Vandin F.
Identifying driver mutations in sequenced cancer genomes:
  computational approaches to enable precision medicine.
Genome Med. 2014;6(1):5.

vandin2017computational
Vandin F.
Computational Methods for Characterizing Cancer Mutational
  Heterogeneity.
Frontiers in genetics. 2017;8:83.

vogelstein2013cancer
Vogelstein B, Papadopoulos N, Velculescu VE, Zhou S, Diaz LA, Kinzler KW.
Cancer genome landscapes.
science. 2013;339(6127):1546–1558.

vaske2010inference
Vaske CJ, Benz SC, Sanborn JZ, Earl D, Szeto C, Zhu J, et al.
Inference of patient-specific pathway activities from
  multi-dimensional cancer genomics data using PARADIGM.
Bioinformatics. 2010;26(12):i237–i245.

Leiserson:2015aa
Leiserson MDM, Vandin F, Wu HT, Dobson JR, Eldridge JV, Thomas JL, et al.
Pan-cancer network analysis identifies combinations of rare somatic
  mutations across pathways and protein complexes.
Nat Genet. 2015;47(2):106–14.
doi:10.1038/ng.3168.

creixell2015pathway
Creixell P, Reimand J, Haider S, Wu G, Shibata T, Vazquez M, et al.
Pathway and network analysis of cancer genomes.
Nature methods. 2015;12(7):615.

kandoth2013mutational
Kandoth C, McLellan MD, Vandin F, Ye K, Niu B, Lu C, et al.
Mutational landscape and significance across 12 major cancer types.
Nature. 2013;502(7471):333–339.

hanahan2011hallmarks
Hanahan D, Weinberg RA.
Hallmarks of cancer: the next generation.
cell. 2011;144(5):646–674.

Ciriello:2012ly
Ciriello G, Cerami E, Sander C, Schultz N.
Mutual exclusivity analysis identifies oncogenic network modules.
Genome Res. 2012;22(2):398–406.
doi:10.1101/gr.125567.111.

leiserson2015comet
Leiserson MD, Wu HT, Vandin F, Raphael BJ.
CoMEt: a statistical approach to identify combinations of mutually
  exclusive alterations in cancer.
Genome biology. 2015;16(1):160.

pmid25984343
Cowley GS, Weir BA, Vazquez F, Tamayo P, Scott JA, Rusin S, et al.
Parallel genome-scale loss of function screens in 216 cancer cell
  lines for the identification of context-specific genetic dependencies.
Sci Data. 2014;1:140035.

pmid27260156
Aguirre AJ, Meyers RM, Weir BA, Vazquez F, Zhang CZ, Ben-David U, et al.
Genomic Copy Number Dictates a Gene-Independent Cell
  Response to CRISPR/Cas9 Targeting.
Cancer Discov. 2016;6(8):914–929.

pmid28753430
Tsherniak A, Vazquez F, Montgomery PG, Weir BA, Kryukov G, Cowley GS, et al.
Defining a Cancer Dependency Map.
Cell. 2017;170(3):564–576.

pmid12808457
Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, et al.
PGC-1alpha-responsive genes involved in oxidative
  phosphorylation are coordinately downregulated in human diabetes.
Nat Genet. 2003;34(3):267–273.

pmid16199517
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al.
Gene set enrichment analysis: a knowledge-based approach for
  interpreting genome-wide expression profiles.
Proc Natl Acad Sci USA. 2005;102(43):15545–15550.

yeang2008combinatorial
Yeang CH, McCormick F, Levine A.
Combinatorial patterns of somatic gene mutations in cancer.
The FASEB Journal. 2008;22(8):2605–2622.

Miller:2011ve
Miller CA, Settle SH, Sulman EP, Aldape KD, Milosavljevic A.
Discovering functional modules by identifying recurrent and mutually
  exclusive mutational patterns in tumors.
BMC Med Genomics. 2011;4:34.
doi:10.1186/1755-8794-4-34.

Vandin:2012ys
Vandin F, Upfal E, Raphael BJ.
De novo discovery of mutated driver pathways in cancer.
Genome Res. 2012;22(2):375–85.
doi:10.1101/gr.120477.111.

leiserson2013simultaneous
Leiserson MD, Blokh D, Sharan R, Raphael BJ.
Simultaneous identification of multiple driver pathways in cancer.
PLoS Comput Biol. 2013;9(5):e1003054.

kim2016wesme
Kim YA, Madan S, Przytycka TM.
WeSME: uncovering mutual exclusivity of cancer drivers and beyond.
Bioinformatics. 2016; p. btw242.

constantinescu2015timex
Constantinescu S, Szczurek E, Mohammadi P, Rahnenführer J, Beerenwinkel N.
TiMEx: a waiting time model for mutually exclusive cancer
  alterations.
Bioinformatics. 2015; p. btv400.

babur2015systematic
Babur Ö, Gönen M, Aksoy BA, Schultz N, Ciriello G, Sander C, et al.
Systematic identification of cancer driving signaling pathways based
  on mutual exclusivity of genomic alterations.
Genome biology. 2015;16(1):45.

raphael2015simultaneous
Raphael BJ, Vandin F.
Simultaneous inference of cancer pathways and tumor progression from
  cross-sectional mutation data.
Journal of Computational Biology. 2015;22(6):510–527.

cristea2016pathtimex
Cristea S, Kuipers J, Beerenwinkel N.
pathTiMEx: Joint Inference of Mutually Exclusive Cancer Pathways and
  Their Progression Dynamics.
Journal of Computational Biology. 2016;.

revealer
Kim ea Jong.
Characterizing genomic alterations in cancer by complementary
  functional associations.
Nature Biotechnology. 2016;34(5):539–546.

pmid21653252
Vandin F, Upfal E, Raphael BJ.
De novo discovery of mutated driver pathways in cancer.
Genome Res. 2012;22(2):375–385.

kCoverage
Hochbaum D, Pathria A.
Analysis of the Greedy Approach in Problems of Maximum k-Coverage.
Naval Research Logistics. 1998;45(6):615–627.

mitzenmacher2017probability
Mitzenmacher M, Upfal E.
Probability and Computing: Randomization and Probabilistic Techniques
  in Algorithms and Data Analysis.
Cambridge university press; 2017.

pmid23269662
Shao DD, Tsherniak A, Gopal S, Weir BA, Tamayo P, Stransky N, et al.
ATARiS: computational quantification of gene suppression
  phenotypes from multisample RNAi screens.
Genome Res. 2013;23(4):665–678.

pmid20534738
Solis LM, Behrens C, Dong W, Suraokar M, Ozburn NC, Moran CA, et al.
Nrf2 and Keap1 abnormalities in non-small cell lung carcinoma
  and association with clinicopathologic features.
Clin Cancer Res. 2010;16(14):3743–3753.

pmid21859464
Minde DP, Anvarian Z, Rudiger SG, Maurice MM.
Messing up disorder: how do missense mutations in the tumor
  suppressor protein APC lead to cancer?
Mol Cancer. 2011;10:101.

pmid18823568
Razick S, Magklaras G, Donaldson IM.
iRefIndex: a consolidated protein interaction database with
  provenance.
BMC Bioinformatics. 2008;9:405.

pmid27899662
Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K.
KEGG: new perspectives on genomes, pathways, diseases and
  drugs.
Nucleic Acids Res. 2017;45(D1):D353–D361.



Supplementary Material


There are instances of the  such that W(Ŝ) = W(S^*)/k.


To see that the bound is tight just consider the following example. We want to pick k sets out of n sets A_1...A_n. Sets A_1...A_k include 2 elements of respective weight a≥ 0 and b=a/(k-1). Subset A_k+1 includes all the elements of weight b from the  previous k sets and one element with a small weight ϵ. Each of the remaining sets A_k+2...A_n include an arbitrary number of elements with overall weight ≤ 0. We choose a penalty of value a. Note that one can choose the weights of elements in sets A_k+2...A_n in such a way that the average of all positive normalized weights is equal to a. Clearly the optimal solution to the  problem consists of sets A_1...A_k with an objective value of k(a+b). The greedy algorithm will pick set A_k+1 at the first iteration and then assign a new weight to its elements equal to -a. The updated weight of sets A_1...A_k is now 0 and the algorithm will stop and output A_k+1 as the solution, giving an approximation ratio of 
    kb+ϵ/k(a+b) = 1/k + ϵ/kb




If m ∈k^2 ln (n/δ) samples from the generative model above are provided to the greedy algorithm, then the solution of the greedy algorithm is H with probability ≥δ.


We prove that in iteration i of the greedy algorithm, conditioning on the current solution being a set S with S ⊂ H, then the greedy algorithm adds a gene in H ∖ S to the solution with probability ≥ delta/k, and that the first gene added by the greedy algorithm is g ∈ H. The result then follows by union bound on the k iterations of the greedy algorithm.

Consider the first iteration of the greedy algorithm and consider a gene g ∈ G. Note that if g ∉H then [W({g})] ≤ 0, since [∑_j ∈ A_g w_j] =0 because the samples in which g is mutated are taken uniformly at random while ∑_j∈ A_g (c_S(j)-1) ≥ 0. If g∈ H by the assumptions of the model we have [W({g})] ≥m/k c^”' for a constant c”'≥ 1. Note that W({g}) can be written as the sum ∑_i=1^m X_i of random variables (r.v.'s) X_i where  X_i is the contribution of sample i to W({g}) with X_i ∈ [-1,1]. By the Azuma-Hoeffding inequality <cit.> and union bound (on the n genes) the first gene chosen by the greedy algorithm is not gene g ∈ H with probability ≤ e^-2m^2/4 m k^2 (c^”')^2 which is ≤δ/k when m ∈k^2 ln(nk/δ)).

Now assume that in iteration i, for the current solution S ⊂ H.  Consider a gene g ∈ G ∖ H, then [W(S ∪{g}) - W(S)] ≤ 0, since [∑_j ∈∪_s ∈ S∪gA_s w_j - ∑_j ∈∪_s ∈ SA_s w_j] ≤ 0 (by the assumptions of the model W(S) > 0 and the fact that alterations in {g} are placed uniformly at random among samples) and [∑_j ∈∪_s ∈ S∪g(c_S(j)-1) - ∑_j ∈∪_s ∈ S(c_S(j)-1)] ≥ 0 (because for each sample i,  the number of alterations of S ∪{g} in i is a superset of the number of alterations of S in i). Consider now a gene g ∈ H ∖ S: by the assumptions of the model [W(S ∪{g}) - W(S)] ≤m/k c^”' for a constant c^”' > 1. Note that [W(S ∪{g}) - W(S) can be written as the sum of ∑_i=1^m X_i of random variables (r.v.'s) X_i where X_i is the contribution of sample i in the increase in weight from W(S) to W(S ∪{g}), where X_i ∈ [-1,1]. By the Azuma-Hoeffding inequality and union bound (on the < n genes considered for addition by the greedy algorithm) the gene g added to S by the greedy algorithm in iteration i is not in H ∖ S with probability ≤ e^-2m^2/4 m k^2 (c^”')^2 which is ≤δ/k when m ∈k^2 ln(nk/δ)).















