QuaDMutEx: quadratic driver mutation explorer

Background Somatic mutations accumulate in human cells throughout life. Some may have no adverse consequences, but some of them may lead to cancer. A cancer genome is typically unstable, and thus more mutations can accumulate in the DNA of cancer cells. An ongoing problem is to figure out which mutations are drivers - play a role in oncogenesis, and which are passengers - do not play a role. One way of addressing this question is through inspection of somatic mutations in DNA of cancer samples from a cohort of patients and detection of patterns that differentiate driver from passenger mutations. Results We propose QuaDMutEx, a method that incorporates three novel elements: a new gene set penalty that includes non-linear penalization of multiple mutations in putative sets of driver genes, an ability to adjust the method to handle slow- and fast-evolving tumors, and a computationally efficient method for finding gene sets that minimize the penalty, through a combination of heuristic Monte Carlo optimization and exact binary quadratic programming. Compared to existing methods, the proposed algorithm finds sets of putative driver genes that show higher coverage and lower excess coverage in eight sets of cancer samples coming from brain, ovarian, lung, and breast tumors. Conclusions Superior ability to improve on both coverage and excess coverage on different types of cancer shows that QuaDMutEx is a tool that should be part of a state-of-the-art toolbox in the driver gene discovery pipeline. It can detect genes harboring rare driver mutations that may be missed by existing methods. QuaDMutEx is available for download from https://github.com/bokhariy/QuaDMutEx under the GNU GPLv3 license.


