 METHOD
 Open Access
 Published:
A novel independence test for somatic alterations in cancer shows that biology drives mutual exclusivity but chance explains most cooccurrence
Genome Biology volume 17, Article number: 261 (2016)
Abstract
In cancer, mutually exclusive or cooccurring somatic alterations across genes can suggest functional interactions. Existing tests for such patterns make the unrealistic assumption of identical gene alteration probabilities across tumors. We present Discrete Independence Statistic Controlling for Observations with Varying Event Rates (DISCOVER), a novel test that is more sensitive than other methods and controls its false positive rate. A pancancer analysis using DISCOVER finds no evidence for widespread cooccurrence, and most cooccurrences previously detected do not exceed expectation by chance. Many mutual exclusivities are identified involving wellknown genes related to cell cycle and growth factor signaling, as well as lesser known regulators of Hedgehog signaling.
Background
Tumor development emerges from a gradual accumulation of somatic alterations that together enable malignant growth. As has been revealed by recent genomic profiling efforts, an immense diversity exists in the alterations that tumors acquire [1, 2]. Whether by e.g., copy number aberration, point mutation, or DNA methylation, alterations of many genes may potentially trigger transformation. Often though, the fate of a cell acquiring a certain alteration depends on other alterations already present [3]. Therefore, with an everexpanding catalog of cancer genes, a need arises to establish how alterations in those genes interact to transform healthy cells to cancer cells. This task can be approached by statistical analyses aiming to uncover more complex, combinatorial patterns in somatic alterations.
Two such patterns are cooccurrence and mutual exclusivity. In the former, alterations of certain combinations of genes tend to coexist in the same tumor, whereas in the latter, mostly only one out of a group of genes is altered in a single tumor. Mutual exclusivity is frequently observed in cancer genomics data [4, 5]. Individual alterations targeting similar biological processes are believed to be mutually redundant, with one alteration being sufficient to deregulate the affected process. Identifying mutual exclusivity can therefore help in finding unknown functional interactions. With this in mind, several statistical methods have been proposed to identify significant patterns of mutual exclusivity [6–12].
Just as mutual exclusivity is interpreted as a sign of redundancy, cooccurrence is often held to entail synergy. Alteration of only one of the two genes would be relatively harmless, whereas cells with alterations in both progress to malignancy. If such synergy exists, cancer genomes should be enriched for these coalterations; i.e., tumors harboring alterations in both genes should be more frequent than expected by chance. Several studies have reported an abundance of cooccurring somatic alterations in various types of cancer [13–19]. For somatic copy number changes, however, it has also been suggested that cooccurring alterations emerge from tumors’ overall levels of genomic disruption [20]. Indeed, tumors display a wide diversity in genomic instability, both across and within cancer types. In tumors harboring many alterations, one should not be surprised to see simultaneous alterations in any pair of genes. In contrast, two genes altered in a tumor carrying a small number of alterations might instead have resulted from a purifying selective process. Suggesting synergy as an explanation for observed cooccurrence is only reasonable if a simpler explanation like tumorspecific alteration rates can be rejected.
In this paper, we address the statistical implications of heterogeneous alteration rates across tumors for cooccurrence and mutual exclusivity detection. With extensive analyses of simulated data, we show how commonly used statistical tests are not equipped to deal with the mismatch between what is assumed by the test and what is encountered in the data. In the presence of heterogeneous alteration rates, countless spurious cooccurrences are picked up in data that are controlled not to contain any. At the same time, many instances of true mutual exclusivity are missed. Based on these observations, we introduce DISCOVER, a novel statistical independence test that incorporates the overall alteration rates of tumors to successfully solve the issues encountered with existing tests. We compared the performance of DISCOVER to that of several other published mutual exclusivity tests: MEMo [6], muex [8], mutex [9], CoMEt [10], MEGSA [11], and TiMEx [12]. Across the whole range of significance levels, DISCOVER is more sensitive while controlling the false positive rate at the specified level.
We also applied DISCOVER to a selection of more than 3000 tumors across 12 different cancer types. Only one cooccurrence was detected that is not explained by overall rates of alteration alone. On the other hand, many more cases of mutual exclusivity were detected than would have been possible with traditional tests. The genes targeted by these alterations cover many of the core cancer pathways known to display such exclusivity. However, we also identified exclusivity among less canonical actors in the cell cycle, and among regulators of Hedgehog signaling.
Results
Common tests for cooccurrence or mutual exclusivity assume homogeneous alteration rates
A commonly used test for both cooccurrence and mutual exclusivity is Fisher’s exact test applied to a 2×2 contingency table [16–18]. The test is used to support cooccurrence when the number of tumors with alterations in both genes is significantly higher than expected by chance. Likewise, it suggests mutual exclusivity when the number of tumors with alterations in both genes is significantly lower. The validity of this test depends on the assumption that genes’ alterations across tumors are independent and identically distributed (i.i.d.). Identical distribution implies that the probability of an alteration in a gene is the same for any given tumor. With cancer’s heterogeneity in mind, this assumption may prove problematic. Surely, a gene is more likely found altered in tumors with many somatic alterations overall, than in tumors with only few such changes.
Other tests used for cooccurrence or mutual exclusivity depend on the same i.i.d. assumption as described for Fisher’s exact test. This is the case for permutation tests that estimate the expected number of tumors altered in both genes by randomly reassigning gene alterations across tumors [7, 13]. It is also true for a simple binomial test that we will use to illustrate the consequences of violating the i.i.d. assumption. This test is depicted in Fig. 1 c. The alteration probability p _{ i } of a gene is estimated to be the proportion of tumors altered in that gene. For example, gene 3 in Fig. 1 a is altered in 2 of the 5 tumors, resulting in p _{3}=0.4 (Fig. 1 c). If alterations targeting two genes are independent, the probability of a tumor altered in both genes equals the product p _{1}·p _{2} of those genes’ alteration probabilities. Hence, out of m tumors, m·p _{1} p _{2} tumors are expected to harbor alterations in both genes. In the example in Fig. 1 a, the probability of alterations in both genes 3 and 5 would be p _{3}·p _{5}=0.4·0.4=0.16. Therefore, if alterations of genes 3 and 5 were independent, we would expect 5·0.16=0.8 tumors with alterations in both. Observing more such tumors suggests cooccurrence, whereas observing fewer suggests mutual exclusivity (Fig. 1 b).
Assuming homogeneous alteration rates leads to invalid significance estimates
To illustrate the effect of the i.i.d. assumption on the detection of mutual exclusivities and cooccurrences, we performed analyses on simulated data. Genomic alterations were generated such that the alteration frequencies both per gene and per tumor resemble those observed in real tumors, but without any designed relation between the genes’ alterations; i.e., genes were simulated to be independent. As these simulated data do not contain cooccurrences or mutual exclusivities, all identified departures from independence are by definition spurious. We can therefore use these data to check the validity of the binomial test. When testing many pairs of independently altered genes, a valid statistical test should produce P values that approximately follow a uniform distribution. In contrast, when we test for cooccurrence in these data, the Pvalue distribution shows a large skew towards extremely low values (Fig. 2 a). Even highly conservative significance levels will mark the majority of gene pairs as significant hits. Given that no true cooccurrences exist in the simulated data, all these hits are false positives. If we test for mutual exclusivities instead, we observe a skew towards the high end of the Pvalue spectrum (Fig. 2 c).
We next evaluated the sensitivity of the binomial test. For this, we tested simulated cooccurrences and mutual exclusivities, which we added to the data. A sensitive test should produce only low P values for these positive cases, and so the resulting Pvalue distribution should be heavily skewed towards zero. If we test for cooccurrences, this is indeed the case (Fig. 2 b). Testing for mutual exclusivity, however, reveals a distribution that, although skewed towards lower P values, is much more stretched out across the [0, 1] interval (Fig. 2 d). Even highly liberal significance levels will only recover a small part of the positive cases.
We conclude that the binomial test is anticonservative as a cooccurrence test. In contrast, as a mutual exclusivity test, it is conservative. While we used the binomial test for this illustration, we found the same to be true for Fisher’s exact test (Additional file 1: Figure S1). To confirm our hypothesis that the i.i.d. assumption is causal to this incorrect behavior, we generated additional simulated data, making sure that the overall alteration rate was similar across the tumors. Using the binomial test to detect cooccurrence and mutual exclusivity of independent genes results in Pvalue distributions that are much closer to uniform (Additional file 1: Figure S2). This confirms that statistical tests that rely on the i.i.d. assumption are not suited for cooccurrence analysis, and have reduced sensitivity for mutual exclusivity analysis.
A novel statistical test for cooccurrence and mutual exclusivity
Our new method, which we call Discrete Independence Statistic Controlling for Observations with Varying Event Rates (DISCOVER), is a statistical independence test that does not assume identically distributed events. The main ingredients of the method are depicted in Fig. 1 d. Unlike the method in the simpler binomial test, we allow different tumors to have different alteration probabilities for the same gene—the alteration probabilities for genes 3 and 5 in Fig. 1 d now vary per tumor, in contrast to Fig. 1 c. For tumors with many altered genes, this probability is higher than for tumors with only few alterations. To estimate these alteration probabilities, we solve a constrained optimization problem that ensures that the probabilities are consistent with both the observed number of alterations per gene and the observed number of alterations per tumor. The probability of concurrent alterations in two independent genes is then obtained for each tumor individually, by multiplying the tumorspecific gene alteration probabilities, as indicated in the right panel of Fig. 1 d. With these probabilities, an analytical test based on the Poissonbinomial distribution can be performed to decide whether the number of tumors altered in both genes deviates from the expectation.
We repeated the simulation study performed for the binomial test, this time applying the DISCOVER test. First, our data only contained independently generated alterations. Testing for cooccurrence (Fig. 2 e) and mutual exclusivity (Fig. 2 g) resulted in Pvalue distributions much closer to uniform, as one would expect. The fact that these distributions are not truly uniform is a property shared by all discrete test statistics [21]; it makes discrete tests slightly more conservative. Most importantly, the anticonservative bias towards cooccurrence of the binomial test is not present in the DISCOVER test. By testing simulated cooccurrences, we established that the removal of the anticonservative bias does not compromise the sensitivity for true cooccurrences (Fig. 2 f). Moreover, the sensitivity for mutual exclusivities is improved when compared with the binomial test (Fig. 2 h).
Extension to a groupbased mutual exclusivity test
Mutual exclusivity is not restricted to pairs of genes. Larger groups of genes may also display alteration patterns in which most tumors only have an alteration in one of the genes. We considered three statistics to assess the mutual exclusivity of groups of genes: coverage, exclusivity, and impurity (Fig. 3 a). For all three of these statistics, its expectation for groups of independent genes can be described by a Poissonbinomial distribution (see Methods), and thus a statistical test can be formulated for determining significance. Based on simulated data, we established that the impuritybased group test has the best balance between sensitivity and specificity (Additional file 1: Figure S3).
Comparison to other mutual exclusivity tests
We compared the performance of the groupbased DISCOVER test to that of several other published mutual exclusivity tests: MEMo [6], muex [8], mutex [9], CoMEt [10], MEGSA [11], and TiMEx [12]. In this comparison, we focused on the statistical tests for mutual exclusivity provided by these methods (see Methods). Although the tests differ in the statistical model upon which they are based, all but MEMo assume identical alteration probabilities across tumors. Like Fisher’s exact test and the binomial test, they are thus examples of tests based on the i.i.d. assumption. MEMo does take into account tumorspecific alteration rates by preserving these rates in a permutation scheme. Unlike DISCOVER, it estimates the alteration rate with respect to a small set of recurrently altered genes as opposed to all genes.
The comparison was performed on simulated data. Groups of genes with mutually exclusive alterations of various degrees of impurity served as positive examples (see Methods). For each such group, we also selected groups of independent genes of the same size and matched to have similar alteration frequencies, to serve as negative examples. In total, 10 data sets of 100 positive and 100 negative groups were generated, and evaluation metrics were averaged across these 10 sets. We evaluated the tests for both specificity and sensitivity.
To evaluate specificity, we considered the extent to which a chosen significance level α predicts the false positive rate obtained when groups with a nominal P value less than α are classified as mutually exclusive. By definition of the P value, rejecting the null hypothesis at a significance level α should guarantee that the false positive rate (or type I error rate in statistical terminology) is at most α. Graphically, if the false positive rate is plotted as a function of the significance level (Fig. 3 b), the resulting curve would ideally follow the diagonal, or it should drop below the diagonal for more conservative tests. With the exception of muex, all methods control their false positive rate below the nominal significance level, but they do so in notably different ways. CoMEt, mutex, and TiMEx only yield false positives at extremely high significance levels. Doing so, they are more conservative than required. In contrast, DISCOVER’s curve follows the diagonal more closely. This is another confirmation that tests based on the i.i.d. assumption—like before with the binomial and Fisher’s exact tests—are more conservative than those that model the varying alteration rates. Indeed, MEMo is also less conservative than CoMEt, mutex, and TiMEx. It is more conservative than DISCOVER though, which may be explained by the different strategies for estimating the tumorspecific alteration rates: based on all genes for DISCOVER, or based on frequently altered genes only for MEMo.
To evaluate sensitivity, we compared the increase of the true positive rate as a function of the significance level (Fig. 3 c). A sensitive test will already attain high true positive rates at low significance levels. Across the whole range of significance levels, DISCOVER was found to be more sensitive than any of the other tests. It identified more mutually exclusive groups at lower significance levels. Only muex initially shows a higher sensitivity, but it does so at the price of many false positives (Fig. 3 b)—we suspect this is partly due to numerical imprecision. At higher significance levels, muex’s sensitivity drops below that of DISCOVER. MEMo only attains a high sensitivity at higher significance levels: it is affected by the limited resolution of its permutation test. We used 10,000 permutations, which makes the lowest possible P value 1×10^{−4}. Again contrasting tests based on their underlying assumption, we conclude that the conservatism caused by the i.i.d. assumption is reflected in a lower sensitivity. The majority of mutually exclusive groups are only identified at relatively high significance levels. If correction for multiple testing is applied, this may render many of them insignificant.
Cooccurrence and mutual exclusivity in pancancer somatic alterations
We analyzed a set of 3386 tumors covering the 12 cancer types studied in the TCGA pancancer initiative [22]. An alteration matrix was constructed from recurrent copy number changes and highconfidence mutational drivers. Copy number changes were analyzed for 118 genes, of which 40 were gains and 78 were losses. In addition, mutation data were added for 286 genes previously classified as highconfidence driver genes [23]. In total 404 genomic alterations were analyzed covering 374 unique genes, as 30 genes are frequently targeted by both copy number changes and mutations.
We tested for pairwise cooccurrence and mutual exclusivity between pairs of genes not located on the same chromosome. These tests were stratified for cancer type to avoid confounding due to cancer typespecific alteration frequencies. Complementing the pairwise tests, we also employed the DISCOVER group test to detect patterns of mutual exclusivity in larger groups of genes. The groups we tested were selected using two different approaches. In the first approach, we extracted gene sets from the canonical pathway collection of MSigDB [24]. We tested 23 such gene sets based on pathway membership. In the second approach, we aimed to detect de novo gene sets purely based on the data. For this, we applied a clustering algorithm to the pairwise mutual exclusivity results to identify groups of genes showing a high degree of interaction.
No evidence for widespread cooccurrence
A remarkable outcome of our analysis is that we found no evidence for widespread cooccurrence of somatic alterations. At a maximum false discovery rate (FDR) of 1%, no significant cooccurrences were identified. Relaxing the FDR threshold to 3%, we could recover one cooccurrence, between mutation of TP53 and amplification of MYC. It was recently suggested that MYCamplified tumors show higher levels of MYC expression in tumors with a TP53 mutation than in tumors without [25]. No further, reasonable relaxation of the significance threshold led to additional hits. Certainly, more gene pairs exist that harbor alterations in overlapping sets of tumors. Yet, the sizes of those overlaps do not exceed what is expected by chance if differences in tumorspecific alteration rates are taken into account. This is in sharp contrast with the significance estimates obtained with the binomial test, which identifies 21,627 significant cooccurrences, almost onethird of all pairs tested.
With the aim of establishing that the DISCOVER test is not overly conservative, we tested for cooccurrence between copy number changes of genes on the same chromosomes. Due to the inherent correlation in copy number of genes situated close to each other, such gene pairs can be considered positive controls. Indeed, all but one of the 112 pairs of tested genes located in the same recurrently altered segment are identified as cooccurring by the DISCOVER test. In addition, 18 pairs of genes situated on the same chromosome arm are detected as cooccurring, as are DDAH1 on 1p22 and MCL1 on 1q21. More generally, pairs within the same segment are assigned lower P values on average than are pairs within the same chromosome arm (P=7×10^{−39}, Additional file 1: Figure S4). The same is true, to lesser extents, for pairs within the same chromosome arm compared to pairs within the same chromosome (P=6×10^{−8}) and for pairs within the same chromosome compared to pairs across chromosomes (P=0.0004).
Mutually exclusive alterations target core cancer pathways
Pairwise mutual exclusivities were found among 181 pairs of genes, at a maximum FDR of 1% (Additional file 2: Table S1). We once more confirmed that detecting mutual exclusivities using the binomial test results in far fewer significant mutual exclusivities—only three pairs were identified. Among the 181 gene pairs, there were 107 unique genes. Many of these are significantly mutually exclusive with only one or a few other genes. For some, reduced statistical power due to low alteration frequency may be the reason for not detecting more associations. However, alteration frequency is not the dominant factor in how often mutual exclusivity is detected (Fig. 4 a). For example, mutations of KRAS are far less frequent than TP53 or PIK3CA mutations. Yet, KRAS was found mutually exclusive with more genes than were the latter two genes.
Since mutual exclusivity is believed often to occur between functionally related genes, we determined the overlap of the identified gene pairs with the STRING functional interaction network [26]. Thirtyone of the identified gene pairs have a highconfidence functional interaction in STRING (Fig. 4 b). This overlap is significantly higher than the 5 overlapping pairs expected by chance (P<1×10^{−4}), as determined using a permutation test. Moreover, 121 of the mutually exclusive gene pairs share a common interactor in the STRING network. By chance, this is only expected to be the case for 80 gene pairs (P=0.003). This suggests that the mutual exclusivities identified are indeed for a large part driven by biological factors. Another confirmation of this is found in the results of the MSigDB gene set tests (Additional file 1: Figure S5). Twelve gene sets representing several cancerrelated pathways show significant mutual exclusivity. The mutual exclusivities that overlap with STRING interactions revolve around three commonly deregulated processes in cancer: growth factor signaling, cell cycle control, and p53 signaling.
Growth factor signaling
Genes coding for proteins involved in growth factor signaling are frequently altered in cancer. These alterations display a high degree of mutual exclusivity. Mutations targeting the receptor EGFR are mutually exclusive with mutations in its downstream mediator KRAS. In turn, KRAS mutations are mutually exclusive with mutations in its family member NRAS, its negative regulator NF1, and its downstream effector BRAF. All of these alterations are able to deregulate RAS signaling, and one is sufficient. Mutual exclusivity of mutations in KRAS and mutations in both PIK3R1 and PIK3CG may be driven by the known crosstalk between RAS signaling and phosphoinositide 3kinase (PI3K) signaling [27].
The PI3K signaling cascade itself is also characterized by many mutually exclusive alterations. Mutations in the PIK3CA and PIK3R1 genes—both coding for components of the PI3K complex—are mutually exclusive. Alterations in the PTEN gene—a negative regulator of the downstream activation of AKT by PI3K—are mutually exclusive with mutations in PIK3CA, but also with alterations in the upstream activator of the cascade ERBB2. PI3K signaling is also the central biological process in several of the gene sets found mutually exclusive with the groupbased test (Fig. 5 a, Additional file 1: Figure S5). Central genes in PI3K signaling such as SOS1, AKT1, and AKT3 were not found as mutually exclusive with other pathway members in the pairwise analysis, yet the groupwise test correctly detects it.
Cell cycle control
Many tumors harbor alterations that disable the cell cycle control present in healthy cells. This control arises from a tightly regulated interplay between cell cycleactivating cyclins and CDKs, and CDK inhibitors, linked together by the master cell cycle regulator RB1. Alterations in these genes are also mutually exclusive. For example, copy number gains in Cyclins D1 and E1 are mutually exclusive, as are CDKN2A copy number loss and both mutation and copy number loss of RB1. The transcriptional activation of CCND1 by MYC is also reflected in the mutual exclusivity between copy number gains in the two genes. Also as a group, cyclins, CDKs, and CDK inhibitors show a clear pattern of mutual exclusivity (Fig. 5 b, Additional file 1: Figure S5). CDK4 and CDKN1B, central players in the regulation of the cell cycle, did not appear in the pairwise results, but are highly exclusive with the other genes involved.
p53 signaling
p53 plays a pivotal role in deciding on cell fate after cellular stresses common in cancer development. For this reason, p53 mutations are the most common alterations in cancer. However, not all tumors disable p53 function genetically. Alterations in regulators of p53 provide an alternative way to deregulate p53 function in p53wildtype tumors, but are likely redundant in tumors that already have a dysfunctional p53 protein. Indeed, we found alterations in several regulators of p53 to be mutually exclusive with TP53 mutation. For example, mutations in its positive regulator ATM, but also mutations in its negative regulator HUWE1 are mutually exclusive with TP53 mutations. MDM2 and MDM4, highly similar negative regulators of p53, have a mutually exclusive pattern of copy number gains. Mutations in CASP8, a downstream mediator of p53induced apoptosis, tend also not to overlap with TP53 mutations.
De novo gene set detection
As a final step in our analysis, we detected de novo gene sets purely based on observed patterns of mutual exclusivity, without input based on recorded biological knowledge. To this end, we applied correlation clustering to a network derived from pairwise mutual exclusivities (see Methods). This identified 120 candidate mutually exclusive gene sets. Testing these gene sets with DISCOVER, 43 were found to be mutually exclusive at a maximum FDR of 1%. The full results are presented in the online Jupyter notebooks (see Availability of data and materials). Below, we discuss two interesting examples.
One of the most significant gene sets includes RB1 and CDKN2A, two pivotal players in cell cycle control (Fig. 5 c). PARK2 [28], WWOX [29], FHIT [30], PTPRD [31, 32], and MAPK12 [33] have also all been linked to a regulating role in various phases of the cell cycle. They have been found to do so by regulating cyclins, CDKs, or CDK inhibitors. This functional similarity may explain these genes’ mutual exclusivity with RB1 and CDKN2A. As of yet, LRP1B and CSMD1 have not been linked to cell cycle control. Their mutual exclusivity with respect to several regulators of the cell cycle may instigate further study in this direction.
Another group of genes with a high degree of mutual exclusivity (P=7×10^{−8}) consists of genes that have been implicated in the regulation of Hedgehog signaling (Fig. 5 d). With the exception of ARHGAP35, all genes in this group have experimentally been linked to a regulatory role in Hedgehog signaling. GNAS [34, 35], TBX3 [36], and WT1 [37] were found to directly regulate the pathway. ARID1A, coding for a component of the SWI/SNF complex, is likely to play a similar role, since loss of another component of this complex, Snf5, was found to lead to activation of the Hedgehog pathway [38]. Besides these two examples, several other gene sets were identified that combine known interaction partners with interesting leads for undiscovered interactions.
Discussion
The recent growth in the number of large genomics data sets gives rise to a parallel increase in statistical power to detect ever more complex associations. However, as another consequence of larger sample sizes, poorly matched assumptions will have an increasing impact on the results. A central assumption behind commonly used statistical tests for cooccurrence and mutual exclusivity is that a gene’s alteration probability is identical across all tumors. Using simulated data, we have shown that this assumption is not only unjustified, but that it leads to a full reversal of the associations. The binomial test we used for illustration is but a representative of a larger class of independence tests based on the same assumption. This class includes analytical approaches such as Fisher’s exact test, CoMEt [10], and MEGSA [11], but also permutation tests where gene alterations are uniformly shuffled across the tumors.
We have presented a novel independence test based on assumptions that better match the reality of cancer genomics data. With this new test, we analyzed tumors across 12 different cancer types for the presence of cooccurrence and mutual exclusivity. Only one case of cooccurrence was found, whereas numerous cases of mutual exclusivity were detected. Performing the same analysis with the binomial test led to the detection of many cooccurrences and almost no mutual exclusivity. Many of the mutual exclusivities missed by the binomial test can be related to central processes in cancer biology. We found strong mutual exclusivity between genes involved in growth factor signaling and cell cycle control. Also, lesser known players in the regulation of cell cycle and Hedgehog signaling were identified. Based on the results of our simulation study, we are confident that most of the cooccurrences detected by the binomial test are spurious.
The absence of widespread cooccurrence contradicts what was found in previous genomewide studies. Besides, it seems counter to our expectation of positive selection for synergy that led us to look for cooccurrence in the first place. It is true that synergy resulting from the alteration of multiple genes has been observed. Comutation of genes has been reported to act on a tumor’s response to chemotherapy, or more generally on patient survival [39, 40]. None of these phenotypes, however, has been the subject of the selection from which the original tumor emerged. Only after selective pressure for that particular phenotype has taken place—for example, by treating patients—would enrichment for such cooccurrences be detected. There is no doubt that cancerdriving alterations often act in concert. Yet if statistical results are to serve as support for, or even meant to identify synergy, other possible explanations for the observed cooccurrence should be accounted for. In our pancancer analysis, overall alteration rates explained most if not all cooccurrence.
The need to take into account higher level structural features of samples is not unique for cooccurrence and mutual exclusivity analysis. In testing the relationship between highdimensional gene expression data and phenotypes of interest, latent sources of heterogeneity can have a profound effect on the results. Approaches like surrogate variable analysis [41] have been developed to adjust analyses appropriately. Similarly, genomewide association studies face the issue of latent population substructure. Again, if ignored, such substructure can drastically alter the findings. Linear mixed models have gained popularity as a method to prevent confounding [42]. Both of these examples have become standard methodologies in many biomedical analyses.
Conclusions
Cooccurrence and mutual exclusivity of somatic alterations are helpful concepts for the interpretation of cancer genomics data. For example, hypotheses about functional interactions between genes are often supported by suggested cooccurrence or mutual exclusivity of their alterations. Alarmingly, we have found that the statistical tests most commonly used for this purpose are not appropriate for testing the significance of cooccurrence. Many gene pairs that are believed to be coaltered more often than expected by chance do not exceed this expectation if the confounding effect of tumorspecific alteration rates is taken into account. Hypotheses formulated based on the results of those tests will therefore have limited support from the data. For this reason, we discourage the use of Fisher’s exact test or simple permutation methods for detecting cooccurrence. We have presented DISCOVER as a better alternative. Mutual exclusivity analysis using existing tests does not suffer from high false positive rates, but the sensitivity is low. DISCOVER identifies more significant mutual exclusivities without increasing the false positive rate. Thus, for both cooccurrence and mutual exclusivity analyses, we expect future cancer genomics studies to benefit from DISCOVER.
Methods
Independence statistic
We assess both cooccurrence and mutual exclusivity by counting how many tumors have an alteration in both genes and comparing this to the number of tumors expected to have such an overlap by chance if these alterations were independent. Importantly, the overlap expected by chance should factor in the fact that tumors with many alterations have a higher chance of such overlap than tumors with fewer alterations. Our null distribution modeling this overlap therefore takes into account both the alteration rate per gene and the alteration rate per tumor. To this end, let p _{ ij } denote the probability of an alteration in gene i and tumor j. We assume that the alteration probability of a gene is higher in tumors with many alterations overall than in tumors with fewer alterations. Therefore, p _{ ij } may be different from p _{ ik } for the same gene i in two different tumors j and k. Then, for two independent genes with alteration probabilities p _{1j } and p _{2j }, the probability of an alteration in both genes in tumor j is p _{1j } p _{2j }, while for tumor k it is p _{1k } p _{2k }. Given such probabilities for a set of tumors, the number of tumors that have an alteration in both genes follows a Poissonbinomial distribution.
The Poissonbinomial distribution [43] describes the sum of independent, nonidentically distributed Bernoulli random variables that have success probabilities p _{1},p _{2},…,p _{ n }. Its probability mass function is defined as follows:
Here, \(\mathcal {F}_{x}\) contains all subsets of size x of {1,2,…,n}, and A ^{c} denotes the complement of A.
Based on this distribution, we can estimate the probability of observing a number of tumors with alterations in two genes as extreme—as high for cooccurrence, or as low for mutual exclusivity—as the one observed.
If, for a given gene i, all probabilities p _{ ij } are equal for every tumor j, then the Poissonbinomial distribution reduces to a binomial distribution. However, estimating an individual alteration probability for every single tumor ensures that the heterogeneity in alteration rates across tumors is taken into account.
Estimating gene and tumorspecific alteration probabilities
To apply the DISCOVER test, we need estimates of the alteration probabilities p _{ ij } for all genes i and all tumors j. Let \(\mathcal {X} \in \{0, 1\}^{n \times m}\) denote the n×m binary alteration matrix where an entry x _{ ij } is 1 in case of an alteration in gene i and tumor j, and 0 otherwise. We use the notation x _{ i∙} and x _{∙j } for the marginal sums of the ith row and jth column, respectively. Furthermore, let X _{ ij } denote the random variable for x _{ ij }, and X _{ i∙} and X _{∙j } the corresponding marginal sums. If we were to assume that the alteration of a gene is equally likely across all tumors, then the alteration probability only depends on the number of altered tumors x _{ i∙} and the total number of tumors m:
Estimating the alteration probabilities this way ensures that the expected number of alterations \(\mathrm {E}_{p}(X_{i\bullet }) = \sum _{j} p_{ij}\) for a gene matches the observed number x _{ i∙}. In fact, the familiar expression above is the one that maximizes the likelihood of the observed alterations under the constraint that the expected number of alterations per gene matches the observed number. To make this more explicit, we can reformulate the probability estimation as a constrained optimization problem:
All of the above is based on the assumption that alteration probabilities for a gene are equal across tumors. Symptomatic for this assumption are probability estimates such that the expected number of alterations per tumor \(\mathrm {E}_{p}(X_{\bullet {}j}) = \sum _{i} p_{ij}\) generally does not match the observed number x _{∙j }. To take into account tumorspecific alteration rates, the above optimization problem can be extended such that this expectation is also matched:
With this new formulation, the number of parameters to fit is increased by a factor m. As a consequence, optimizing the likelihood \(\mathrm {L}_{p}(\mathcal {X})\) of the model risks overfitting the data. Therefore, instead of optimizing the likelihood, we choose to optimize the information entropy \(\mathrm {H}_{p}(\mathcal {X})\). It can be shown that in the optimal solution to this reformulated problem, each alteration probability can be written in terms of two parameters (Additional file 1: Parameter estimation):
Here, each parameter μ _{ i } for gene i is shared by all tumors, and each parameter λ _{ j } for tumor j is shared by all genes. Because of this, while the original optimization problem aims to estimate n×m alteration probabilities, we can obtain the optimal solution by estimating only n+m parameters. Moreover, all genes with the same number of altered tumors share the same value for μ _{ i }. Likewise, all tumors with the same number of altered genes share the same value for λ _{ j }. This sharing of parameters leads to an even larger reduction in the effective dimensionality of the optimization.
Unlike for the binomial case, there is no closedform solution for estimating the μ _{ i } and λ _{ j } parameters. Instead, we use the quasiNewton numerical optimization algorithm LBFGS [44].
Stratified analysis
When the data consist of clearly separate groups of tumors, such as is the case in the pancancer analysis with its different cancer types, it is preferable to stratify the analysis on these groups. For example, in the mutual exclusivity analysis, if group structure is not taken into account, the detected mutual exclusivities may be little more than markers for the underlying cancer types, rather than biologically related genes. An example of this type of confounding is presented in Additional file 1: Stratification in pancancer analysis. The DISCOVER test is easily stratified for different groups by solving the constrained optimization problem separately for the tumors of each group. The groupspecific background matrices can then be concatenated to construct a single global, but stratified, parameter matrix.
More formally, the binary alteration matrix \(\mathcal {X}\) can be seen as a concatenation of several n×m _{ c } submatrices \(\mathcal {X}_{c}\), where c∈{1,2,…,C} refers to one of C possible subgroups—e.g., a cancer type in the pancancer analysis—and m _{ c } is the number of tumors in that group:
To illustrate this, Additional file 1: Figure S6a shows an alteration matrix with tumors of two different subtypes. The parameter estimation procedure described in the previous section is then applied to each submatrix \(\mathcal {X}_{c}\) individually, resulting in subgroupspecific probability matrices P _{ c } (Additional file 1: Figure S6b). The global, stratified probability matrix is obtained by concatenating these matrices:
As in the nonstratified case, the expected number of alterations for each gene matches the observed number. However, unlike for the nonstratified probabilities, the expected numbers also match the observed numbers within each subgroup. With this stratified probability matrix, the Poissonbinomial test is applied in the same way as in the nonstratified setting (Additional file 1: Figure S6c).
False discovery rate control
Commonly used procedures for multiple testing correction assume that the P values are distributed uniformly under the null hypothesis. This is the case for, e.g., Bonferroni correction and the BenjaminiHochberg procedure. However, hypothesis tests that are based on a discrete test statistic, such as our DISCOVER test, are known to lead to nonuniform Pvalue distributions under the null hypothesis. In fact, pooling the P values across tests with a large set of different parameters results in a Pvalue distribution that is skewed towards 1.0. This complicates the application of the standard procedures for multiple testing correction. While these procedures would still control the familywise error rate or false discovery rate at the specified threshold, they will be more conservative because of the nonuniformity caused by the discrete test statistic. For the analyses in this paper, we used an adaptation of the BenjaminiHochberg procedure for discrete test statistics [45]. Further details on this procedure are provided in Additional file 1: False discovery rate control for discrete tests.
Groupbased mutual exclusivity test
We have defined a family of groupbased mutual exclusivity tests. The following statistics can be used to assess groupwise mutual exclusivity. Each of these statistics can be shown to follow a Poissonbinomial distribution, which we make use of to estimate significance.

