Identifying restrictions in the order of accumulation of mutations during tumor progression: effects of passengers, evolutionary models, and sampling
- Ramon Diaz-Uriarte^{1}Email author
https://doi.org/10.1186/s12859-015-0466-7
© Diaz-Uriarte; licensee BioMed Central. 2015
Received: 15 July 2014
Accepted: 15 January 2015
Published: 12 February 2015
Abstract
Background
Cancer progression is caused by the sequential accumulation of mutations, but not all orders of accumulation are equally likely. When the fixation of some mutations depends on the presence of previous ones, identifying restrictions in the order of accumulation of mutations can lead to the discovery of therapeutic targets and diagnostic markers. The purpose of this study is to conduct a comprehensive comparison of the performance of all available methods to identify these restrictions from cross-sectional data. I used simulated data sets (where the true restrictions are known) but, in contrast to previous work, I embedded restrictions within evolutionary models of tumor progression that included passengers (mutations not responsible for the development of cancer, known to be very common). This allowed me to assess, for the first time, the effects of having to filter out passengers, of sampling schemes (when, how, and how many samples), and of deviations from order restrictions.
Results
Poor choices of method, filtering, and sampling lead to large errors in all performance measures. Having to filter passengers lead to decreased performance, especially because true restrictions were missed. Overall, the best method for identifying order restrictions were Oncogenetic Trees, a fast and easy to use method that, although unable to recover dependencies of mutations on more than one mutation, showed good performance in most scenarios, superior to Conjunctive Bayesian Networks and Progression Networks. Single cell sampling provided no advantage, but sampling in the final stages of the disease vs. sampling at different stages had severe effects. Evolutionary model and deviations from order restrictions had major, and sometimes counterintuitive, interactions with other factors that affected performance.
Conclusions
This paper provides practical recommendations for using these methods with experimental data. It also identifies key areas of future methodological work and, in particular, it shows that it is both possible and necessary to embed assumptions about order restrictions and the nature of driver status within evolutionary models of cancer progression to evaluate the performance of inferential approaches.
Keywords
Background
Cancer progression is caused by the sequential accumulation of somatic mutations, including changes in copy number (structural variants), single nucleotides (SNP variants) and DNA methylation patterns during the life of an individual [1-3]. Among the mutations causally responsible for the development of cancer (drivers) not all possible orders of accumulation seem equally likely, and the fixation of some mutations can depend on the presence of other mutations. For example, in colorectal cancer APC mutations are an early event that precedes mutations in KRAS [4-6]. Understanding the restrictions in the temporal order of accumulation of driver mutations not only provides insights into cancer biology, but can help identify early markers of disease as well as therapeutic targets [5-9], and can be an instrumental tool in the search for the “Achilles’ Heel” of oncogene addiction [3,10,11]. In addition, understanding the correct order of events is necessary for the assessment of the validity of the genetic context of cell lines and animal models of human cancer [7,8].
In this context, a variety of methods have been developed to try to infer the possible restrictions in the order of accumulation of driver mutations from cross-sectional data. Longitudinal data would be better suited for this problem but it is much harder to obtain and cross-sectional data is (and will remain) the main source of data (e.g., the growing number of genomes available through international sequencing projects) for addressing these and similar problems [5,12]. I provide next a brief review of the main methods, including recent developments, but see more extensive reviews in [13,14]. The oncogenetic tree (OT) model [15] was introduced as an extension of the linear path model [16]: in OTs progression starts from a common (non-altered) root, and branches out, so that there are several mutational pathways that can be observed simultaneously. OTs, by virtue of being trees, can only model order restrictions where an event depends on its single parent. Another early model are distance-based trees [17,18], but their meaning is rather different, since the observed mutations are only placed in the leaves or terminal nodes of the tree, and the internal nodes are unobserved and unknown events, which precludes an interpretation in terms of order restrictions like “mutation A is required for mutation B”. Distance-based trees and other models [19] that do not try to infer order restrictions will not be considered further in this paper.
Conjunctive Bayesian Networks (CBNs) [20] were developed as a generalization of OTs: these are graphs where the occurrence of a mutation can depend on the occurrence of two or more parents (i.e., a conjunction). The disease progression models of OTs and CBNs assume that a mutation can only occur with non-negligible probability if the preceding parent mutation(s) in the graph have occurred, which has been called monotonicity [12]. Thus, for driver genes, under strict OT and CBN models it would be impossible to observe a genotype that is not compatible with the relations specified in the graph. Less restrictive models for tumor progression were suggested early on, including general Markov models and Bayesian Networks which allow for mutations to occur even if no other aberrations have occurred [14,21-23]. Progression Networks [12] have been proposed for learning models that include OTs, CBNs, as well as several other special types of Bayesian Networks, and can explicitly incorporate deviations from monotonicity. Retracing the Evolutionary Steps in Cancer (RESIC) [7,8] differs from other methods because it attempts to find the order of events taking into account the evolutionary dynamics of mutation accumulation. CBNs, OTs, and Progression Networks can be directly applied to module/pathway data, provided those data are partitioned into predefined pathways before the analysis (e.g., [6,7]), although some recent work [5,22], not the focus of this paper, simultaneously tries to find modules or pathways and their order restrictions.
Having a single graph means having a single set of restrictions that is common to all individuals, but that does not mean that all cells follow the same path (so the actual genotypes and their paths can be quite diverse under one graph). Mixtures of OTs [24] and mixtures of Hidden-variable OTs [25] are a further generalization of OTs where disease progression is modeled allowing for different order restrictions in different subsets of individuals, each one modeled as a (Hidden-variable) OT. By using a star as one of the trees in the mixture, these models can also account for any mutation occurring without its parent(s) having occurred. In this paper I restrict attention to finding a single graph, the approach most widely used in the literature (but see Discussion).
Most of the above are general methods, and can be applied to different kinds of data including cytogenetic, gene mutation, and pathway alteration data [6,15]. This versatility, coupled with the increasing wealth of cross-sectional data available, provides an excellent opportunity to try to understand the still largely unknown details of the order of mutations. However, in spite of the relevance of the problem for both diagnostic and therapeutic purposes, there are very few systematic comparisons of method performance, and they do not provide a clear and robust answer to the question of method choice.
Applied usage of the above methods faces at least three additional major problems. First, most of the mutations present in cancer cells are not driver mutations, but passenger mutations not responsible for the development of cancer [26-29]. Passenger mutations can show a non-negligible frequency because they “hitchhike” on drivers [1,30]. Unless we know what mutations are drivers, the presence of passengers in our data sets forces us to use some filtering procedure to select which mutations (or, generally, alterations) to use with (or to pass on to) the methods to infer order restrictions. However, the simulations in the only comparison of methods available [13], as well as in the original descriptions of new methods [12,31,32], have all been conducted assuming that the identity of the driver mutations is known.
Virtually all papers that try to infer order restrictions, including methodological papers, rely on simple frequency-based selection or filtering procedures to select which genes to use [6,15,18,22,32-36], but the effects of these filtering approaches on the performance of the methods to infer order restrictions are completely unknown.
Second, attention to sampling decisions is largely missing from the literature. OTs, CBNs, and Progression Networks are generative models, and simulations that examine method performance [12,13,31,32] obtain genotypes directly from these generative models. But, except when we use single cell sampling, our experimental data are from samples that aggregate over many cells and the joint and marginal frequencies of mutations of those aggregates can depend not only on the aggregation per se but also on when we sample (due, for example, to the clonal expansion episodes), and differ greatly from distributions obtained from the generative models.
Finally, development and evaluation of methods of reconstruction of order restrictions are conducted without consideration for the evolutionary model of tumor progression (but see [7,8] and Discussion). This problem is highlighted by Sprouffske et al. [37]: referring to oncogenetic tree models they say (p. 1136) “This is not an evolutionary model because the oncogenetic tree does not represent ancestral relationships within a neoplasm but rather a summary of the observed co-occurrences of mutations across independent neoplasms”. This lack of consideration for the evolutionary model is also unfortunate since it does not provide a clear mechanistic interpretation of (nor a simple mechanistically-based procedure for generating) deviations from the restrictions encoded in the graph. Of particular interest is monotonicity (a mutation in a driver gene can only be observed if the preceding parent mutations in the graph have occurred), because deviations from it can easily arise when a mutation behaves as a driver or as a passenger depending on the genetic context —i.e., depending on which other genes are mutated [30,38].
Factors considered and their levels or possible values, together with acronyms used through the text
Factor | Description | Values |
---|---|---|
Model | Evolutionary model of cancer progression | exp, Bozic, McF_4, McF_6 |
sh | Penalization of deviations from monotonicity | 0, Inf (for ∞) |
True graph | The true graph: the structure that encodes the order restrictions. All possible combinations of Number of nodes and Conjunction | 11-A, 11-B, 9-A, 9-B, 7-A, 7-B |
Number of nodes (NumNodes) | Number of genes or alterations | 11, 9, 7 |
Conjunction | Whether or not the graph has conjunctions | Yes, No |
Sample size (S.Size) | Number of samples used for reconstructing the graph | 100, 200, 1000 |
Sampling time (S.Time) | When the sample is taken | Last, unif (for uniform) |
Sampling type (S.Type) | How tissue is collected | singleC (for single cell), wholeT_0.5 (whole tumor, detection threshold=0.5), wholeT_0.01 (whole tumor, detection threshold=0.1) |
Filtering | Method for selecting drivers, or filtering passengers, when the true drivers are not known | S1, S5, J1, J5 (for frequency of Single event and Joint frequency of events, with thresholds 1% and 5% respectively) |
Method | Method for inferring the order restrictions | CBN, CBN-A, DiP, DiP-A, OT, OT-A |
In this paper I incorporate the order restrictions into evolutionary models to address how the performance of all available methods for inferring order restrictions is affected by: a) passenger mutations that lead to uncertainty about the identity of the true drivers and the need to use filtering approaches; b) sampling choices (when and how and how many to sample); c) type of underlying true graph, including presence/absence of conjunctions; d) deviations from the order restrictions encoded in the graphs (deviations from monotonicity); e) evolutionary models of tumor progression.
Methods
Table 1 provides an overview of the main factors considered in this study, and Figure 1 a schema of all the steps. We will deal separately with two different scenarios: one where we know the true identity of the drivers, which we will call “Drivers Known”, and another scenario that replicates the common situation where the data includes both passengers and drivers and we do not know with certainty which is which, a scenario we will refer to as “Drivers Unknown”. The “Filtering” factor is only relevant for the “Drivers Unknown” scenario (as shown in Figure 1 by the “No” path from the decision diamond “Drivers known?”).
Details about the experimental design are provided in Additional file 1. Briefly (see also Figure 1) each simulation run produces the observed genotype for a subject, and a data set is made from the genotypes of multiple subjects. We can analyze the same data set with different Methods (or different Filtering by Method combinations) but for other factors (e.g., Model) different settings of the factor produce different data sets. I have used use an experimental design with among- and within-data set factors (see Table 1) so as to examine the effect of Method and Filtering, controlling for possible among-data set variation.
For factors Model, sh (see details below in section “Deviations from monotonicity and genetic context dependence of driver status: sh ”), True Graph (= Number of Nodes * Conjunction), S.Size, S.Type, and S.Time, the among-data set factors, I used a full factorial design (thus, 4∗2∗3∗2∗3∗3∗2=864 among-data set factor combinations). For every combination of the among-data set factors I used twenty independent replicate data sets. Each of the twenty replicate data sets was analyzed with every Method or every Filtering by Method combination (the within-data set factors) to infer a graph from the data (i.e., to try to infer the order restrictions among events). Therefore, a total of 864∗6=5184 or 864∗6∗4=20736 factor combinations for the Drivers Known and Drivers Unknown scenarios, respectively, were examined.
Evolutionary models and simulation
Main parameters for each of the tumor progression models
Model | Birth rate | Death rate | Mutation | Cancer |
---|---|---|---|---|
rate (per | reached if | |||
gene per | ||||
unit time) | ||||
Bozic | 1 | (1−s)^{ j }(1+s _{ h })^{ p } | 10^{−6} | >10^{9} cells |
exp | (1+s)^{ j }(1−s _{ h })^{ p(+)} | 1 | b _{ j }∗10^{−7} | >10^{9} cells |
McF_4 | (1+s)^{ j }/(1+s _{ h })^{ p } | log(1+N/K) | 5∗10^{−7} | Number of |
drivers ≥4 | ||||
McF_6 | (1+s)^{ j }/(1+s _{ h })^{ p } | log(1+N/K) | 5∗10^{−7} | Number of |
drivers ≥6 |
Details about the models, choice of parameters, simulations, and examples of simulated trajectories are provided in Additional file 1.
Deviations from monotonicity and genetic context dependence of driver status: sh
Genetic context dependence of driver/passenger status [30,38] and deviations from monotonicity (i.e., from the order of events implied by the graph of the oncogenetic model) can be closely related, and affect the performance of methods to infer order restrictions. A mutation in a driver gene for which all the preceding required mutations have occurred (i.e., a mutation in a gene that has its dependencies satisfied) will lead to an increase in fitness (through its increase of the value j, for number of drivers, as shown in Table 2). What about mutations in driver genes that do not have their dependencies satisfied? Enforcing monotonicity is equivalent to considering such a mutation as a mutation in an essential housekeeping gene, which can be modeled setting s _{ h }, in the notation of [44], to ∞ (so fitness of such clones is zero). Deviations from monotonicity can arise, however, if such mutation is similar to a passenger mutation: it confers no fitness benefits (and if it has no deleterious effect s _{ h }=0, similar to setting s _{ p }=0 in [42]). Of course, in none of these two cases (restrictions not satisfied) would the value of j be increased (because a driver only increases fitness if its dependencies are satisfied). In the simulations reported here I considered two extreme scenarios: a) no deviations from the graph of the oncogenetic model are allowed, which I will refer to as “sh=Inf” (from ∞) ; b) drivers without dependencies satisfied are equivalent to passengers with no deleterious effects, which I will refer to as “sh=0”. Note that the implementations I used to infer CBN and OT incorporate errors [6,31,32] and the OT model explicitly allows for errors due to the occurrence of genetic events outside the model implied by the graph of the oncogenetic model [45]. DiProg (DiP), the method to infer Progression Networks models [12], explicitly incorporates deviations from monotonicity with the parameter ε.
True graphs, number of drivers, and number of passengers
Six different true graphs have been used, three of them with conjunctions (i.e., graphs that could only be perfectly inferred with either CBN or Progression Networks) and three of them without conjunctions (i.e., trees that could be perfectly reconstructed by all methods compared). The trees (the graphs without conjunctions) are derived from the graphs with conjunctions by removing conjunctions. The graphs have 7, 9, and 11 nodes. The number of nodes of the graphs targets the range of nodes commonly considered in studies that try to reconstruct graphs from real cancer data: 10, 11, and 12 in [31], 7 to 11 (including gene and core pathways) in [6], 7 in [15], 11 and 12 in [33], 8 in [46], 7 (modules) in [22], 9 in [47], 17 in [9] and [23], 12 in [48], 6 and 13 in [21], 12 in [14]. The size of graphs was limited to 11 because CBN cannot deal with more than 14 nodes and we need to allow for the possible selection of more than 11 nodes when Drivers are Unknown. The graphs are shown in Additional file 1, and we will refer to them by the number of events, using post-fix A for conjunction and B for no conjunction.
Simulations that used graphs with 11, 9, and 7 seven drivers generated clones with between two and six drivers (see also Additional file 1), a range which is well within the range of drivers considered in the literature: although some authors [44,49] examine scenarios with 20 drivers, most studies deal with much smaller numbers of drivers [37,41,42] and recent reviews suggest that the number of drivers in the cells of most tumors lies between two and six [26,42]. Regarding number of passengers, it is now widely accepted that most mutations in cancer cells are passengers [26-29,41,50]. In the Drivers Unknown scenario I set the proportion of passengers to drivers constant, so that there are four passengers for every driver, a range within that seen in the literature. Our scenario is also relevant if the actual number or fraction of passengers is much larger, but many of those passengers can be excluded a priori based on other information, so that they are never considered as candidates for the process of filtering data and inferring graphs (i.e., they are never passed on to step 3 in Figure 1).
Sample size and sampling type and time
Sample size (S.Size) was set to three possible values: 100, 200, 1000. The values considered in other studies vary widely (100, 400, and 800 in [33]; 83 to 95, plus a pool of 268 in [6]; 50, 100, 200, 500, and 1000 in [13]; 971 in [46]; 887 in [9]). The set of 100, 200, and 1000 covers a realistic range of sample sizes and will allow us to compare the effects of sample size with those of other factors. Sampling time (S.Time) refers to when sampling occurs. S.Time = last means that samples were collected at the end of the simulation (at the end of cancer progression). S.Time=unif (for uniform) means that sampling time was uniformly distributed between the time of appearance of the first mutated driver and the end of the simulation. Uniform sampling is a very simple model for obtaining cross-sectional samples of patients at different stages of the disease. Sampling uniformly between the time of appearance of the first mutated driver and the final stages of the disease is, of course, unrealistic, but I used this type of sampling because it provides a stark contrast with sampling at the end of the disease: S.Time=unif and S.Time=last can be regarded as two extremes of sampling tissue that harbors at least one mutated driver (i.e., cancerous or pre-cancerous tissue). Sampling type (S.Type) refers to whether single cell or whole tumor sampling was used, and three values have been used for this factor. When S.Type = singleC (single cell), a simulation provided the genotype of one single cell (or, equivalently, one single clone), where the probability of selecting a clone was proportional to its abundance. When using whole tumor sampling, and as in [37], a biopsy was the entire tumor, but whether a gene was considered mutated or not depended on the detection threshold, and here I used two levels: 0.5 (like [37]) and 0.01, meaning that a gene was considered mutated if it was mutated in 50% or 1%, respectively, of the cells. Of course, it is unlikely that a study using single cell would take a single sample from a patient but we focus on cross-sectional data, and single cell sampling is the type of sampling that leads to data most similar to the type of data obtained when we simulate using the generative model of the underlying graph. Moreover, single cell sampling and whole tumor sampling, as used here, can be considered two extremes in the range of sampling possibilities. Likewise, a detection threshold of 0.01 is probably unrealistically low but that setting is used because it combines the capacity of detecting very low frequency co-occurring events (as in single cell sampling) with summing over distinct cells (which could lead to problems similar to the ecological fallacy). The sampling schemes used here ignore any possible spatial structure and tissue architecture [51,52], not because they are considered irrelevant, but because none of the evolutionary models considered here incorporate them (but note that the uniform sampling scheme can sometimes be equivalent to incorporating spatial structure, if that spatial structure is correlated with time).
Filtering
When the identity of the driver genes is not known, it is often necessary to select genes before trying to infer the order restrictions. Some studies that deal with chromosomal abnormalities have used the methods of Brodeur and collaborators or Taetle and collaborators, to try to locate non-random breakpoints (see discussion and references in [14,21]) but these methods are not directly applicable to other types of data. Other authors that deal with chromosomal abnormalities, or that use mutation data, have used one of the following general approaches to decide which alterations to analyze: a) selecting the most frequent mutations, either by setting a minimal number as in [32,34], and [22], where the seven, 13, or 25, respectively, most frequent alterations are used, or setting a minimal frequency such as in [6,18,35,36] where the threshold is set at 5%, 5%, 10%, 10%, respectively; b) selecting the largest set of events so that every pair of events is observed at least k times, as in [45] and [15] where the threshold is five times, out of 124 and 117 cases, respectively. The key difference between these two filtering procedures is that the second uses the joint occurrences of pairs. To comprehensively incorporate common uses, I used four filtering procedures: two of them only consider the marginal frequency of each single event, and I use an “S” to denote “frequency of Single event”, and the other two take into account joint occurrences, and a I use “J” to denote “Joint frequency of pairs of events”. The procedures are S1, that selects any mutation with a frequency larger than 1%, S5 where the threshold is 5%, J1 that selects the largest set of events so that every pair of events is observed at least in 1% of the cases and J5, where the threshold is 5%. In the rare case where a filtering procedure returned more than 12 mutations, the 12 most common were selected.
Our focus here is in the effect that a filtering procedure has on the reconstruction of the order restrictions. Filtering, by itself, can introduce errors (true drivers can be missed or passengers can be considered drivers) but these errors can have different impact on the reconstructions of the order restrictions depending on the methods; colloquially, different methods are not necessarily equal when trying to make the best of a bad situation. An example and further details are provided in Additional file 2, “Commented example of filtering + method effects”.
Inferring order restrictions: CBN, DiP, OT
I have used three types of methods to infer order restrictions from data: methods that infer OTs, methods that infer CBNs (which should also be able to reconstruct OTs), and methods that infer Progression Networks (and, thus, should be able to reconstruct both OTs and CBNs). Each method, when applied to a data set, returns what we will refer to as an “inferred graph” (see Figure 1). For OTs I used the R package Oncotree [53] with its default settings. Some of the analysis were rerun (see Discussion and Additional file 3) with the implementation available in the BioConductor package Rtreemix [54,55]. For CBNs as detailed in [31] I used the software from [6]. I used the same default settings for temp (1) and steps (n u m b e r o f n o d e s ^{2}) and started the simulated annealing search for the best poset from an initial linear poset as in [6]. For Progression Networks I used the DiProg program (the method we call DiP) from [12] to fit monotone networks (option “MPN”), choosing the best k from 1 to 3 (and the results reported here have ε=0.05). Further details about software versions and parameters used for all methods are provided in Additional file 1.
Other methods have been described in the literature, but I have not been able to use them here. The method in [33] is too slow (analysis of data sets of 200 cases exceeding 4 hours; see further details in Additional file 2) if we need to do more than 100,000 analysis, as in this paper. The methods in [7,25] have no software available. Therefore, this paper includes all currently existing approaches for which software is available. It is worth emphasizing the crucial role of the availability of free and open source software both in the growth and development of bioinformatics and computational biology [56-58] and for implementing reproducible research [59]. Moreover, the lack of public implementations precludes comparison of otherwise promising approaches, which ultimately hurts practitioners [60].
Finally, it is important to mention that the methods used differed greatly in speed: the median and mean execution times, over all 172800 analysis performed by each family of methods, were 0.045 and 0.07 seconds for OT, 3.89 and 12.60 seconds for DiP, but 31 and 1127 seconds for CBN. In addition, DiProg (DiP) currently depends on IBM’s CPLEX ILOG library, which not only is not open source but has a severely restrictive license. Further details of execution times are provided in Additional file 2.
CBN-A, OT-A, DiP-A
In some cases one or more mutations were present in all or almost all of the subjects. Even if these are driver mutations on which all other events depend, events with a frequency of 1 are often removed from the graph (e.g., by the OT method) or placed as nodes that descend directly from Root and that have no descendants. To try to minimize this problem, we can augment the data by adding “pseudosamples” that have no mutations in any gene. Adding “pseudosamples” does not amount to knowing anything about the order of events, nor the truth about which genes are drivers or not (and in the Drivers Unknown scenario I always augmented after the filtering step). Data augmentation only requires being able to differentiate between presence and absence of a genetic alteration, mutation, or aberration, which is always assumed in these analyses. In this paper, “CBN-A”, “DiP-A”, and “OT-A”, refer to using CBN, DiP, or OT on data that has been augmented by adding to it another 10% of samples filled with zeroes (0 is the code that denotes no alteration).
Analysis
We want to address two questions: a) what procedures (choice of Method, Filtering, S.Time, and S.Type) are “best” (for reconstructing the underlying true graph from the data), so that we can choose a course of action when faced with new data; b) what factors have an important effect on performance, including interactions with other factors, even if they are not under user control, so that future research can focus on them. The first question (what method is best?) is most straightforwardly addressed by ranking Method(s) and Method by Filter combinations and by finding the best Method(s) (or Method by Filtering combination(s)) using the Multiple Comparisons with the Best procedure.
The second question (what factors affect performance? are there interactions among them?), is best addressed with statistical modeling that focuses on identifying factors with relevant effects; for instance, here the question would not be whether the Bozic model or the McF_4 model lead to better performance, but rather whether Evolutionary Model affects performance and shows interactions with other factors. The approaches used reflect these two questions and are based on very different procedures and assumptions. Of course, results from the different approaches complement each other (see further comments in section “Why GLMs, MCB, and ranking?” in Additional file 1). Below I detail the different analyses, after explaining how performance was measured.
Performance measures
I consider here that the main goal of most studies is the reconstruction of the topology of the graph, which is what captures the order restrictions [12,15,32]. There is no single performance measure that can fully characterize the performance in this task, and therefore I have used four performance measures that capture performance along different dimensions. One is a global score of the difference between the inferred graph and the true graph. The other three are measures of classification or diagnostic performance common in medical testing and machine learning [61,62] that focus on the fractions or proportions across specific rows or columns of the confusion matrix (where entries in that matrix are commonly called “true positives”, “false positives”, “false negatives”, and “true negatives”). Thus, the dimensions measured by each of these four performance measures relate to concepts already familiar to researchers, and arguably capture the key features of the methods’ behavior. As we will see below, using these four different performance measures is also key to understanding some of the major differences between methods.
Diff is the sum of the absolute value of the entries in the matrix of the Difference between the adjacency matrices of the true (T) and inferred (F) graphs; this is the square of the “usual” Frobenius norm [63] of that matrix difference, and is the same as the “graph edit distance” of [13]. The Proportion of False Discoveries (PFD) is defined as \(\frac {{\# of relations in } F { but not in}T}{{\# of relations in } F}\). Following [6,31], we define “relations” as the transitive closure of “cover relations”. For instance, suppose a graph with A→B→C (where A→B means that A needs to occur for B to occur); the cover relations are A→B and B→C, but we also include A→C in the relation. As in [6,31], we do not include the root node when finding cover relations and their transitive closure (in contrast to what is done in the computation of Diff). The numerator is, therefore, the number of false positives (FP). The Proportion of Negative Discoveries (PND) is defined as \(\frac {{\# of relations in } T { but not in } F}{{\# of relations in } T}\) and the numerator is, therefore, the number of false negatives (FN). The False Positive Fraction (FPF) is defined as \(\frac {{\# of relations in } F { but not in } T}{{\# of relations not in } T}\) and its numerator is the number of FP; it should be noted, however, that the FPF is of minor value compared to PFD and PND. Further details for all performance measures are provided in Additional file 1.
Overall ranking of filtering, method, and sampling scheme
Ranking of all 36 combinations of Method and Sampling scheme (time, type) when drivers are known with respect to each performance measure
Method and sampling | Conjunction | No conjunction | |||||||
---|---|---|---|---|---|---|---|---|---|
Diff | PFD | PND | FPF | Diff | PFD | PND | FPF | ||
OT-A, last, singleC | 1 | 2 | 14 | 10 | 1 | 3 | 2 | 15 | |
OT-A, last, wholeT_0.5 | 2 | 1 | 15 | 12 | 2 | 4 | 3 | 14 | |
OT-A, last, wholeT_0.01 | 3 | 3 | 7 | 22 | 3 | 6 | 1 | 22 | |
OT-A, unif, singleC | 4 | 6 | 19 | 13 | 4 | 9 | 12.5 | 9.5 | |
OT, unif, singleC | 5 | 5 | 20 | 11 | 5 | 7 | 12.5 | 9.5 | |
OT-A, unif, wholeT_0.01 | 8 | 11 | 16 | 24 | 8 | 11 | 5 | 24 | |
OT, last, singleC | 10 | 9 | 23 | 1 | 10 | 1 | 17 | 2 | |
OT, last, wholeT_0.01 | 11 | 4 | 18 | 18 | 12 | 5 | 14 | 18 | |
OT, last, wholeT_0.5 | 12 | 7 | 24 | 2 | 11 | 2 | 23 | 1 | |
CBN-A, unif, wholeT_0.01 | 13 | 13 | 1 | 26 | 13 | 13 | 4 | 26 | |
CBN-A, unif, singleC | 14 | 16 | 2 | 28 | 15 | 15 | 9 | 29 | |
CBN-A, unif, wholeT_0.5 | 15 | 18 | 3 | 34 | 14 | 16 | 8 | 31 | |
CBN, unif, singleC | 16 | 17 | 5 | 29 | 17 | 19 | 10 | 34 | |
CBN, unif, wholeT_0.01 | 17 | 14 | 4 | 27 | 18 | 14 | 6 | 27 | |
DiP-A, unif, singleC | 31 | 28 | 31 | 4 | 31 | 30 | 31 | 6 | |
DiP, last, wholeT_0.5 | 33 | 35 | 34 | 5 | 34 | 34 | 36 | 5 | |
DiP, unif, singleC | 35 | 31 | 33 | 3 | 35 | 32 | 33 | 3 | |
DiP, unif, wholeT_0.5 | 36 | 36 | 36 | 6 | 36 | 36 | 35 | 4 |
This procedure does not take into account the repeated usage of the same data set for each of the four/sixteen methods, but we are using it simply to rank alternatives. The advantage of this procedure is that it provides an overall view of the results that is equivalent to examining all possible interactions of Filtering by Method by Sampling scheme, marginalizing over all other terms, and this is done with a simple procedure that does not depend on additional modeling assumptions. The disadvantage is that it does not allow us to judge the relative size of different effects.
Best method and filtering: multiple comparisons with the best (MCB)
Frequencies of most common confidence sets using multiple comparisons with the best with a coverage of 0.90, when drivers are known
Confidence sets | Conjunction | No conjunction | |||||||
---|---|---|---|---|---|---|---|---|---|
Diff | PFD | PND | FPF | Diff | PFD | PND | FPF | ||
OT, OT-A | 0.60 | 0.47 | 0.03 | 0.05 | 0.65 | 0.53 | 0.23 | 0.05 | |
DiP, DiP-A, OT, OT-A | 0.02 | 0.08 | - | 0.57 | 0.03 | 0.16 | 0.05 | 0.59 | |
CBN, CBN-A | 0.02 | 0.09 | 0.61 | - | - | 0.04 | 0.32 | 0.01 | |
DiP, DiP-A | - | 0.04 | 0.01 | 0.17 | 0.01 | 0.02 | - | 0.16 | |
CBN, CBN-A, OT, OT-A | 0.02 | 0.02 | 0.07 | - | 0.02 | 0.03 | 0.12 | 0.01 | |
DiP-A, OT, OT-A | - | 0.06 | - | 0.02 | 0.02 | 0.05 | 0.02 | 0.03 | |
OT-A | 0.16 | 0.09 | 0.06 | - | 0.15 | 0.04 | 0.13 | - | |
OT | 0.07 | 0.05 | - | 0.01 | 0.03 | 0.03 | - | - | |
DiP, OT | - | 0.04 | - | 0.08 | 0.01 | 0.05 | - | 0.08 | |
CBN-A | 0.03 | - | 0.06 | - | - | - | 0.01 | - |
Generalized linear modeling of performance measures
The procedures above do not provide a simple and direct way to compare the relative magnitude of the effects of different factors. We can approach this problem using a statistical model for each of the performance measures. I used generalized linear mixed models (GLMM), where data set was a random effect and the rest of the factors in the design were regarded as fixed effects.
Data for the Drivers Known and Drivers Unknown scenarios were modeled separately. All models were fitted using INLA [67,68] (and, to compare results, a subset also with R package MCMCglmm [69] — there were no relevant differences in results). Models were fitted using sum-to-zero contrasts: each main effect parameter is to be interpreted as the (marginal) deviation of that level from the overall mean, and the interaction parameter as the deviation of the linear predictor of the cell mean (for that combination of levels) from the addition of the corresponding main effect parameters. As explained in Additional file 1, we will focus on models with two-way interactions. We will refer to these analyses as the GLMM fits. Further details of the statistical modeling and interpretation of coefficients are provided in Additional file 1.
Results
We first examine the results when Drivers are Known. When the identity of drivers is not known (Drivers Unknown), we need to add the step of filtering or selecting mutations before inferring the restrictions.
Drivers known
There was wide variation in performance: under some Models and with some Methods perfect results were obtained but, for those same Models and S.Sizes, there were choices of Method and Sampling scheme that led to many incorrect decisions, with PFD and PND of 0.7 to 0.9 (see examples of inferred graphs in Additional file 5 and median values for all performance measures for all combinations of factors in Additional file 6): even for the easiest models and largest sample sizes, careful choice of Method can be crucial. There was also a large difference in the number of relations inferred (the transitive closure of the cover relations), whether correct or incorrect: the mean values were 19.1, 18.5, 1.7, 2.6, 7.0 and 8.4 for CBN, CBN-A, DiP, DiP-A, OT, and OT-A, respectively (see Additional file 6).
Table 3 shows the overall ranking of Method and Sampling scheme. OT and OT-A were the best methods according to Diff and PFD. CBN and CBN-A were among the best methods according to PND (in graphs with conjunctions) and DiP and DiP-A according to FPF. This is coherent with the patterns of number of edges (number of dependency relations) inferred: CBN and CBN-A inferred more edges and thus the number of false negatives (FN) decreased, so they had larger sensitivity or recall. But this was done at the cost of increasing the false positives (FP) and, thus, increasing PFD and FPF: a larger fraction of the discoveries were false (precision was smaller) and a larger fraction of the non-existing relationships were regarded as being present (specificity was smaller). DiP and DiP-A showed the opposite trend: these were the methods that inferred the smallest number of relations (in many cases no edges, beyond those from Root, were inferred), leading to a smaller number of false positives (FP), so that a smaller fraction of non-existing relations were regarded as being present, but this was done at the cost of a very large number of FN that affected not only PND but also Diff.
Figure 2, as well as Table 4 and Table 3 (full tables available in Additional file 2) show that OT-A was generally superior to OT, and similar but weaker patterns affected CBN and DiP. The magnitude of the differences between augmented and non-augmented alternatives, however, depended strongly on performance measure and showed interactions with other effects (e.g., S.Time). If we focus on Diff, with OT and DiP there was little to loose from always using the augmented alternative, even if no mutation had a frequency of one, but this was not the case with CBN (see Additional file 2).
The effect of S.Time was marginally very small (Figure 2), and there were also small interactions with Method, which depended on performance measure (see Additional file 2). These differences in performance of each Method with different S.Time in each performance measure can also be observed in the overall ranking of methods (Table 3), as well as in the change in the relative frequencies of OT, OT-A and CBN in confidence sets as a function of S.Time (see Additional file 2). All of these interactions, however, were much smaller than the main effect of Method and should rarely affect which method to choose.
Although the marginal effects of S.Type were generally small (and single cell sampling provided no benefit), overall it seems best to avoid whole tumor sampling with very small detection thresholds. However, as with many other effects, this was reverted with PND. This result is intuitively reasonable: whole tumor sampling at very low thresholds can lead to obtaining samples where we observe together two low frequency events that rarely occur together in the same individual clone (i.e., that do not correspond to a pattern encoded in the true graph), leading to the observation of possible artifacts, but also allowing the detection of co-occurring events of very low frequency. The effects of S.Type were amplified by method (the interaction of Method by S.Type in Figure 2), and this effect was strongly performance measure-dependent (see also Additional file 2). The marginal effect of S.Size was as expected: larger sample size led to better performance with all performance measures. However, the effect in performance was small compared to the effects of choosing a bad method or even the effects of S.Type and S.Time for some of the performance measures (or the effects of non user-controllable factors such as Model). Moreover, the effects of increases in sample size depended on method and performance measure (see also Additional file 2): DiP/DiP-A were, comparatively, the methods that benefited the most from increasing S.Size (except for FPF).
Model also showed interactions with Method, and Model determined the effects of sh (the strength of enforcement of monotonicity, see section “Deviations from monotonicity and genetic context dependence of driver status: sh ” and Table 2) and affected its interaction with Method. If we knew that the true model was McF_6, even when minimizing PND we should then choose OT-A or OT, not CBN or CBN-A (Figure 5(a) and Figure 3). At the same time, for PND the performance of CBN/CBN-A is much less affected by evolutionary model variation.
Finally, Number of Nodes had an effect, which for some performance measures was large, and also presented large interactions with Model, Method and S.Time for some performance measures; these effects will not be examined here further (although a specific one will be addressed in the Discussion), but it is important to notice that the graphs used in this study differed in many ways that are not simply summarized by the number of nodes (see Additional file 1).
Drivers unknown
Many results here were similar to the ones in the previous section, so we will focus on the main differences as well as the added factor of filtering. Overall, performance was worse when Drivers were Unknown.
Ranking of all 144 combinations of method, filtering, and sampling scheme (time, type) when drivers are unknown with respect to each performance measure
Method and sampling | Conjunction | No conjunction | |||||||
---|---|---|---|---|---|---|---|---|---|
Diff | PFD | PND | FPF | Diff | PFD | PND | FPF | ||
S5, OT-A, last, singleC | 1 | 7 | 15 | 60 | 1 | 13 | 13 | 60 | |
S5, OT-A, last, wholeT_0.5 | 2 | 6 | 23 | 59 | 2 | 7 | 23 | 59 | |
J5, OT-A, last, wholeT_0.01 | 3 | 1 | 26 | 71 | 4 | 2 | 29 | 74 | |
S5, OT-A, last, wholeT_0.01 | 4 | 2 | 2 | 94 | 3 | 5 | 4 | 94 | |
S5, OT-A, unif, wholeT_0.01 | 5 | 17 | 33 | 69 | 10 | 17 | 32.5 | 68.5 | |
S5, OT, unif, singleC | 6 | 41 | 61.5 | 28 | 9 | 38 | 55.5 | 20.5 | |
S5, OT-A, unif, singleC | 7 | 50 | 61.5 | 31 | 8 | 49 | 55.5 | 20.5 | |
J1, OT-A, last, singleC | 8 | 15 | 38 | 80 | 6 | 12 | 37 | 75 | |
J1, OT-A, last, wholeT_0.5 | 9 | 10 | 39 | 79 | 7 | 6 | 39 | 71 | |
S5, OT-A, unif, wholeT_0.5 | 10 | 47 | 67 | 29.5 | 14 | 45 | 60.5 | 27.5 | |
S5, OT, unif, wholeT_0.01 | 11 | 11 | 34 | 68 | 11 | 10 | 32.5 | 68.5 | |
J1, OT-A, last, wholeT_0.01 | 13 | 21 | 9 | 109 | 15 | 19 | 8 | 108 | |
S1, OT-A, last, wholeT_0.5 | 18 | 39 | 3 | 117 | 5 | 27 | 3 | 116 | |
S1, OT-A, last, singleC | 21 | 38 | 4 | 120 | 12 | 32 | 2 | 120 | |
S5, OT, last, singleC | 23 | 16 | 51 | 21 | 20 | 8 | 42 | 19 | |
S5, OT, last, wholeT_0.5 | 24 | 3 | 55 | 16 | 23 | 3 | 49 | 17 | |
J1, OT, unif, wholeT_0.01 | 29 | 9 | 36 | 91 | 30 | 14 | 36 | 91 | |
S5, CBN-A, unif, wholeT_0.01 | 30 | 4 | 37 | 89 | 40 | 20 | 40 | 103 | |
J5, OT, last, wholeT_0.01 | 31 | 5 | 59 | 32 | 36 | 4 | 62 | 50 | |
S1, OT-A, last, wholeT_0.01 | 36 | 48 | 1 | 134 | 24 | 50 | 1 | 134 | |
S5, OT, last, wholeT_0.01 | 38 | 8 | 29 | 83 | 34 | 1 | 24 | 82 | |
J1, OT, last, wholeT_0.5 | 47 | 14 | 89 | 61 | 38 | 9 | 83 | 43 | |
J5, OT, last, singleC | 48 | 35 | 101 | 4 | 58 | 43 | 109 | 4.5 | |
J5, DiP-A, unif, singleC | 49 | 125 | 139 | 10.5 | 51 | 128 | 138 | 9 | |
J5, OT, last, wholeT_0.5 | 57 | 20 | 102 | 4 | 75 | 23 | 112 | 4.5 | |
J5, DiP, unif, singleC | 61 | 131 | 143 | 8 | 60 | 132 | 144 | 4.5 | |
J5, DiP, unif, wholeT_0.5 | 64 | 144 | 144 | 4 | 64 | 144 | 143 | 4.5 | |
J5, DiP, last, singleC | 79 | 115 | 140 | 4 | 81 | 119 | 140 | 4.5 | |
J5, DiP, last, wholeT_0.5 | 82 | 137 | 138 | 4 | 83 | 136 | 139 | 4.5 | |
S5, DiP, last, wholeT_0.5 | 83 | 134 | 120 | 9 | 87 | 130 | 110 | 10 | |
J1, DiP, last, wholeT_0.5 | 91 | 132 | 134 | 4 | 89 | 134 | 134 | 4.5 | |
J1, DiP, last, singleC | 92 | 110 | 135 | 4 | 93 | 115 | 132 | 4.5 | |
S1, OT-A, unif, wholeT_0.01 | 102 | 65 | 5 | 136 | 70 | 59 | 5 | 136 | |
S1, OT, unif, wholeT_0.01 | 109 | 62 | 6 | 135 | 74 | 55 | 6 | 135 | |
S1, CBN-A, unif, wholeT_0.01 | 137 | 82 | 7 | 143 | 140 | 81 | 7 | 143 | |
S1, CBN-A, last, wholeT_0.01 | 142 | 93 | 10 | 137 | 139 | 94 | 12 | 137 | |
S1, CBN, unif, wholeT_0.01 | 143 | 84 | 8 | 144 | 143 | 90 | 9 | 144 |
Frequencies of most common confidence sets using multiple comparisons with the best with a coverage of 0.90, when drivers are unknown
Confidence sets | Conjunction | No conjunction | |||||||
---|---|---|---|---|---|---|---|---|---|
Diff | PFD | PND | FPF | Diff | PFD | PND | FPF | ||
S1:OT-A, S5:OT-A | 0.01 | - | 0.04 | - | 0.01 | - | 0.07 | - | |
S1:OT, S1:OT-A | 0.07 | 0.07 | 0.11 | - | 0.08 | 0.07 | 0.22 | - | |
S1:OT-A | 0.02 | - | 0.06 | - | 0.02 | - | 0.05 | - | |
S5:OT-A | 0.05 | - | - | - | 0.06 | - | - | - | |
S5:OT, S5:OT-A | 0.10 | 0.03 | - | - | 0.12 | 0.02 | - | - | |
S5:DiP, S5:DiP-A, S5:OT, S5:OT-A | 0.04 | - | - | - | 0.05 | - | - | - | |
S1:OT, S1:OT-A, S5:OT, S5:OT-A | 0.02 | - | 0.02 | - | 0.05 | - | 0.07 | - | |
J5:OT, J5:OT-A, S5:OT, S5:OT-A | 0.01 | 0.06 | - | - | - | 0.04 | - | - | |
S1:CBN, S1:CBN-A | 0.01 | - | 0.05 | - | 0.01 | 0.01 | 0.01 | - | |
S1:CBN, S1:CBN-A, S1:OT, S1:OT-A | - | 0.01 | 0.25 | - | 0.02 | - | 0.14 | - | |
J1:DiP, J1:OT, J5:DiP, J5:OT, S1:DiP, S5:DiP, S5:OT | - | 0.02 | - | 0.04 | - | 0.03 | - | 0.05 |
As in the Drivers Known scenario, DiP and DiP-A led to the smallest FPF (again, FPF is of minor relevance compared to PND and PFD), but did so at the cost of the other performance measures. However, when Drivers were Unknown DiP and DiP-A were better performers than CBN and CBN-A for Diff.
There were interactions of Filter by Method (see Figure 9; also figures in Additional file 2), and their direction and magnitude depended on performance measure, but these interactions were not large enough to revert our preferences of methods. S.Type followed similar patterns as seen for the Drivers Known scenario, although the effects were of larger magnitude, and its interactions with Method or Filter (Figure 9; see also Additional file 2), even when present, produced no change in the ordering of preferences of Filtering, Method, or S.Type. The effects of S.Size or its interactions with Method were similar to those in the Drivers Known scenario, although S.Size seemed to be more important for decreasing PFD when Drivers are Unknown (Figure 9 and figures in Additional file 2).
Other effects and interactions were different with respect to the Drivers Known scenario, illustrating that the need to filter passengers can lead to counterintuitive and much more complicated interpretation of results. There was an interaction of Filter by Model (Figure 10): the best choice of S5 vs. J5 for performance measures Diff and PFD depended on the Model, as shown in Figure 11. Moreover, we can see that the very poor results in measure FPF of model McF_6 were largely due to its terrible results with Filter S1 (see also below). There was also an interaction Model by Method, and OT/OT-A were now superior in all performance measures (except FPF) not only overall, but also virtually for every one of the four models, as we can see from panel b) in Figure 5.
The effect of S.Time for all performance measures except Diff was reverted and had larger magnitude compared to the Drivers Known scenario. The cause is that S.Time affected the number of genes that were selected and thus the number of false positives (FP) and false negatives (FN). As seen in panel (a) of Figure 8, under S.Time = uniform fewer genes were selected for all filtering methods (the number of genes selected is independent of Method but depends on Filtering). Under uniform sampling many of the samples have few mutations (those corresponding to the early stages of the disease), and thus fewer genes are above the filtering thresholds. As fewer genes are selected under S.Time = uniform, the number of FP goes down, and thus both PFD and FPF showed a decrease (and the opposite pattern was observed for PND). None of these phenomena were present when Drivers are Known, because there was no filtering step there. For Diff, however, the marginal effect of S.Time=uniform was still slightly better than that of last-period sampling, as was when Drivers are Known. These differences of the effect of S.Time between the Drivers Known and Drivers Unknown scenarios explain also the differences in the patterns of interactions of Method by S.Time and Model by S.Time (Figures 9 and 10; also figures in Additional file 2).
Another difference with the Drivers Known scenario are some of the interactions between Model and sh in several performance measures, probably related to the filtering step. Specially with PFD and FPF, sh = 0 led to a large decrease in error rates as seen in Figures 7(b) and 10. Panel (b) of Figure 8 shows that, for all filtering procedures, sh=0 led to smaller number of genes being selected with models Bozic and Exp (but had no effect on the McF models), and Bozic and exp were the two models where sh=0 led to a larger relative improvement in performance in PFD and FPF (Figure 10 and panel b of Figure 7). Interestingly, the difference between the two levels of sh is more pronounced in S1, and suggests that sh=Inf allows the accumulation of larger numbers of mutated genes that pass the less stringent filters, but that have too little signal for the inference of the restrictions (in contrast, sh=0 leads to the accumulation of mutations with too low a frequency to even pass the filters). That sh should have virtually no effect on number of genes in the McF models has been explained before. These interactions, therefore, highlight the complex and counterintuitive relationships between Model, sh, and Filtering, that then cascade onto the overall performance differences between Methods or Sampling schemes.
Discussion
Most of the above results cannot be compared to any other studies, since these factors have not been considered before. It has been previously discussed [15,32] that very good reconstructions of oncogenetic trees are achievable with realistic sample sizes, and we have seen that, at least under several scenarios, it is possible to obtain perfect reconstructions even with sample sizes as small as 100.
The results concerning the superiority of OT with respect to CBN contrast with those of Hainke et al. [13], who find that CBN outperforms OT. A more detailed discussion is provided in Additional file 3, but these differences are attributable to [13] basing their conclusion on a single graph with a small number of nodes for each scenario. Additional file 3 shows the within-data set difference in the Diff measure between CBN and OT separately for the different combinations of Model, sh (see section “Deviations from monotonicity and genetic context dependence of driver status: sh ” and Table 2), True Graph, and Sampling scheme (under the Drivers Known scenario only, since [13] do not consider passengers): OT systematically outperformed CBN except for Graph 7A when sampling last (under all models except McF_6), and Graph 7B when sampling last under the Bozic and exponential models. This pattern was reproduced when I fitted the OT models with the Rtreemix package [54] (see Additional file 3). The single graph used by [13] for the non-conjunction case (their Figure A1) contains five nodes, and the single graph they used for the case with conjunctions contains four nodes (their Figure A7), a very small number of nodes compared to the graphs that are seen on the literature (see Methods section). Our graphs 7A and 7B, those where CBN outperformed OT in certain scenarios, are the closest in number of nodes to the graphs in [13].
Given that the results of Hainke et al. [13] are, thus, probably not really contradictory with ours, is the recommendation that practitioners generally use OT instead of CBN still valid? Yes, since most graphs in the literature (including studies involving citogenetic bands and genes) contain many more than four or five nodes, and we cannot be sure if the evolutionary model is one that would favor using CBN. In addition, we have focused only on the comparison between CBN and OT, since those are the only variants used in [13]: if we included OT-A in the comparison, CBN (or CBN-A) would then very rarely be better alternatives (see, e.g., Additional file 4 or Additional file 3). However, this apparent difference in results emphasizes the need for considering at least a few different scenarios with regards to potentially key variables, and suggests that a through examination of the impact of graph characteristics (and its interaction with evolutionary model and sampling scheme) on method performance is warranted.
This paper is also one of the first to explicitly connect evolutionary models with restrictions on the order of mutations. Recently, in [39] a simulation tool has been described where restrictions are incorporated into the evolutionary model of [41]; our approach is more general as [39] are limited to four drivers and no passengers whereas, in addition to passengers and other evolutionary models, we can specify restrictions in the order of mutations using arbitrary graphs and allowing for a range of deviations from monotonicity. In fact, one of the attractive features of OTs, CBNs, and Progression Networks is their mechanistic interpretation as graphs that encode restrictions in the order in which driver mutations can accumulate [12,13,15,32]. And one major result of this paper is that inferring those restrictions can be strongly affected by evolutionary model (including deviations from monotonicity) and sampling scheme, and that the relative effect of these factors depends on the performance measure used. Yet restrictions in the order of driver mutations and evolutionary models are virtually always examined separately.
There is a rich literature about tumor progression models that focuses on the consequences of drivers, passengers, and variation in selection pressures [20,41,42,44,70,71], and a largely separate body of work [13-15,19,21,31,32,45] that deals with understanding the restrictions and order of accumulation of mutations (but see [6] for a connection between the λ _{ i } of CBNs and selection coefficients, in the context of the Fisher-Wright model of tumor progression in [49]). The work of Cheng et al. [7,8] tries to infer the order of mutations within a explicit evolutionary model of tumor progression; unfortunately, no software is available, and thus comparisons are not possible.
Focusing on methods with available software, the actual values (and, thus, interpretation) of the conditional probabilities inferred by OTs or the λ _{ i } parameter for the waiting time to event i, in CBNs, will be a complex interplay between the restrictions encoded in the graphs and the details of the tumor progression model as well as the sampling scheme used. For both OTs and CBNs we should expect estimated λs and conditional probabilities to vary by node level or depth (where level or depth refers here to how many edges there are in the path to the root): deeper nodes will show smaller values and, for a given depth, λs/conditional probabilities should be larger for those nodes than “unlock” more downwards mutations. The strength of this effect will increase with the number of nodes along the largest path along the graph, especially when the evolutionary model and sampling scheme result in strong selection for clones with many drivers, and we should see competition between multiple nodes that descend from the same parent.
As an example of the impact of evolutionary model and sampling scheme on the observable consequences of the restrictions encoded in the graphs, Additional file 7 shows inferred oncogenetic trees with their estimated conditional probabilities for each of the three true graphs without conjunctions (trees), under several scenarios, including inferred oncogenetic trees with a sample size of 26000 (to minimize sampling variation). The inferred graphs are perfect, or almost perfect, reconstructions of the topology of the true graph, but the estimates of the conditional probabilities show large differences. The variations are in directions we would expect both between S.Time and among models, as well as among nodes. Even if the above results are intuitively reasonable, they highlight that whereas the topologies of the graphs (the partial orders) encode constraints in the order of mutations, the conditional probabilities (or λs) we estimate and, most importantly, the patterns of co-occurrence of mutations and the sets of clones we observe, will depend crucially on the evolutionary model and sampling scheme. Since the topology reconstructions depend on the patterns observed our inferences will be strongly impacted by evolutionary model (and sampling scheme), as we have seen repeatedly in the results. Moreover, examining the consequences of sampling scheme (S.Time and S.Type) and the detrimental effects of having to separate drivers from passengers on the quality of our inferences can only be meaningfully considered with respect to an evolutionary tumor progression model that generates the data.
As mentioned in the Introduction, the above interaction between order restrictions and evolutionary model, and the unavoidable need to interpret parameters in the context of a given evolutionary model, are coherent with the limitations pointed out by Sprouffske et al. [37]: that oncogenetic tree models (and related models) are not really evolutionary models and do not represent ancestral relationships, but only summarize patterns of co-occurrences of mutations across samples. Virtually all studies of methods for inferring order restrictions are susceptible to this criticism, since they simulate data directly from the generative (but non-evolutionary) OT/CBN/Progression Network model. However, the design I have used here completely overcomes this limitation: I have simulated the data using plausible evolutionary models that incorporate the restrictions in the order of mutations via a straightforward effect on the fitness of clones. Moreover, deviations from monotonicity are not added to the model just as an unexplained error term, but are an integral part of the evolutionary model that can be related, for instance, to the genetic context-dependence of the driver/passenger status.
Sprouffske et al. [37] conclude that cross-sectional data can be misleading if we try to infer the order of mutations. But this conclusion is based on a design where a single OT is inferred from a cross-sectional sample where mutations are not restricted to obey a pre-specified set of restrictions. Thus, it is not surprising that the OT fit does not do well. The results of [37] of course highlight that if different subjects have different sets of order restrictions, then no single OT will capture these patterns, a limitation that is already recognized in the early literature on ongenetic trees [14,15], and that has prompted the development of mixtures of oncogenetic trees [24,25,72]. But, by themselves, their results do not show that OTs (or CBNs or DiPs) from cross-sectional data cannot fare well if there is a true underlying set of restrictions that can be represented as a single graph. Quite to the contrary, I have shown here, embedding the restrictions in evolutionary models, that they can do very well and even recover the exact underlying graph (at least under certain scenarios). Moreover, [37] do not show that any particular within-subject method is actually capable of recovering the true paths from their data (they sidestep that problem altogether).
The two key remaining questions to be answered regarding the usage of cross-sectional data, then, are two: 1) whether the accumulation of mutations in cancer progression can be reasonably represented by a single graph that encapsulates restrictions; 2) if 1) does not hold, whether cross-sectional methods such as mixtures of oncogenetic trees can recover the set of different restrictions. If the answer to 1) is positive, the results of this paper indicate that we have methods that can recover those relationships, and these results also highlight possible avenues to improve them. But question 1) is one that neither this study nor the one of [37] can answer (I simulated data assuming a Yes to that question and [37] assuming a No). Question 2) remains to be thoroughly addressed, and neither this paper nor [37] shed light on the matter. If the answer is negative, then we need to start focusing on within-individual data, which are much harder to obtain. Nevertheless, the approaches in this paper provide a principled and general way to address that question by simulating data under scenarios were there is no single set of restrictions in common to all subjects, and examining the consequences both for our methods of inferring trajectories and for the data patterns themselves (so as to try to infer, from them, whether or not there is a single set of restrictions).
Conclusions
- 1.
Method and sampling choice should be guided by the performance measure considered most important: no combination of Method, Filtering, and Sampling scheme excels in all performance measures. This is not unexpected, but it is worth emphasizing that each performance measure is sensitive to a different kind of deviation, and the results in this paper show that characterizing behavior with only one or two performance measures could have been deceptive. Moreover, as we have seen, the relative strengths of each method are better captured and understood using different performance measures.
- 2.
In terms of method choice, a very simple summary is (see also examples of inferred graphs in Additional file 5): CBN tends to return graphs with too many edges, including too many conjunctions, DiP tends to return graphs with too few edges between non-root nodes, and OT does a good overall job even though it will fail, by construction, to return any conjunction. In more detail, OT and OT-A are the best methods, except if we are particularly interested in minimizing PND and we suspect conjunctions are present (when we might want to consider CBN) or FPF (when we might want to consider DiProg —recalling that FPF is generally of minor value compared to measures PND and PFD). Since it is impossible for OT to return any conjunctions, further research on computationally efficient methods to recover conjunctions is sorely needed.
- 3.
Using frequency-based statistics when we do not know which mutations are true passengers can lead to a heavy performance penalty mainly in the form of failure to discover existing restrictions. In addition, having to filter genes makes it much harder to intuitively understand, and reason about, what is likely to happen in any scenario, and this in turn makes interpretation of results and reconciliation of output from different methods much harder. Thus, it probably pays off to try to use other approaches that incorporate information about non-silent mutation rates, pathway information together with combinatorial properties of drivers in pathwayws, or functional consequences of mutations to differentiate drivers from passengers [73-78]. It might not always be possible to use these other methods. If we need to rely on frequency-based approaches, selecting those mutations with a frequency larger than 5% is an overall reasonable choice (but not the single best choice for any performance measure other than Diff).
- 4.
Sampling time and type, by themselves, had minor effects compared to, say, filtering or method choice (and we will rarely have control over these factors when we use data already available in databases). However, we might have information about characteristics of the tumor that indicate that it is in an exponential (the Bozic or exp models) or logistic-like growth (the McF models) phase; sampling as late as possible is to be preferred for the first cases, whereas trying to obtain samples distributed over different stages of the disease is best in the second. The best choice of sampling time, however, will depend on performance measure and whether or not we are certain about which are the driver genes.
- 5.
Single cell sampling is about as good as whole tumor sampling, unless we use whole tumor sampling with extremely small detection thresholds (which leads to poorer performance, except for performance measure PND).
- 6.
Although with S.Size the larger the better, its effect is relatively minor for OT and CBN (not so for DiP), a result that agrees with those in [13], specially for some performance measures. In particular, resources might be better spent trying to be certain about the identity of the true drivers than increasing sample size from 100 to 1000.
- 7.
Data augmentation is not always the right choice, although with OT and DiP there is little to loose from always using data augmentation but potentially a lot to loose from not using augmentation (see Additional file 2). Unfortunately, simple rules of thumb like “always use data augmentation when at least one gene has a frequency of one, and never otherwise” do not work well, especially across all methods, and appropriate choice warrants further study.
- 8.
Evolutionary model had a strong impact on method performance, both directly and indirectly through its interactions with other factors (such as sampling time or filtering). Assessment of method performance via simulations in this and related problems should thus be done incorporating order restrictions within plausible models of tumor evolution, which also allows us to naturally examine the effects of other factors, such as sampling. In particular, a framework similar to the one used here could be applied to scenarios where order restrictions differ between subgroups of individuals.
Availability of source code and data
Declarations
Acknowledgements
C. Lazaro-Perea, L. del Peso and J. Poyatos for comments on the ms. C. Lazaro-Perea for help with diagrams and tables. I. B. Diaz for rooting many graphs. R. Salomón for discussion about CBNs and OTs. H. S. Farahani, T. Graham and R. S. Datta, M. Gerstung, K. Hainke, C. D. McFarland, W. Mather and L. Tsimring, T. Sakoparnig, and K. Sprouffske for answering questions about their papers and methods and/or providing code of the implementations of their algorithms. Two anonymous reviewers for their comments that have substantially improved the ms. Work partially supported by Project BIO2009-12458 from the Spanish MINECO; computer hardware partially supported by Project BIO2009-12458 from the Spanish MINECO. The funding body had no role in study design, data collection, analysis and interpretation, preparation of the manuscript, or decision to submit the manuscript for publication.
Authors’ Affiliations
References
- Merlo LMF, Pepper JW, Reid BJ, Maley CC. Cancer as an evolutionary and ecological process. Nat Rev Cancer. 2006; 6(12):924–35. doi:10.1038/nrc2013.View ArticlePubMedGoogle Scholar
- Hanahan D, Weinberg R. Hallmarks of Cancer: The Next Generation. Cell. 2011; 144(5):646–74. doi:10.1016/j.cell.2011.02.013.View ArticlePubMedGoogle Scholar
- Weinstein IB, Joe A. Oncogene addiction. Cancer Res. 2008; 68(9):3077–803080. doi:10.1158/0008-5472.CAN-07-3293.View ArticlePubMedGoogle Scholar
- Fearon E, Vogelstein B. A genetic model for colorectal tumorigenesis. Cell. 1990; 61:759–67.View ArticlePubMedGoogle Scholar
- Raphael B, Vandin F. Simultaneous Inference of Cancer Pathways and Tumor Progression from Cross-Sectional Mutation Data In: Sharan R, editor. Heidelberg, Germany: Springer: 2014. p 250–64. http://link.springer.com/chapter/10.1007/978-3-319-05269-4_20.
- Gerstung M, Eriksson N, Lin J, Vogelstein B, Beerenwinkel N. The Temporal Order of Genetic and Pathway Alterations in Tumorigenesis. PLoS ONE. 2011; 6(11):27136. doi:10.1371/journal.pone.0027136.View ArticleGoogle Scholar
- Cheng Y-K, Beroukhim R, Levine RL, Mellinghoff IK, Holland EC, Michor F. A mathematical methodology for determining the temporal order of pathway alterations arising during gliomagenesis. PLoS Comput Biol. 2012; 8(1):1002337. doi:10.1371/journal.pcbi.1002337.View ArticleGoogle Scholar
- Attolini C, Cheng Y, Beroukhim R, Getz G, Abdel-Wahab O, Levine RL, et al.A mathematical framework to determine the temporal sequence of somatic genetic events in cancer. Proc Nat Acad Sci. 2010; 107(41):17604–9. doi:10.1073/pnas.1009117107/-/DCSupplemental.www.pnas.org/cgi/doi/10.1073/pnas.1009117107.View ArticlePubMedPubMed CentralGoogle Scholar
- Longerich T, Mueller MM, Breuhahn K, Schirmacher P, Benner A, Heiss C. Oncogenetic tree modeling of human hepatocarcinogenesis. Int J Cancer J Int du Cancer. 2012; 130(3):575–83. doi:10.1002/ijc.26063.View ArticleGoogle Scholar
- Roulston A, Muller WJ, Shore GC. BIM, PUMA, and the achilles’ heel of oncogene addiction. Sci signaling. 2013; 6(268):12. doi:10.1126/scisignal.2004113.View ArticleGoogle Scholar
- Luo J, Solimini NL, Elledge SJ. Principles of cancer therapy: oncogene and non-oncogene addiction. Cell. 2009; 136(5):823–37. doi:10.1016/j.cell.2009.02.024.View ArticlePubMedPubMed CentralGoogle Scholar
- Farahani HS, Lagergren J. Learning oncogenetic networks by reducing to mixed integer linear programming. PloS one. 2013; 8(6):65773. doi:10.1371/journal.pone.0065773.View ArticleGoogle Scholar
- Hainke K, Rahnenführer J, Fried R. Cumulative disease progression models for cross-sectional data: A review and comparison. Biometrical J Biometrische Zeitschrift. 2012; 54(5):617–40. doi:10.1002/bimj.201100186.View ArticleGoogle Scholar
- Radmacher MD, Simon R, Desper R, Taetle R, Schaffer AA, Nelson MA. Graph Models of Oncogenesis with an Application to Melanoma. J Theor Biol. 2001; 212(4):535–48. doi:10.1006/jtbi.2001.2395.View ArticlePubMedGoogle Scholar
- Desper R, Jiang F, Kallioniemi OP, Moch H, Papadimitriou CH, Schäffer AA. Inferring tree models for oncogenesis from comparative genome hybridization data. J Comput Biol. 1999; 6(1):37–51.View ArticlePubMedGoogle Scholar
- Vogelstein B, Fearon ER, Hamilton SR, Kern SE, Preisinger AC, Leppert M, et al. Genetic Alterations during Colorectal-Tumor Development. New England J Med. 1988; 319(9):525–32. doi:10.1056/NEJM198809013190901.View ArticleGoogle Scholar
- Desper R, Jiang F, Kallioniemi OP, Moch H, Papadimitriou CH, Schäffer AA. Distance-based reconstruction of tree models for oncogenesis. J Comput Biol. 2000; 7(6):789–803. doi:10.1089/10665270050514936.View ArticlePubMedGoogle Scholar
- von Heydebreck A, Gunawan B, Füzesi L. Maximum likelihood estimation of oncogenetic tree models. Biostatistics. 2004; 5(4):545–56.View ArticlePubMedGoogle Scholar
- Youn A, Simon R. Estimating the order of mutations during tumorigenesis from tumor genome sequencing data. Bioinformatics (Oxford, England). 2012; 28(12):1555–61. doi:10.1093/bioinformatics/bts168.View ArticlePubMed CentralGoogle Scholar
- Beerenwinkel N, Eriksson N, Sturmfels B. Conjunctive Bayesian networks. Bernoulli. 2007; 13(4):893–909. doi:10.3150/07-BEJ6133.View ArticleGoogle Scholar
- Simon R, Desper R, Papadimitriou CH, Peng A, Alberts DS, Taetle R, et al. Chromosome abnormalities in ovarian adenocarcinoma: III. Using breakpoint data to infer and test mathematical models for oncogenesis. Genes, Chromosomes and Cancer. 2000; 28(1):106–20. doi:10.1002/(SICI)1098-2264(200005)28:1%3C106::AID-GCC13%3E3.0.CO;2-S.View ArticlePubMedGoogle Scholar
- Hjelm M, Höglund M, Lagergren J. New probabilistic network models and algorithms for oncogenesis. J Comput Biol. 2006; 13(4):853–65. doi:10.1089/cmb.2006.13.853.View ArticlePubMedGoogle Scholar
- Bulashevska S, Szakacs O, Brors B, Eils R, Kovacs G. Pathways of urothelial cancer progression suggested by Bayesian network analysis of allelotyping data. Int J Cancer. 2004; 110(6):850–6. doi:10.1002/ijc.20180.View ArticlePubMedGoogle Scholar
- Beerenwinkel N, Rahnenführer J, Däumer M, Hoffmann D, Kaiser R, Selbig J, et al. Learning multiple evolutionary pathways from cross-sectional data. J Comput Biol: J Comput Mol Cell Biol. 2005; 12(6):584–98. doi:10.1089/cmb.2005.12.584.View ArticleGoogle Scholar
- Tofigh A, Sjolund E, Hoglund M, Lagergren J. A global structural EM algorithm for a model of cancer progression In: Shawe-Taylor J, Zemel RS, Bartlett PL, Pereira F, Weinberger KQ, editors. Advances in Neural Information Processing Systems. Red Hook, NY, USA: Curran Associates, Inc: 2011. p. 163–71. https://papers.nips.cc/paper/4397-a-global-structural-em-algorithm-for-a-model-of-cancer-progression.pdf.Google Scholar
- Vogelstein B, Papadopoulos N, Velculescu VE, Zhou S, Diaz LA, Kinzler KW. Cancer Genome Landscapes. Science. 2013; 339(6127):1546–58. doi:10.1126/science.1235122.View ArticlePubMedPubMed CentralGoogle Scholar
- Tomasetti C, Vogelstein B, Parmigiani G. Half or more of the somatic mutations in cancers of self-renewing tissues originate prior to tumor initiation. Proc Nat Acad Sci USA. 2013; 110(6):1999–2004. doi:10.1073/pnas.1221068110.View ArticlePubMedPubMed CentralGoogle Scholar
- Wood LD, Parsons DW, Jones S, Lin J, Sjoblom T, Leary RJ, et al. The Genomic Landscapes of Human Breast and Colorectal Cancers. Science. 2007; 318(5853):1108–13. doi:10.1126/science.1145720.View ArticlePubMedGoogle Scholar
- Greenman C, Stephens P, Smith R, Dalgliesh GL, Hunter C, Bignell G, et al.Patterns of somatic mutation in human cancer genomes. Nature. 2007; 446(7132):153–8. doi:10.1038/nature05610.View ArticlePubMedPubMed CentralGoogle Scholar
- Greaves M, Maley CC. Clonal evolution in cancer. Nature. 2012; 481(7381):306–13. doi:10.1038/nature10762.View ArticlePubMedPubMed CentralGoogle Scholar
- Gerstung M, Baudis M, Moch H, Beerenwinkel N. Quantifying cancer progression with conjunctive Bayesian networks. Bioinf (Oxford, England). 2009; 25(21):2809–15. doi:10.1093/bioinformatics/btp505.View ArticleGoogle Scholar
- Szabo A, Boucher KM. Oncogenetic trees In: Tan W-Y, Hanin L, editors. Handbook of Cancer Models with Applications. Singapore: World Scientific Publishing Co. Pte. Ltd.: 2008. p. 1–24. Chap. 1. http://www.worldscibooks.com/lifesci/6677.html.Google Scholar
- Sakoparnig T, Beerenwinkel N. 2012. Efficient sampling for Bayesian inference of conjunctive Bayesian networks, Vol. 28. doi:10.1093/bioinformatics/bts433.Google Scholar
- Gunawan B, von Heydebreck A, Sander B, Schulten HJ, Haller F, Langer C, et al. An oncogenetic tree model in gastrointestinal stromal tumours (GISTs) identifies different pathways of cytogenetic evolution with prognostic implications. J Pathol. 2007; 211(4):463–70. doi:10.1002/path.2128.View ArticlePubMedGoogle Scholar
- Diep CB, Kleivi K, Ribeiro FR, Teixeira MR, Lindgjaerde OC, Lothe RA. The order of genetic events associated with colorectal cancer progression inferred from meta-analysis of copy number changes. Genes Chromosomes Cancer. 2006; 45(1):31–41. doi:10.1002/gcc.20261.View ArticlePubMedGoogle Scholar
- Rieker RJ, Penzel R, Aulmann S, Blaeker H, Morresi-Hauf A, Hecker E, et al.Oncogenetic tree models based on cytogenetic data: new insights into the development of epithelial tumors of the thymus. Cancer Genet Cytogenet. 2005; 158(1):75–80. doi:10.1016/j.cancergencyto.2004.08.026.View ArticlePubMedGoogle Scholar
- Sprouffske K, Pepper JW, Maley CC. Accurate reconstruction of the temporal order of mutations in neoplastic progression. Cancer Prev Res (Philadelphia, Pa.) 2011; 4(7):1135–44. doi:10.1158/1940-6207.CAPR-10-0374.View ArticleGoogle Scholar
- Yap TA, Gerlinger M, Futreal PA, Pusztai L, Swanton C. Intratumor heterogeneity: seeing the wood for the trees. Sci Translational Med. 2012; 4(127):127–10. doi:10.1126/scitranslmed.3003854.View ArticleGoogle Scholar
- Reiter J, Bozic I, Chatterjee K, Nowak M. TTP: tool for tumor progression In: Sharygina N, Veith H, editors. Computer Aided Verification, Lecture Notes in Computer Science. Berlin, Heidelberg, Germany: Springer Berlin Heidelberg: 2013. p. 101–6. http://link.springer.com/chapter/10.1007/978-3-642-39799-8_6 http://dx.doi.org/10.1007/978-3-642-39799-8_6.Google Scholar
- Eberwine J, Sul J-Y, Bartfai T, Kim J. The promise of single-cell sequencing. Nat Methods. 2013; 11(1):25–7. doi:10.1038/nmeth.2769.View ArticleGoogle Scholar
- Bozic I, Antal T, Ohtsuki H, Carter H, Kim D, Chen S, et al. Accumulation of driver and passenger mutations during tumor progression. Proc Nat Acad Sci USA. 2010; 107:18545–50. doi:10.1073/pnas.1010978107.View ArticlePubMedPubMed CentralGoogle Scholar
- McFarland CD, Korolev KS, Kryukov GV, Sunyaev SR, Mirny LA. Impact of deleterious passenger mutations on cancer progression. Proc Nat Acad Sci USA. 2013; 110(8):2910–5. doi:10.1073/pnas.1213968110.View ArticlePubMedPubMed CentralGoogle Scholar
- Mather WH, Hasty J, Tsimring LS. Fast stochastic algorithm for simulating evolutionary population dynamics. Bioinf (Oxford, England). 2012; 28(9):1230–8. doi:10.1093/bioinformatics/bts130.View ArticleGoogle Scholar
- Datta RS, Gutteridge A, Swanton C, Maley CC, Graham TA. Modelling the evolution of genetic instability during tumour progression. Evolutionary App. 2013; 6(1):20–33. doi:10.1111/eva.12024.View ArticleGoogle Scholar
- Szabo A, Boucher K. Estimating an oncogenetic tree when false negatives and positives are present. Math Biosci. 2002; 176(2):219–36. doi:10.1016/S0025-5564(02)00086-X.View ArticlePubMedGoogle Scholar
- Sweeney C, Boucher KM, Samowitz WS, Wolff RK, Albertsen H, Curtin K, et al. Oncogenetic tree model of somatic mutations and DNA methylation in colon tumors. Genes Chromosomes Cancer. 2009; 48(1):1–9. doi:10.1002/gcc.20614.View ArticlePubMedPubMed CentralGoogle Scholar
- Jiang HY, Huang ZX, Zhang XF, Desper R, Zhao T. Construction and analysis of tree models for chromosomal classification of diffuse large B-cell lymphomas. World J Gastroenterol. 2007; 13(11):1737–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Pathare S, Schäffer AA, Beerenwinkel N, Mahimkar M. Construction of oncogenetic tree models reveals multiple pathways of oral cancer progression. Int J Cancer J Int du Cancer 2003. 2009; 124(12):2864–71. doi:10.1002/ijc.24267.View ArticleGoogle Scholar
- Beerenwinkel N, Antal T, Dingli D, Traulsen A, Kinzler KW, Velculescu VE, et al. Genetic progression and the waiting time to cancer. PLoS Comput Biol. 2007; 3(11):225. doi:10.1371/journal.pcbi.0030225.View ArticleGoogle Scholar
- Sjoblom T, Jones S, Wood LD, Parsons DW, Lin J, Barber TD, et al. The Consensus Coding Sequences of Human Breast and Colorectal Cancers. Science. 2006; 314(5797):268–74. doi:10.1126/science.1133427.View ArticlePubMedGoogle Scholar
- Vermeulen L, Morrissey E, van der Heijden M, Nicholson AM, Sottoriva A, Buczacki S, et al. Defining stem cell dynamics in models of intestinal tumor initiation. Sci (New York, N.Y.) 2013; 342(6161):995–8. doi:10.1126/science.1243148.View ArticleGoogle Scholar
- Korolev KS, Xavier JB, Gore J. Turning ecology and evolution against cancer. Nat Rev Cancer. 2014; 14(5):371–80. doi:10.1038/nrc3712.View ArticlePubMedGoogle Scholar
- Szabo A, Pappas L. Oncotree: Estimating oncogenetic trees. 2013. http://cran.r-project.org/package=Oncotree.
- Bogojeska J. Rtreemix: Mutagenetic tree mixture models. 2014. http://master.bioconductor.org/packages/release/bioc/html/Rtreemix.html.
- Bogojeska J, Alexa A, Altmann A, Lengauer T, Rahnenfuhrer J. Rtreemix: an R package for estimating evolutionary pathways and genetic progression scores. Bioinformatics. 2008; 24(20):2391–2. doi:10.1093/bioinformatics/btn410.View ArticlePubMedPubMed CentralGoogle Scholar
- Dudoit S, Gentleman RC, Quackenbush J. Open source software for the analysis of microarray data. BioTechniques. 2003; Suppl:45–51.PubMedGoogle Scholar
- Vincent A, Charette S. Freedom in bioinformatics. Frontiers in genetics. 2014; 75(23):7537–41. doi:10.1128/AEM.01541-09.Google Scholar
- Barnes N. Publish your computer code: it is good enough. Nature. 2010; 467(7317):753. doi:10.1038/467753a.View ArticlePubMedGoogle Scholar
- Peng RD. Reproducible Research in Computational Science. Science. 2011; 334:1226–7. doi:10.1126/science.1213847.View ArticlePubMedPubMed CentralGoogle Scholar
- Smith R, Ventura D, Prince JT. Novel algorithms and the benefits of comparative validation. Bioinf (Oxford, England). 2013; 29(12):1583–5. doi:10.1093/bioinformatics/btt176.View ArticleGoogle Scholar
- Davis J, Goadrich M. The relationship between Precision-Recall and ROC curves. Proc 23rd Int Conference Machine Learning - ICML’06. New York, NY, USA: ACM; 2006, pp. 233–40. doi:10.1145/1143844.1143874.Google Scholar
- Pepe M. The Statistical Evaluation of Medical Tests for Classification and Prediction. Oxford, UK: Oxford University Press; 2003.Google Scholar
- Gentle JE. Matrix Algebra. New York: Springer; 2007.View ArticleGoogle Scholar
- Hsu JC. Multiple Comparisons: Theory and Methods: Chapman and Hall/CRC Press; 1996.Google Scholar
- Hsu J. Simultaneous inference with respect to the best treatment in block designs. J Am Stat Assoc. 1982; 77(378):461–7.View ArticleGoogle Scholar
- Hsu J, Qiu P, Hin L, Mutti D, Zadnik K. Multiple comparisons with the best ROC curve. Lecture Notes-Monograph Series. 2004; 47:65–75. http://projecteuclid.org/euclid.lnms/1196285619.View ArticleGoogle Scholar
- Rue H. v., Martino S, Chopin N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J R Stat Soc: Ser B (Stat Methodology). 2009; 71(2):319–92. doi:10.1111/j.1467-9868.2008.00700.x.View ArticleGoogle Scholar
- Fong Y, Rue H. v., Wakefield J. Bayesian inference for generalized linear mixed models. Biostatistics (Oxford, England). 2010; 11(3):397–412. doi:10.1093/biostatistics/kxp053.View ArticlePubMed CentralGoogle Scholar
- Hadfield J. MCMC Methods for Multi-response Generalized Linear Mixed Models: The MCMCglmm R Package. J Stat Software. 2010; 33(2):1–22.View ArticleGoogle Scholar
- Durrett R, Foo J, Leder K, Mayberry J, Michor F. Evolutionary dynamics of tumor progression with random fitness values. Theor Population Biol. 2010; 78(1):54–66. doi:10.1016/j.tpb.2010.05.001.View ArticleGoogle Scholar
- Durrett R, Foo J, Leder K, Mayberry J, Michor F. Intratumor heterogeneity in evolutionary models of tumor progression. Genetics. 2011; 188(2):461–77. doi:10.1534/genetics.110.125724.View ArticlePubMedPubMed CentralGoogle Scholar
- Bogojeska J, Lengauer T, Rahnenführer J. Stability analysis of mixtures of mutagenetic trees. BMC Bioinf. 2008; 9:165. doi:10.1186/1471-2105-9-165.View ArticleGoogle Scholar
- Fischer A, Greenman C, Mustonen V. Germline fitness-based scoring of cancer mutations. Genetics. 2011; 188(2):383–93. doi:10.1534/genetics.111.127480.View ArticlePubMedPubMed CentralGoogle Scholar
- Youn A, Simon R. Identifying cancer driver genes in tumor genome sequencing studies. Bioinf (Oxford, England). 2011; 27(2):175–81. doi:10.1093/bioinformatics/btq630.View ArticleGoogle Scholar
- 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.View ArticlePubMedPubMed CentralGoogle Scholar
- Leiserson MDM, Blokh D, Sharan R, Raphael BJ. Simultaneous Identification of Multiple Driver Pathways in Cancer. PLoS Comput Biol. 2013; 9(5):1003054. doi:10.1371/journal.pcbi.1003054.View ArticleGoogle Scholar
- Khurana E, Fu Y, Chen J, Gerstein M. Interpretation of Genomic Variants Using a Unified Biological Network Approach. PLoS Comput Biol. 2013; 9(3):1002886. doi:10.1371/journal.pcbi.1002886.View ArticleGoogle Scholar
- Watson IR, Takahashi K, Futreal PA, Chin L. Emerging patterns of somatic mutations in cancer. Nat Rev Genet. 2013; 14(10):703–18. doi:10.1038/nrg3539.View ArticlePubMedPubMed CentralGoogle Scholar
Copyright
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. 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.