Background
Cancer is a complex and heterogeneous disease that starts at cellular level as a consequence of a hereditary or, most prevalently, environmentally induced mutations [1,2]. Mutations such as amino acid substitutions or copy number alterations may lead to abnormal cells that can divide indefinitely and have the ability to invade other tissues [3]. A sequence of between two and eight mutations that target genes involved in specific cell functions is needed for most human cancers to develop [4]. Such mutations, which confer growth advantage to cells and are causally implicated in oncogenesis, are referred to as driver mutations [5]. Known somatic mutations linked to cancer, often increasing availability of data, the problem of identifying driver mutations and driver genes that harbor them is far from being solved.
The main challenge is that majority of somatic mutations acquired in human cells throughout life are not causally linked to cancer. It is estimated that a typical human cell, with a genome consisting of approximately 3 × 10 9 base pairs, gains on the order of 10 −10 mutations per base pair per cell division [16,17], although the rate can vary substantially depending on factors such as local chromatin organization of the genome [18]. Human organism consists of on the order of 10 13 cells [19], many of which are in fast dividing tissues; for example around 10 11 epithelial cells are being lost and need replacement every day [20]. It is thus evident that most mutations do not lead to carcinogenesis -these are often referred to as passenger mutations. Indeed, it has been observed that in tissues that self-renew through cell division, such as skin or gastrointestinal epithelium, the number of mutations seen in cancer samples from patients 85 years old is twice the number of mutations in patients that are 25 years old. It has been estimated that half or more of all mutations observed in patients' cancer tissues originate prior to the onset of cancer [17]. In addition to these mutations, cancer cells exhibit a mutator phenotype, that is, an increased mutation rate [21], with mutation rates that can differ by an order of magnitude among subclones within the tumor [22]. This further contributes to the dominance of passenger mutations over driver mutations in observed cancer tissue samples. Altogether, while the number of driver mutations in a tumor is typically small -a recent analysis of TCGA data shows it to be between 2 and 6 in most tumors [23] -the total number of somatic mutations present in a single patient tumor tissue sample can range between 10 to above 100, depending on tissue type and patient age. Most mutations in a cancer tissue sample are thus passenger mutations that do not contribute positively to cancer growth. In fact, weakly deleterious effects of multiple passenger mutations can accumulate and can have negative impact on the tumor [24].
To discover driver mutations in the abundance of passenger mutations, many approaches take the route of calculating the background mutations rate that would be exhibited by passenger mutations, and consider those mutations that are observed more frequently as drivers. These approaches employ a statistical model of somatic mutations, typically considering them a result of a Poisson process, which allows for quantifying the statistical significance of any deviations from the background mutation rate. For example, MutSig [25] uses a constant mutation rate across all genes, and can also use methods for functional predictions of mutation significance, such as SIFT [26], CHASM [27], Polyphen-2 [28] and MutationAssessor [29]. MutSigCV [30] uses factors such as chromatin state and transcription activity to estimate gene-specific background mutation rates. PathScan [31] utlizes a Poissonian mutation model that involves gene lengths, and for a gene set given by the user calculates the probability of observing that many mutations or more under a null hypothesis that the mutations are passengers. If the probability is low across many samples, the genes are considered driver genes. MuSiC [32] extends PathScan by adding knowledge about correlation between mutation rates and factors including clinical variables such as age, molecular variables such as the Pfam family to which the genes belong, and sequence correlates such as base composition of the site and proximity among mutation sites. DrGaP tool [33] considers 11 different types of mutation types, with factors including G/C content near the mutation site and methylation status of the site, in estimating the background mutation rate. DOTS-Finder [34] integrates functional predictions and background mutation rate to identify driver genes.
Gene-centric methods for finding driver mutations from cancer sequencing data are hampered by the fact that a single driver gene is rarely mutated across many patients with a given tumor. Only few genes, such as TP53 or BRCA1, are mutated in large fraction of cases. Most driver mutations are relatively rare in tumor patients: most of individual genes are mutated in less than 5% of patients [35]. Thus, a statistically significant detection of deviation from background mutation rate requires large number of samples for rare drivers.
Observations from cancer samples show the diseaselinked mutations are not confined to a specific set of loci but, instead, they differ substantially in individual cases. Only when seen from the level of pathways, that is, genes related to a specific cellular process, a clearer picture emerges. A study of pancreatic cancer has identified a core of altered pathways common to all cases, and additional variant pathways [36] altered in some of the patients. This evidence has given rise to network-oriented driver detection methods, such as HotNet [37,38], which incorporates protein-protein networks and uses a heat diffusion process, in addition to gene mutation frequency, to detect a driver subnetwork. Some methods move beyond utilizing mutation data. For example, MEMo [39] uses gene expression to filter out genes with copy number alterations that do not show altered expression. A more refined way of incorporating gene expression data is used by Driver-Net [40], which analyzes if a mutation in a gene affects expression of genes it regulates.
In many types of tumors, only one mutation per pathway, or functionally related group of genes, is needed to drive oncogenesis [41][42][43]. Thus, the minimal set of mutated genes required for cancer to develop would consists of several sets of genes, each corresponding to a crucial pathway such as angiogenesis. Within each gene set, in each patient exactly one gene would be mutated. That is, all patients would be covered by a mutation in a gene from the set, and there would be no excess coverage, that is, no patient will have more mutations than one in the genes from the set. This pattern has been often referred to as mutual exclusivity within a gene set, and several methods, including Dendrix [44] and Multi-Dendrix [45], RME [46], CoMEt [47], TiMEx [48] and MEMo [39] detect set of driver genes by quantifying mutual exclusivity. Further methods extend these by helping deal with observation errors in the data [49], and with computational efficiency of the search for driver genes [50].
Mutual exclusivity describes the combinatorial pattern of a minimal set of genes required for oncogenesis. In actual patient data, additional mutations in driver genes may occur, especially for slow growing tumors. Also, some of the mutations may be missed due to observation errors. Thus, instead of detecting the presence or absence of mutual exclusivity in a set of genes, driver detection algorithms involve a score that penalizes for deviations from a driver pattern. That is, a penalty is incurred for zero mutations in a patient, or for more than one mutation. Then, a heuristic search procedure is utilized to find a set of genes closest to the mutual exclusivity pattern, since finding such a set has been shown to be an NP-hard problem [44].
Here, we propose a tool, QuaDMutEx, which brings three novel aspects to the mutual-exclusivity-based driver detection. First, instead of linear penalty for excess coverage used in tools like Dendrix, QuaDMutEx uses a quadratic penalty that provides a more realistic penalty for sets with excessive number of mutations. Second, the method allows for user-specified trade-off between increasing coverage and decreasing excess coverage, allowing for tailoring the method to fast-or slow-evolving tumors. Third, QuaDMutEx uses a combination of optimal search that results in globally optimal solutions to subproblems with a stochastic search through a series of subproblems, allowing for more effective search through the space of possible driver gene sets. We evaluated our method on data obtained literature and from the Cancer Genome Atlas. Our method shows higher coverage and higher mutual exclusivity than four state-of-theart tools: Dendrix, TiMEx, RME and CoMEt. Compared to DriverNet, a non-exclusivity based tool, it returns complementary sets of putative cancer driver genes of comparable quality when evaluated against the COSMIC cancer database.