Coverage: the number of tumors that have an alteration in at least one of the genes. Significance is based on the probability of observing a coverage at least as high in independent genes. The Poissonbinomial parameters for a group of genes {g _{ i }∣i∈I} can be derived from the individual gene alteration probabilities as follows:
$$ p_{j} = 1  \prod_{i \in I} (1  p_{ij}) \;, \qquad 1 \leq j \leq m $$That is, the probably of at least one alteration is one minus the probability of not having any alteration.

Exclusivity: the number of tumors that have an alteration in exactly one of the genes. Significance is based on the probability of observing exclusivity at least as high in independent genes. The Poissonbinomial parameters can be derived from the gene alteration probabilities as follows:
$$p_{j} = \sum_{i \in I} p_{ij} \prod_{k \in I\setminus{}\{i\}} (1  p_{kj}) \;, \qquad 1 \leq j \leq m $$ 
Impurity: the number of tumors that have an alteration in more than one gene. Significance is based on the probability of observing impurity at least as low in independent genes. The Poissonbinomial parameters can be derived from the gene alteration probabilities as follows:
$$\begin{aligned} p_{j} = 1 \! \prod_{i \in I} (1  p_{ij}) \! \sum_{i \in I} p_{ij} \prod_{k \in I\setminus{}\{i\}} \!(1  p_{kj}) \; \!, \qquad 1 \!\leq \!j\! \leq m \end{aligned} $$That is, the probability of more than one alteration is one minus the probabilities of no alterations and exactly one alteration. As a special case of this, if a group of only two genes is tested, the above expression reduces to p _{ j }=p _{1j } p _{2j }. This is the same parameterization as was used for the pairwise test.
Simulation data
An alteration matrix was constructed such that alteration frequencies across both genes and tumors resembled those of real tumors. For this, we used the copy number data of the TCGA breast cancer study as a reference. Based on the copy number matrix for 24,174 genes and 1044 tumors, we constructed two sequences of marginal counts corresponding to the number of amplifications across genes and across tumors. These two sequences were used as degree sequences to construct a random bipartite graph following the configuration model. The adjacency matrix of this bipartite graph was then used as the alteration matrix for the simulated data analyses. Because of the way this matrix was constructed, the alteration frequencies across both genes and tumors resemble those of the breast cancer tumors used for reference, yet there is no dependence between alterations across genes. For the analyses, only genes with at least 50 alterations were tested.
Mutually exclusive and cooccurring gene pairs, as well as mutually exclusive gene sets, were generated based on two parameters: coverage, the number of tumors altered in at least one of the genes; and impurity or overlap, the proportion of covered tumors altered in more than one of the genes. To generate pairs of mutually exclusive genes, we used quantile regression to relate the coverage of independent gene pairs to their impurity. Simulated mutually exclusive gene pairs were generated such that their impurity was below the first percentile predicted by the quantile regression model based on their coverage. Likewise, pairs of cooccurring genes were generated such that the number of tumors altered in both genes exceeded the 99th percentile based on the coverage of independent gene pairs.
Mutually exclusive gene sets were generated by first constructing sets of purely mutually exclusive gene alterations and then adding additional, nonexclusive alterations to obtain a prespecified degree of impurity. For the former, the percentage of covered tumors was randomly sampled from a truncated normal distribution with mean 0.4 and standard deviation 0.2, truncated on the interval [0.2, 0.8]. Next, individual gene alteration frequencies were sampled from the empirical distribution of alteration frequencies in the TCGA breast cancer matrix. Gene alteration frequencies were sampled until their sum reached the coverage of the group. The number of genes thus depends on the coverage in a way that is based on realistic cancer data. As some of the mutual exclusivity tests we compared with become intractable with larger numbers of genes, we restricted the maximum number of genes to 6. In addition, we also used a minimum gene set size of 3. Finally, the impurity was sampled from the set {0.02,0.05,0.08}. Impure alterations, i.e., additional alterations in an already covered tumor, were assigned to tumors with a probability proportional to the tumor’s overall alteration frequency.
For all analyses, the background matrix for the DISCOVER test was estimated on the complete alteration matrix, including genes with fewer than 50 alterations, and including simulated cooccurrences or mutual exclusivities.
Comparison to other mutual exclusivity tests
We compared the performance of the groupbased DISCOVER test to that of MEMo [6], muex [8], mutex [9], CoMEt [10], MEGSA [11], and TiMEx [12]. Some of these methods do more than just test for mutual exclusivity. They combine a statistical test for mutual exclusivity with an algorithm that identifies groups of genes to test. In our comparison, we were interested in comparing the performance of the statistical tests only. We therefore evaluated the mutual exclusivity tests by applying them to preidentified groups of genes.
For muex, MEGSA, and TiMEx, we used the R implementations provided with their respective publications. For CoMEt, we used a modified version of the official software implementation. Due to the computational complexity of the CoMEt test, it became intractable for some of the gene sets in the comparison. For this reason, the CoMEt publication suggests a set of heuristics to decide between the exact test and a faster binomial approximation, but we found those to be inadequate in our comparison. Instead, we changed the implementation such that it interrupts the CoMEt exact test after 1 minute and returns the P value obtained with the binomial approximation. For the MEMo and mutex tests, we used our own implementations, which we verified to give the same results as their original Java implementations.
Pancancer alteration data
Preprocessed somatic mutation and copy number data for the 12 cancer types studied in the TCGA pancancer initiative [22] were obtained via Firehose (analysis run 2014_07_15 at http://gdac.broadinstitute.org/runs/analyses__2014_07_15/). Mutations were extracted from the input of the MutSig 2CV analysis. Mutations for genes that have previously been identified as highconfidence mutational drivers [23] were included in the analysis. Discretized copy number changes were extracted from the output of GISTIC2. We considered genes altered if GISTIC2 qualified their copy number change as high level. Pancancer recurrently altered regions were obtained via Synapse (syn2203662 at https://www.synapse.org/#!Synapse:syn2203662). For each region, we selected their most likely driver genes for inclusion in the analysis. If a region contained only one gene, this gene was assumed its driver. In the case of more genes, genes were selected if they overlapped with the list of highconfidence mutational driver genes, or with a curated list of cancer genes (http://www.bushmanlab.org/links/genelists).
Background matrices for the DISCOVER test were estimated for each type of alteration—mutation, amplification, and deletion—separately, and based on the genomewide alteration matrices before gene selection. Stratification for the 12 different cancer types was applied as described before. The background matrix used in the analysis was subsequently composed from the relevant rows in the three alteration typespecific background matrices.
Overlap with the STRING functional interaction network
Version 10.0 of the STRING network [26] was used to determine overlap of detected mutual exclusivities and functional interactions. We constructed a functional interaction graph by connecting genes with an edge if they had a highconfidence STRING interaction, defined by a combined score greater than 800. A mutual exclusivity graph was constructed by connecting genes with an edge if alterations in these genes were found mutually exclusive at a maximum FDR of 1%. The overlap corresponds to the number of edges appearing in both graphs. To determine the enrichment of this overlap, we estimated a null distribution by randomly shuffling the gene labels of the mutual exclusivity graph 10,000 times and computing the overlap of these shuffled mutual exclusivity graphs with the unshuffled functional interaction graph.
De novo gene set detection
Our algorithm for detecting de novo sets of mutually exclusive genes combines two ideas from community detection. Its goal is to detect gene sets with a high likelihood of being mutually exclusive based on the results of a pairwise mutual exclusivity analysis. There are three main steps. First, a mutual exclusivity graph is constructed where genes are connected by an edge if their alterations have been identified as mutually exclusive by the pairwise test. For this step, we used a permissive significance criterion—a maximum FDR of 10%—so as not to exclude potentially interesting gene pairs that may simply not have reached significance due to the limited sample size. Second, groups of genes with a high density of mutual exclusivity edges between them are identified using a graph partitioning algorithm. Finally, these groups are subjected to the groupwise mutual exclusivity test to retain only those groups that are mutually exclusive as a group.
The graph partitioning step is based on overlapping correlation clustering. In correlation clustering, nodes in a graph are clustered such that the combined weight of edges within clusters is maximized and the combined weight of edges between clusters is minimized. The particular algorithm we used [46] allows nodes to be assigned to multiple clusters. Moreover, we modified the original algorithm such that groups of nodes can be designated that should always share the same cluster assignments. We used this for two situations. First, genes in the same copy number segment have highly correlated copy number alterations and, consequently, highly similar mutual exclusivities. Purely based on genomic data, there is no reason to prefer one gene over the other, which is why we always assign all such genes to the same cluster. Second, we assume that copy number alterations and mutations targeting the same gene serve the same function, and therefore add the constraint that these are always assigned to the same cluster.
The edge weights of the mutual exclusivity graph play an important role in the objective function of correlation clustering. A common phenomenon in pairwise associations is that one gene is found mutually exclusive with many other genes, but those genes are not all mutually exclusive with each other. The edges connecting the former gene may therefore not be indicative of gene set membership. They should be assigned a lower weight than edges that more specifically connect genes with a high degree of internal connectivity. To this aim, we selected the edge weights to optimize a modularity objective. In modularity optimization, a graph is compared with random graphs having the same number of nodes, edges, and degree distribution. Edges that are specific to the graph being partitioned are preferably kept within clusters, whereas edges that also appear in many of the random graphs will often span two clusters. We used a modularity measure based on conditional expected models [47]. This measure ensures that edges connecting sets of nodes with high node degrees receive a lower weight than edges that connect sets of nodes with low node degrees. It also allows for the covariance between the mutual exclusivity tests to be taken into account.
Abbreviations
 DISCOVER:

Discrete Independence Statistic Controlling for Observations with Varying Event Rates
 FDR:

False discovery rate
 i.i.d.:

Independent and identically distributed
 TCGA:

The Cancer Genome Atlas
References
 1
Vogelstein B, Papadopoulos N, Velculescu VE, Zhou S, Diaz Jr LA, Kinzler KW. Cancer genome landscapes. Science. 2013; 339(6127):1546–58. doi:10.1126/science.1235122.
 2
Lawrence MS, Stojanov P, Mermel CH, Robinson JT, Garraway LA, Golub TR, Meyerson M, Gabriel SB, Lander ES, Getz G. Discovery and saturation analysis of cancer genes across 21 tumour types. Nature. 2014; 505(7484):495–501. doi:10.1038/nature12912.
 3
Ashworth A, Lord CJ, ReisFilho JS. Genetic interactions in cancer progression and treatment. Cell. 2011; 145(1):30–8. doi:10.1016/j.cell.2011.03.020.
 4
Thomas RK, Baker AC, Debiasi RM, Winckler W, Laframboise T, Lin WM, Wang M, Feng W, Zander T, MacConaill L, Macconnaill LE, Lee JC, Nicoletti R, Hatton C, Goyette M, Girard L, Majmudar K, Ziaugra L, Wong KK, Gabriel S, Beroukhim R, Peyton M, Barretina J, Dutt A, Emery C, Greulich H, Shah K, Sasaki H, Gazdar A, Minna J, Armstrong SA, Mellinghoff IK, Hodi FS, Dranoff G, Mischel PS, Cloughesy TF, Nelson SF, Liau LM, Mertz K, Rubin MA, Moch H, Loda M, Catalona W, Fletcher J, Signoretti S, Kaye F, Anderson KC, Demetri GD, Dummer R, Wagner S, Herlyn M, Sellers WR, Meyerson M, Garraway LA. Highthroughput oncogene mutation profiling in human cancer. Nat Genet. 2007; 39(3):347–51. doi:10.1038/ng1975.
 5
Yeang CH, McCormick F, Levine A. Combinatorial patterns of somatic gene mutations in cancer. FASEB J. 2008; 22(8):2605–22. doi:10.1096/fj.08108985.
 6
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.
 7
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.
 8
Szczurek E, Beerenwinkel N. Modeling mutual exclusivity of cancer mutations. PLoS Comput Biol. 2014; 10(3):1003503. doi:10.1371/journal.pcbi.1003503.
 9
Babur Ö, Gönen M, Aksoy BA, Schultz N, Ciriello G, Sander C, Demir E. Systematic identification of cancer driving signaling pathways based on mutual exclusivity of genomic alterations. Genome Biol. 2015; 16:45. doi:10.1186/s1305901506126.
 10
Leiserson MD, Wu HT, Vandin F, Raphael BJ. CoMEt: a statistical approach to identify combinations of mutually exclusive alterations in cancer. Genome Biol. 2015; 16:160. doi:10.1186/s1305901507007.
 11
Hua X, Hyland PL, Huang J, Song L, Zhu B, Caporaso NE, Landi MT, Chatterjee N, Shi J. MEGSA: A powerful and flexible framework for analyzing mutual exclusivity of tumor mutations. Am J Hum Genet. 2016; 98(3):442–55. doi:10.1016/j.ajhg.2015.12.021.
 12
Constantinescu S, Szczurek E, Mohammadi P, Rahnenführer J, Beerenwinkel N. TiMEx: a waiting time model for mutually exclusive cancer alterations. Bioinformatics. 2016; 32(7):968–75. doi:10.1093/bioinformatics/btv400.
 13
Bredel M, Scholtens DM, Harsh GR, Bredel C, Chandler JP, Renfrow JJ, Yadav AK, Vogel H, Scheck AC, Tibshirani R, Sikic BI. A network model of a cooperative genetic landscape in brain tumors. JAMA. 2009; 302(3):261–75. doi:10.1001/jama.2009.997.
 14
Gorringe KL, George J, Anglesio MS, Ramakrishna M, Etemadmoghadam D, Cowin P, Sridhar A, Williams LH, Boyle SE, Yanaihara N, Okamoto A, Urashima M, Smyth GK, Campbell IG, Bowtell DD, Australian Ovarian Cancer Study. Copy number analysis identifies novel interactions between genomic loci in ovarian cancer. PLoS ONE. 2010; 5(9). doi:10.1371/journal.pone.0011408.
 15
Klijn C, Bot J, Adams DJ, Reinders M, Wessels L, Jonkers J. Identification of networks of cooccurring, tumorrelated DNA copy number changes using a genomewide scoring approach. PLoS Comput Biol. 2010; 6(1):1000631. doi:10.1371/journal.pcbi.1000631.
 16
Milosevic JD, Puda A, Malcovati L, Berg T, Hofbauer M, Stukalov A, Klampfl T, Harutyunyan AS, Gisslinger H, Gisslinger B, Burjanivova T, Rumi E, Pietra D, Elena C, Vannucchi AM, Doubek M, Dvorakova D, Robesova B, Wieser R, Koller E, Suvajdzic N, Tomin D, Tosic N, Colinge J, Racil Z, Steurer M, Pavlovic S, Cazzola M, Kralovics R. Clinical significance of genetic aberrations in secondary acute myeloid leukemia. Am J Hematol. 2012; 87(11):1010–6. doi:10.1002/ajh.23309.
 17
Kandoth C, McLellan MD, Vandin F, Ye K, Niu B, Lu C, Xie M, Zhang Q, McMichael JF, Wyczalkowski MA, Leiserson MD, Miller CA, Welch JS, Walter MJ, Wendl MC, Ley TJ, Wilson RK, Raphael BJ, Ding L. Mutational landscape and significance across 12 major cancer types. Nature. 2013; 502(7471):333–9. doi:10.1038/nature12634.
 18
Remy E, Rebouissou S, Chaouiya C, Zinovyev A, Radvanyi F, Calzone L. A modeling approach to explain mutually exclusive and cooccurring genetic alterations in bladder tumorigenesis. Cancer Res. 2015; 75(19):4042–52. doi:10.1158/00085472.CAN150602.
 19
Pereira B, Chin SF, Rueda OM, Vollan HK, Provenzano E, Bardwell HA, Pugh M, Jones L, Russell R, Sammut SJ, Tsui DW, Liu B, Dawson SJ, Abraham J, Northen H, Peden JF, Mukherjee A, Turashvili G, Green AR, McKinney S, Oloumi A, Shah S, Rosenfeld N, Murphy L, Bentley DR, Ellis IO, Purushotham A, Pinder SE, BorresenDale AL, Earl HM, Pharoah PD, Ross MT, Aparicio S, Caldas C. The somatic mutation profiles of 2,433 breast cancers refine their genomic and transcriptomic landscapes. Nat Commun. 2016; 7:11479. doi:10.1038/ncomms11479.
 20
Zack TI, Schumacher SE, Carter SL, Cherniack AD, Saksena G, Tabak B, Lawrence MS, Zhang CZ, Wala J, Mermel CH, Sougnez C, Gabriel SB, Hernandez B, Shen H, Laird PW, Getz G, Meyerson M, Beroukhim R. Pancancer patterns of somatic copy number alteration. Nature Genet. 2013; 45(10):1134–40. doi:10.1038/ng.2760.
 21
Lancaster HO. Significance tests in discrete distributions. J Am Stat Assoc. 1961; 56(294):223–34. doi:10.2307/2282247.
 22
The Cancer Genome Atlas Research Network, Weinstein JN, Collisson EA, Mills GB, Shaw KR, Ozenberger BA, Ellrott K, Shmulevich I, Sander C, Stuart JM. The Cancer Genome Atlas PanCancer analysis project. Nat Genet. 2013; 45(10):1113–20. doi:10.1038/ng.2764.
 23
Tamborero D, GonzalezPerez A, PerezLlamas C, DeuPons J, Kandoth C, Reimand J, Lawrence MS, Getz G, Bader GD, Ding L, LopezBigas N. Comprehensive identification of mutational cancer driver genes across 12 tumor types. Sci Rep. 2013; 3:2650. doi:10.1038/srep02650.
 24
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledgebased approach for interpreting genomewide expression profiles. Proc Natl Acad Sci U S A. 2005; 102(43):15545–50. doi:10.1073/pnas.0506580102.
 25
Ulz P, Heitzer E, Speicher MR. Cooccurrence of MYC amplification and TP53 mutations in human cancer. Nat Genet. 2016; 48(2):104–6. doi:10.1038/ng.3468.
 26
Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, HuertaCepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, von Mering C. STRING v10: proteinprotein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015; 43(Database issue):447–52. doi:10.1093/nar/gku1003.
 27
RodriguezViciana P, Warne PH, Dhand R, Vanhaesebroeck B, Gout I, Fry MJ, Waterfield MD, Downward J. Phosphatidylinositol3OH kinase as a direct target of Ras. Nature. 1994; 370(6490):527–32. doi:10.1038/370527a0.
 28
Gong Y, Zack TI, Morris LG, Lin K, Hukkelhoven E, Raheja R, Tan IL, Turcan S, Veeriah S, Meng S, Viale A, Schumacher SE, Palmedo P, Beroukhim R, Chan TA. Pancancer genetic analysis identifies PARK2 as a master regulator of G1/S cyclins. Nat Genet. 2014; 46(6):588–94. doi:10.1038/ng.2981.
 29
Yan H, Tong J, Lin X, Han Q, Huang H. Effect of the WWOX gene on the regulation of the cell cycle and apoptosis in human ovarian cancer stem cells. Mol Med Rep. 2015; 12(2):1783–8. doi:10.3892/mmr.2015.3640.
 30
Sard L, Accornero P, Tornielli S, Delia D, Bunone G, Campiglio M, Colombo MP, Gramegna M, Croce CM, Pierotti MA, Sozzi G. The tumorsuppressor gene FHIT is involved in the regulation of apoptosis and in cell cycle control. Proc Natl Acad Sci U S A. 1999; 96(15):8489–92. doi:10.1073/pnas.96.15.8489.
 31
Wang D, Wang L, Zhou J, Pan J, Qian W, Fu J, Zhang G, Zhu Y, Liu C, Wang C, Jin Z, He Z, Wu J, Shi B. Reduced expression of PTPRD correlates with poor prognosis in gastric adenocarcinoma. PLoS ONE. 2014; 9(11):113754. doi:10.1371/journal.pone.0113754.
 32
Solomon DA, Kim JS, Cronin JC, Sibenaller Z, Ryken T, Rosenberg SA, Ressom H, Jean W, Bigner D, Yan H, Samuels Y, Waldman T. Mutational inactivation of PTPRD in glioblastoma multiforme and malignant melanoma. Cancer Res. 2008; 68(24):10300–6. doi:10.1158/00085472.CAN083272.
 33
Zarubin T, Han J. Activation and signaling of the p38 map kinase pathway. Cell Res. 2005; 15(1):11–8. doi:10.1038/sj.cr.7290257.
 34
Regard JB, Malhotra D, GvozdenovicJeremic J, Josey M, Chen M, Weinstein LS, Lu J, Shore EM, Kaplan FS, Yang Y. Activation of Hedgehog signaling by loss of GNAS causes heterotopic ossification. Nat Med. 2013; 19(11):1505–12. doi:10.1038/nm.3314.
 35
He X, Zhang L, Chen Y, Remke M, Shih D, Lu F, Wang H, Deng Y, Yu Y, Xia Y, Wu X, Ramaswamy V, Hu T, Wang F, Zhou W, Burns DK, Kim SH, Kool M, Pfister SM, Weinstein LS, Pomeroy SL, Gilbertson RJ, Rubin JB, Hou Y, WechslerReya R, Taylor MD, Lu QR. The G protein α subunit G αs is a tumor suppressor in Sonic hedgehogdriven medulloblastoma. Nat Med. 2014; 20(9):1035–42. doi:10.1038/nm.3666.
 36
Takabatake Y, Takabatake T, Sasagawa S, Takeshima K. Conserved expression control and shared activity between cognate Tbox genes Tbx2 and Tbx3 in connection with Sonic hedgehog signaling during Xenopus eye development. Dev Growth Differ. 2002; 44(4):257–71. doi:10.1046/j.1440169X.2002.00640.x.
 37
Kann M, Bae E, Lenz MO, Li L, Trannguyen B, Schumacher VA, Taglienti ME, Bordeianou L, Hartwig S, Rinschen MM, Schermer B, Benzing T, Fan CM, Kreidberg JA. WT1 targets Gas1 to maintain nephron progenitor cells by modulating FGF signals. Development. 2015; 142(7):1254–66. doi:10.1242/dev.119735.
 38
Jagani Z, MoraBlanco EL, Sansam CG, McKenna ES, Wilson B, Chen D, Klekota J, Tamayo P, Nguyen PT, Tolstorukov M, Park PJ, Cho YJ, Hsiao K, Buonamici S, Pomeroy SL, Mesirov JP, Ruffner H, Bouwmeester T, Luchansky SJ, Murtie J, Kelleher JF, Warmuth M, Sellers WR, Roberts CW, Dorsch M. Loss of the tumor suppressor Snf5 leads to aberrant activation of the HedgehogGli pathway. Nat Med. 2010; 16(12):1429–33. doi:10.1038/nm.2251.
 39
Gross AM, Orosco RK, Shen JP, Egloff AM, Carter H, Hofree M, Choueiri M, Coffey CS, Lippman SM, Hayes DN, Cohen EE, Grandis JR, Nguyen QT, Ideker T. Multitiered genomic analysis of head and neck cancer ties TP53 mutation to 3p loss. Nat Genet. 2014; 46(9):939–43. doi:10.1038/ng.3051.
 40
MalinowskaOzdowy K, Frech C, Schönegger A, Eckert C, Cazzaniga G, Stanulla M, zur Stadt U, Mecklenbräuker A, Schuster M, Kneidinger D, von Stackelberg A, Locatelli F, Schrappe M, Horstmann MA, Attarbaschi A, Bock C, Mann G, Haas OA, PanzerGrümayer R. KRAS and CREBBP mutations: a relapselinked malicious liaison in childhood high hyperdiploid acute lymphoblastic leukemia. Leukemia. 2015; 29(8):1656–67. doi:10.1038/leu.2015.107.
 41
Leek JT, Storey JD. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genet. 2007; 3(9):161. doi:10.1371/journal.pgen.0030161.
 42
Yu J, Pressoir G, Briggs WH, Vroh Bi I, Yamasaki M, Doebley JF, McMullen MD, Gaut BS, Nielsen DM, Holland JB, Kresovich S, Buckler ES. A unified mixedmodel method for association mapping that accounts for multiple levels of relatedness. Nat Genet. 2006; 38(2):203–8. doi:10.1038/ng1702.
 43
Wang YH. On the number of successes in independent trials. Statistica Sinica. 1993; 3(2):295–312.
 44
Liu DC, Nocedal J. On the limited memory BFGS method for large scale optimization. Math Prog. 1989; 45(13):503–28. doi:10.1007/BF01589116.
 45
Carlson JM, Heckerman D, Shani G. Estimating false discovery rates for contingency tables. Technical Report MSRTR200953. 2009. http://research.microsoft.com/apps/pubs/default.aspx?id=80571.
 46
Bonchi F, Gionis A, Ukkonen A. Overlapping correlation clustering. Knowl Inf Syst. 2013; 35(1):1–32. doi:10.1007/s1011501205229.
 47
Chang YT, Leahy RM, Pantazis D. Modularitybased graph partitioning using conditional expected models. Phys Rev E. 2012; 85(1):016109. doi:10.1103/PhysRevE.85.016109.
Acknowledgements
We thank Daniel Vis and Evert Bosdriesz for their helpful comments on the manuscript, and John Foekens and Marcel Smid for useful discussions. The results published here are in part based upon data generated by the TCGA Research Network.
Funding
This study was funded by Cancer Genomics Netherlands.
Availability of data and materials
Python and R packages implementing DISCOVER are available at https://github.com/NKICCB/DISCOVER, and are licensed under the Apache License 2.0. The version of the DISCOVER software used for this manuscript has been assigned doi:10.5281/zenodo.167332. Jupyter notebooks containing all the code required to reproduce the figures and results of this manuscript are available at https://github.com/scanisius/discovernotebooks. These have been assigned doi:10.5281/zenodo.167469.
Authors’ contributions
SC, JM, and LW conceived the study. SC developed and implemented the method and performed the analyses. LW supervised the study. SC and LW wrote the manuscript. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Author information
Affiliations
Corresponding author
Additional files
Additional file 1
Supplementary methods and Figures S1–S6. (PDF 897 kb)
Additional file 2
Table S1. Significantly mutually exclusive alterations found in the pancancer data at a maximum FDR of 1%. (XLSX 16 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Canisius, S., Martens, J. & Wessels, L. A novel independence test for somatic alterations in cancer shows that biology drives mutual exclusivity but chance explains most cooccurrence. Genome Biol 17, 261 (2016). https://doi.org/10.1186/s130590161114x
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s130590161114x
Keywords
 Mutual exclusivity
 Cooccurrence
 Computational biology