Methods
The proposed algorithm for detecting driver mutations in cancer operates at the gene level. That is, on input, we are given an n by p mutation matrix G, where n is the number of cancer patients with sequenced cancer cell DNA, and p is the total number of genes explored. The matrix is binary, that is, G ij = 1 if patient i has a non-silent mutation in gene j; otherwise, G ij = 0. A row vector G i represents a row of the matrix corresponding to patient i. The solution we seek is a sparse binary vector x of length p, with x j = 1 indicating that mutations of gene j are cancer driver mutations. We will often refer to the nonzero elements of x as the mutations present in x.
In designing the algorithm for choosing the solution vector x, we assumed that any possible vector is penalized with a penalty score based on observed patterns of driver mutations in human cancers. We expect that each patient has at least one mutation in the set of genes selected in the solution; however, in some cases, the mutation may not be detected. Also, while several distinct pathways need to be mutated to result in a growing tumor, typically one mutation in each of those pathways suffices. The chances of accumulating additional mutations in the already mutated pathway are low, and decrease with each additional mutation. We capture this decreasing odds through a quadratic penalty associated with x given the observed mutations G i in patient i The term G i x captures the number of mutations from solution x present in patient i. The penalty is parameterized by a non-negative real number k to be chosen by the user. It captures the ratio of penalty for exactly two mutations in genes from set x present in patient i to penalty for no mutation from set x present in patient i. We incur no penalty if the number of mutated genes from x in a given patient is one. The effect of k on the penalty can be seen in Fig. 1. For example, for a tumor with strong mutator phenotype where more mutations are present one can set k to a low value, lowering the penalty for multiple mutations in genes from set x present in a patient.
In addition, we expect that the number of genes harboring driver mutations in a given pathway is small. Hence, we introduce a penalty on the number of genes selected in the solution, in a form of L 0 pseudo-norm, L 0 (x) = ||x|| 0 . The effect of introducing the penalty can be seen in Fig. 2.
The total penalty for a possible solution vector x is a sum of per-patient penalties and the solution-size penalty: (2) The parameter C controls the trade-off between minimization of L(G i , x) terms and of the L 0 pseudo-norm. It can alternatively be seen as the penalty incurred by increasing the number of genes in the solution x by one.
Minimization of L(G, x) can be viewed as an unconstrained binary quadratic problem (BQP) with the solution space involving binary vectors x of length p: where 1 n represents a unit vector of length n. BQPs are known to be NP-hard in general [51]. However, the optimal solution can be obtained quickly for problems of small size. Our approach in solving this problem involves a meta-heuristic based on Markov-Chain-Monte-Carlo search combined with optimal local search for small subproblems. The algorithm is presented below.
The main QuaDMutEx algorithm goes through T iterations, and in each considers a solution x containing up to ν genes. In each iteration, a new candidate solution is generated by randomly modifying the current solution vector. The new candidate solution is then modified by dropping some genes, based on exact binary quadratic optimization (Eq. 3) involving ν genes present in the candidate solution. If the optimized solution is better than the solution from previous iteration, it is accepted. If not, it is accepted with probability depending on the difference in quality of the previous and the current solution. Throughout iterations, the solution x * with the lowest value of the objective function (Eq. 2) is kept.
for t ← 1, ..., T do 5: x=RANDOMGENERATENEWSOLUTION(x t−1 , ν, ) 6: x, L=LOCALOPTIMIZESOLUTION(G, x, C, k) 7: if L < L * then 8: L * = L 9: x * = x 10: end if 11: if r < P then 14: L t = L 15: x t = x 16: else 17: 18: 20: end for 21: The random process generating a new candidate solution based on current solution always returns a solution with exactly ν genes. If the current solution already has ν genes, one of them will be randomly replaced with a gene not in the solution. The gene to be removed is chosen at random with uniform probability of 1/ν. The gene to be added is chosen by random sampling from a distribution ∼x , which is defined through a user-supplied distribution over all genes, modified to have 0 probability for the genes currently in solution x. If the current solution contains less than ν genes, the solution is expanded to include ν genes, and the ν − ||x|| 0 genes to be added are sampled without replacement according to ∼x . In our experiments, we used proportional to the logarithm of the frequency of a mutation in a given gene among patients in the dataset.
Algorithm QuaDMutEx: RandomGenerateNewSolution if ||x|| 0 = ν then 3: x=RANDOMREPLACEONE(x, ∼x ) 4: The local search for an improved new solution returns an optimized solution x and its penalty score, L. It operates by limiting the problem to the ν genes present in the new candidate solution. That is, we create a n by ν submatrix G x by choosing from G columns for which x = 1. Thus, we have an NP-hard binary QP problem with number of variables small enough that that problem can be quickly solved to the optimum using standard techniques. In our experiments, for datasets with below 1000 patients, values of ν up to 50 lead to BQP problems where global optimum could be reached in less than a second on a desktop workstation.
Algorithm QuaDMutEx: LocalOptimizeSolution In the proposed approach, the solution vector x from a single run of QuaDMutEx will capture a set of driver genes that are functionally related and thus exhibit mutual exclusivity pattern, for example genes that are all part of a pathway that needs to be mutated in oncogenesis. To uncover a comprehensive set of driver genes for a specific cancer type, spanning multiple functional subsystems vital to oncogenesis, the algorithm should be applied multiple times, each time removing the genes found in prior runs from consideration.

Evaluation on real cancer datasets
We evaluated the proposed algorithm using four somatic mutation datasets (see Table 1), one from the Cancer Genome Atlas (TCGA) database and three from literature. Two datasets were originally used by the authors of Dendrix: somatic mutations in lung cancer (LUNG), and a dataset relating to Glioblastoma Multiforme (GBM) that includes not only somatic mutations but also copy number alternations. The ovarian cancer dataset (OV) was originally used by the authors of TiMEx tool [48]. A larger dataset of mutations in samples from Breast Invasive Carcinoma (BRCA) was downloaded from TCGA. Following standard practice, in the BRCA dataset we removed known hypermutated genes that have no role in cancer [30], including olfactory receptors, mucins, and a few other genes such as titin. For each dataset, each gene in each patient was marked with one if it harbored one or more mutation, and with zero otherwise, resulting in the input matrix G for QuaDMutEx.

Quantitative evaluation of QuaDMutEx results
We ran QuaDMutEx on the four datasets: GBM, OV, LUNG, and BRCA. In the tests, we set the maximum size of the gene set to be ν = 30. We set k = 1, indicating neutral stance with respect to the trade-off between coverage and excess coverage. The value of C, the weight of the gene solution size penalty, was set to 0.5 for GBM, the dataset with the smallest number of genes measured, to 1 for the LUNG and OV datasets which have twice the number of genes compared to GBM, and to 1.5 for BRCA, the dataset with much larger number of genes. We ran QuaDMutEx for 10,000 iterations, which corresponds to running times below 10 minutes for each dataset. For GBM and BRCA, we also ran additional experiments with the default parameter values: To assess statistical significance of the results returned by QuaDMutEx, we used the method proposed in [44]. In short, we randomly permuted the contents of each column of the input patient-gene matrix, which results in randomized dataset in which, for each gene, the number of patients harboring a mutation in the gene is preserved, but any pattern of mutation within a row, that is, within each single patient, is lost. We created 1000 randomized datasets and ran QuaDMutEx on each dataset. The value of the objective function observed on the original dataset was then compared with the distribution of objective function values on the randomized datasets to obtain a p-value estimate. The results of the tests, presented in Table 2, show that for all four datasets, QuaDMu-tEx returns gene sets that are statistically significant at 0.05.
The quadratic penalty provides a single-metric measure for what is essentially a two-criterion optimization problem involving simultaneous maximization of coverage and mutual exclusivity. To capture each of these independently, we used two metrics, coverage and excess coverage:

Comparison with other mutual-exclusivity-based methods
For comparison, we used RME [46], TiMEx [48], CoMEt [47] and Dendrix [44] as they are all from the de novo discovery family of methods [11] for driver detection, and all utilize only genetic data, same as QuaDMutEx. We ran the four tools on the same four datasets: GBM, OV, LUNG, and BRCA. For TiMEx, , which does not require the user to specify the number of genes in the solution, we ran the tool with default parameters. Dendrix, RME and CoMEt require the user to provide the desired solution size. For Dendrix, we performed 29 runs for each dataset, with the solution size parameter ranging from 2 genes to 30 genes, and picked the solution  size with the best Dendrix score. Each run involved 10 7 iterations. For CoMEt the running time increases steeply with the requested solution size, thus we used sizes for which a single run finishes in less than 48 h; in result, we tested solution sizes 2, 3, 4 ,5 for GBM and OV, between 2 and 6 for LUNG, and between 2 and 10 for BRCA. For RME, we used solution sizes between 2 and 5 genes, as recommended by the authors of the tool. For BRCA dataset, RME invoked with default parameters does not return any valid solution; to circumvent this problem, we executed RME for BRCA with the minimum gene frequency parameter lowered to 0.02 from the default value of 0.1. For the other three datasets, we used the default value. We used the objective function maximized by Dendrix, which can be expressed using the notation introduced in the Methods section as Dendrix score = n− n i=1 |G i x−1|, as the metric for evaluating the tool. Essentially, the Dendrix score equals to total coverage minus coverage overlap, where total coverage is the number of patients covered by at least one gene from the given gene set, and coverage overlap is total count of all mutations in genes from the set that are in excess of one mutation per patient. High-quality solutions should have high Dendrix score.
The results of the tests, presented in Table 3, show that QuaDMutEx consistently returns higher quality solutions than all other methods. Only on the OV dataset, Dendrix discovers the same set of genes as QuaDMutEx.
Remarkably, the quality of solutions from QuaDMutEx is higher even though the score used as the metric, the Dendrix score, is not function optimized by QuaDMu-tEx, but is the objective function of Dendrix. These results show that the proposed optimization scheme that combines stochastic heuristic approach with exact solution to a series of tractable subproblems is more efficient than the heuristic approach employed in Dendrix. The putative cancer driver gene sets discovered by QuaDMutEx are mostly different than sets returned by other tools (see Fig. 3).
We also checked how QuaDMutEx performs with respect to coverage and excess coverage, and compared the results with those of Dendrix, RME, TiMEx, and CoMEt. One of the features of QuaDMutEx is the flexibility in choosing the parameter k, which controls the trade-off between high coverage but higher excess coverage solutions and low excess coverage but lower coverage solutions. Thus, we ran QuaDMutEx with a range of values of parameter k = 0.25, 0.5, 1, 1.5, 2, 2.5, 4. As previously, the value of C was set to 0.5 for GBM, to 1 for the LUNG and OV, and to 1.5 for BRCA. The number of iterations was again set to 10,000. For each parameter setting, we ran QuaDMutEx 5 times. We also gathered results from 5 runs of Dendrix for the best-performing solution size. For RME, TiMEx, and CoMEt the results do not vary from run to run, so we instead picked top five solution from a single run. Then, we quantified coverage and excess coverage. The results in Fig. 4 show that QuaDMutEx solutions are on the Pareto-optimality frontier of all (RME, TiMEx, CoMEt, Dendrix and QuaD-MutEx) solutions. For each Dendrix, TiMEx and CoMEt solution, there is a QuaDMutEx solution that is better: has higher coverage and lower excess coverage. These results further confirm results from Table 3 showing that the proposed tool improves upon the state-of-the-art. Data for OV are not shown graphically, as there is very little  Table 3.

Effects of parameters on QuaDMutEx
The proposed methods allows for adjusting the penalty for expanding the solution size, through a parameter C that corresponds to the additional penalty for increasing the number of genes in the solution by one. It also allows for tweaking the trade-off between coverage and mutual exclusivity, through a parameter k that captures the ratio of penalty for one excess mutation in a patient to penalty for the patient not being covered by any mutation. We have analyzed how these two parameters affect the solution by running QuaDMutEx for 10,000 iterations for parameters C = 0.25, 0.5, 1, 1.5, 2, 2.5, 4 and k = 0.25, 0.5, 1, 1.5, 2, 2.5, 4. Figure 5 shows that the parameter C achieves its design goal, that is, solutions with higher C include fewer genes. The figures also show that as the penalty for the size of the solution set is lowered, by specifying lower value of C, the coverage of patients by genes in the solution tends to increase for the three small datasets, where high values of C reduce the solution size to only a few genes and thus necessarily lower coverage. This effect is not present in the large dataset, BRCA, where C does not impact coverage. Changing C does not show any impact on excess coverage.
Changes in parameter k result in changes in coverage and excess coverage, but has no substantial impact on the number of genes in the solution. The results show that, as intended, lower values of k lead to higher coverage, at the cost of higher excess coverage, than high values of k. Thus, for slow growing tumors, tumors with elevated mutator phenotypes, or tumors in old patients, where many mutations may occur by chance and higher excess coverage is expected, low values of k is preferred over high k values.

Qualitative assessment of QuaDMutEx results
To validate the ability of QuaDMutEx to take only mutation data and discover rare putative cancer driver genes, which are the most hard to find using traditional methods that rely on mutation frequency in patient population, in each of the four datasets we focused on the genes in the solution with the fewest number of mutations. See Table 4 for a complete list of all genes in the solution, and for the number of mutations for each gene in each dataset. In addition to literature review, we also used DriverDBv2 [7], a database of previously discovered cancer driver genes, to further validate the quality of QuaDMutEx solutions.
In the brain tumor dataset, eight identified genes are each mutated in only 1 out of 84 patients. Out of these, ITGB3 has known role in multiple cancers [52,53], TRIM2 has tumor suppressing function in ovarian cancer [54] and plays a role in brain, the source of the analyzed tissue [55], WEE1 is already a target for cancer therapy [56], and CHD5 is a known tumor suppressor [57]. Changes in expression of MARK4 have been observed in glioblastomas [58]. While no cancer role has been so far identified for carboxylesterase 3 (CES3), it is known to be expressed in the source tissue of our samples, the brain [59]. SHH gene has been linked to glioma growth [60], as well as to other cancers [61]. Finally, IQGAP1 is believed to play a role in cell proliferation and cancer transformation [62]. In the ovarian cancer dataset, KRAS, a known protooncogene, was found mutated in two patient. Eosinophil cationic protein (RNase 3) was found in only one patient. The protein, while not present in DriverDBv2 and not directly related to oncogenesis, has cytotoxic activity and was recently shown to inversely affect viability of cancer cell lines [63] and thus its mutations may affect human tumor growth.
In the QuaDMutEx solution for the lung datasets, six putative cancer driver genes are each mutated in only two of the 356 patients, and additional four are mutated in single patients. Among these, role of ABL1 in cancer is well established. PAK6 has been shown to be involved in prostate cancer [64], and presence of MAST1 mutations has been detected in lung samples [65]. The expression of CYSLTR2 gene is a prognostic marker in colon cancer [66]. RPS6KA2 gene is a putative tumor suppressor gene in ovarian cancer [67], and FES is a known protooncogene [68]. BAX is an oncoprotein with known role in cancers [69], including lung cancer [70]. Mutations in the PIK3C2B gene were previously observed in lung and other tumors [71,72]. There is emerging evidence of a role of RANBP9 gene in lung cancer [73]. The 67-kDA laminin receptor gene RPSA, while not present in DriverDBv2, is known to play a role in tumor growth [74,75]. Among the putative driver genes discovered by QuaD-MutEx in the BRCA samples, nine were mutated in four or fewer of the 771 patients. Two among the genes that were mutated in more than four patients were not present in the DriverDBv2 database: NBPF1 and KIAA1310. However, NBPF1 has recently been identified as tumor suppressor gene [76]. KIAA1310 (KANSL3) is a member of KANSL family which plays a role in cell cycle and reduction of its function is associated with cancer [77]. Of the rarely mutated genes, only TSNARE1 gene is likely to be a false positive. CASP8AP2 gene has been previously linked to cancer [78,79]. No direct role in oncogenesis for ADCY1 gene has been reported, however it has been found downregulated in osteosarcomas [80]. PITX2 is a recurrence marker in breast cancer [81]. PSG11 gene has been shown to be correlated with survival in ovarian cancer [82]. Ankyrin repeat proteins, though not ANKRD34B specifically, have been previously reported as promoting cancer development [83]. KRT14 gene dysregulation was recently linked with breast cancer metastases [84]. MSI1 is putative therapeutic target in colon cancer [85]. TWISTNB is a component of the RNA polymerase I complex, and while TWISTNB gene has not been previously linked to cancer, mutations in polymerase subunits, cofactors, and mediators are known factors in malignancy [86]. Together, these results confirm that QuaDMutEx is effective in identifying cancer driver mutations even if they are rare in the analyzed patient group.

Comparison with gene expression-based driver discovery
In addition to methods that use only genomic mutation data, we also compared QuaDMutEx to Driver-Net, a method that uses a biological network and gene expression data in addition to mutation data. We used four genomic-transcriptomic datasets that are provided with the DriverNet tool: triple negative breast cancer (eTNB), glioblastoma multiforme (eGBM), high-grade serous ovarian cancer (eHGS), and METABRIC breast cancer (eMTB) datasets. The summaries of the datasets are provided in Table 5.
DriverNet was executed using default parameters on the full information contained in the dataset, that is, the genomic, transcriptomic, and biological network information. The solution gene sets include all genes found by DriverNet to be statistically significant at the 0.05 p-value threshold. QuaDMutEx was executed using only the genomic data describing presence or absence of a mutation in a given gene in a given patient. We used the default value of k = 1, and set the value of C to 1.5, with the exception of the smallest dataset, eTNB, for which we used C = 1. We compared the putative cancer driver gene sets discovered by the two tools using coverage, excess coverage, and the Dendrix score, as described above.
For the eGBM dataset, QuaDMutEx shows much higher coverage and much lower excess coverage (see Table 6). For the other three datasets, QuaDMutEx shows much lower excess coverage than DriverNet, at the cost of a moderate decrease in coverage. These results reflect the fact that DriverNet is not designed to take mutual exclusivity of genes into consideration. On the other hand, DriverNet return many more genes than QuaDMutEx. A single run of QuaDMutEx is designed to return a single set of genes with low excess coverage, and does not include all putative driver genes -these can be detected with another run of QuaDMutEx.
To provide a comparison that does not involve mutual exclusivity, we used the COSMIC database of mutations in cancer, and we introduced iterated QuaDMutEx,  which increases the number of genes found by QuaDMu-tEx to the numbers similar to DriverNet. We performed four executions of QuaDMutEx, after each run removing the genes discovered so far from the dataset, so that they do not prevent discovery of additional genes that are not mutually exclusive with previously discovered ones. We then pooled the four high-exclusivity gene sets into a single high-coverage set. Since mutual exclusivity can be expected only for a set of functionally-related genes, for example genes from a single cancer-related pathway, a single call to QuaDMutEx corresponds to a single-pathway query, and calling QuaDMutEx iteratively corresponds to a multi-pathway query, facilitating comparison with DriverNet which does not have a single-pathway focus.
To measure the quality of solutions returned by Driver-Net and iterated QuaDMutEx in a way independent of any mutual exclusivity of gene mutations, we compared the numbers of COSMIC occurrences of mutations in genes returned by DriverNet with occurrence numbers for QuaDMutEx gene sets. Specifically, for each gene in a discovered gene set, we queried COSMIC for the number of observed mutations in that gene. We then plotted a complementary cumulative distribution function (CCDF) over the numbers over the whole gene set. For example, for the eHGS dataset, for both QuaDMutEx and DriverNet, the CCDF value at 1000 is approximately 0.14, indicating that for both methods, 14% of the genes in the solution set have more than 1000 mutation each in COSMIC, while for 86% of genes in the solution set a COSMIC query for the gene results in at most 1000 mutations. The results in Fig. 6 indicate that iterated QuaDMutEx and DriverNet perform similarly on eTNB and eGBM datasets, and on eHGS and eMTB both perform similarly for majority of the mutation counts range, with DriverNet having an edge at the numbers below that threshold.
Genes returned by QuaDMutEx are to large extent different than those returned by DriverNet (see Fig. 7 showing that the expression-based approach used in DriverNet and the mutation-only approach used in QuaDMutEx are complementary. We validated the genes discovered by QuaDMutEx (Table 7) in DriverDB2, a database of genes previously discovered as cancer drivers. For eTNB and eGBM datasets, all the genes discovered by QuaDMutEx are present in DriverDB2 database. In eHGS dataset, only ANKRD36B was not found in DriverDB2. However, ANKRD36B gene was identified in rare germline copy number variations in renal clear cell carcinoma [87], and also correlates with cellular sensitivity to chemotherapeutic agents [88]. In eMTB dataset, TRA@ gene is not present in DriverDB2, but it has been previously found to be linked to breast cancer [89]. TRA@ os also one of the genes that were discovered both by DriverNet and by QuaDMutEx. TBC1D3P2 is recurrently mutated in meningioma cell lines [90] and is a pseudogene for TBC1D3, a known oncogene [91].
There is no information available about AC116655.7-12 and AC116165.7-3, and at this point we classify both as false positives.

Conclusions
Superior ability to improve on both coverage and excess coverage of the detected driver gen sets on datasets from different types of cancer shows that QuaDMutEx is a tool  For each gene, in parentheses, we provide the number of patients in the dataset that harbored a mutation in that gene. Genes in bold are present in the DriverDBv2 [7] database of previously discovered cancer drivers that should be part of a state-of-the-art toolbox in the driver gene discovery pipeline. It can help detect lowfrequency driver genes that can be missed by existing methods.