Inferring clonal evolution of tumors from single nucleotide somatic mutations
Abstract
Highthroughput sequencing allows the detection and quantification of frequencies of somatic single nucleotide variants (SNV) in heterogeneous tumor cell populations. In some cases, the evolutionary history and population frequency of the subclonal lineages of tumor cells present in the sample can be reconstructed from these SNV frequency measurements. However, automated methods to do this reconstruction are not available and the conditions under which reconstruction is possible have not been described.
We describe the conditions under which the evolutionary history can be uniquely reconstructed from SNV frequencies from single or multiple samples from the tumor population and we introduce a new statistical model, PhyloSub, that infers the phylogeny and genotype of the major subclonal lineages represented in the population of cancer cells. It uses a Bayesian nonparametric prior over trees that groups SNVs into major subclonal lineages and automatically estimates the number of lineages and their ancestry. We sample from the joint posterior distribution over trees to identify evolutionary histories and cell population frequencies that have the highest probability of generating the observed SNV frequency data. When multiple phylogenies are consistent with a given set of SNV frequencies, PhyloSub represents the uncertainty in the tumor phylogeny using a “partial order plot.” Experiments on a simulated dataset and two real datasets comprising tumor samples from acute myeloid leukemia and chronic lymphocytic leukemia patients demonstrate that PhyloSub can infer both linear (or chain) and branching lineages and its inferences are in good agreement with ground truth, where it is available. PhyloSub can be applied to frequencies of any “binary” somatic mutation, including SNVs as well as small insertions and deletions. The PhyloSub and partial order plot software, and supplementary material are available from http://morrislab.med.utoronto.ca/phylosub/.
Background
Cancer is a complex disease often associated with a characteristic series of somatic genetic variants [1, 2]. Substantial effort has been devoted to genetic profiling of tumors in hopes of identifying these driver mutations and studying how they drive tumor development and resistance to treatment[3]. Tumors often contain multiple, genetically diverse subclonal populations of cells[4], and in some cases it is possible to reconstruct the evolutionary history of the tumor, thereby aiding in the identification of driver mutations, by computing the population frequencies of mutations that distinguish the subclonal populations [5, 6, 7, 8, 9, 10, 11, 12, 13].
Somatic mutations can be detected, and roughly quantified, using exome and whole genome sequencing of a sample from a bulk tumor[14]. However, recent attempts to reconstruct subclonal phylogenies have employed much deeper targeted sequencing [15] of tumorassociated single nucleotide variants (SNVs) to achieve higher accuracy in estimated SNV frequencies[10, 16, 17, 9]. These SNV frequencies were then used to partially reconstruct the evolutionary history of tumors based on a single[10, 16] or multiple [9] samples of same tumor. However, due to short read sequencing, the frequencies of different SNVs are measured independently, so linkage between the SNVs in subclones is unavailable and standard phylogenetic methodology can not be used to construct evolutionary histories (as done in [18] or [17]). However, if one makes the infinite sites assumption about tumor evolution, namely that every SNV only appeared once, then it is possible to use SNV frequencies to automatically reconstruct full or partial subclonal phylogenies while also inferring the multiple SNV genotypes of the major subclonal lineages in the tumor.
Here we describe a new method that automatically performs this phylogenetic reconstruction. First, we demonstrate that an unambiguous reconstruction is possible by describing topological constraint rules that are sufficient conditions to infer whether a triplet of SNV frequencies is consistent with only a chain or a branching phylogeny. We then describe a new method, PhyloSub, that automatically infers tumor phylogenies from SNV allele frequencies measured in single or multiple tumor samples. PhyloSub is based on a generative probabilistic model, inference in which implicitly implements the two rules by inferring the hidden phylogenies that have high probabilities of generating the observed SNV frequencies. It uses Bayesian inference, based on Markov Chain Monte Carlo (MCMC) sampling, to infer a distribution over phylogenies that incorporates uncertainty due to multiple phylogenies being consistent with the SNV frequencies and also noise in the measurement of the SNV frequencies. PhyloSub uses a Dirichlet process prior over phylogenies [19] to group SNVs into major subclonal lineages.
Model assumptions
We assume that the tumor evolution proceeds according to the clonal evolution theory, namely that all tumor cells are derived from ancestors that gain growth advantages over normal tissue and begins to expand [18]. Subsequent mutations can provide a further fitness or survival advantage to their subclonal lineage[20] which subsequently increases in frequency compared to cells containing only the SNVs in the parental lineage. A given tumor sample is a snapshot of this evolutionary process and may contain, at nonnegligible frequency, cells from multiple major subclonal lineages, each containing a different assortment of these advantageous mutations. We make the infinite sites assumption [21, 22], namely that each SNV appears only once and furthermore that once it appears, it does not revert back to its original state. As we describe below and illustrate in Figure 1, in some circumstances, this assumption highly constrains the phylogenies that are consistent with the SNV allele frequency data, especially if SNV frequencies from multiple samples from the same tumor are available. Finally, to make our model robust to low tumor cellularity, we assume that each tumor is derived from a single clone, however, this assumption is not critical in modeling tumor evolution and we revisit this assumption in the discussion section where we describe how to generalize our model to multicentral tumors (e.g., [23]).
To simplify our initial discussion, we will assume that the exact population frequencies of the cells containing each SNV (i.e., the SNV population frequency) are available before discussing how we estimate these frequencies from deep sequencing data of the SNV locus. Note, we assume that the copy number of a locus is available as per [10]. In the datasets that we considered, most SNVs are heterozygous at a normal copy number locus and the population frequency of other SNVs is easily inferred from their allele frequencies. In more complex situations, a number of tools are available to infer copy number changes associated with specific subclonal lineages from whole genome sequencing data[11, 13].
An important consequence of the infinite sites assumption is that if SNV B occurred in a cell that contained SNV A, then all cells that have B also have A and thus the population frequency of A must always be greater than or equal to that of B, regardless of where and when the tumor sample was taken. However, a given set of three SNV population frequencies can still be consistent with two different phylogenies: a linear phylogeny or a branching phylogeny (see Figure 1A).
Topological constraint rules
One can distinguish linear or branching descent under some circumstances. For example, if we have already established that SNV A is ancestral to both B and C (i.e., that all cells with B or C also contain A), then if the population frequency of B plus the population frequency of C is greater than the population frequency of A, then the phylogeny must be linear. This is true because in a branching phylogeny, there are no cells that contain both B and C, so the population frequency of A must be at least as large as sum of the frequencies of B and C (see Figure 1B). We call this the “sum rule”. However, because a linear phylogeny is consistent with any set of SNV frequencies from a single sample, without making any further assumptions about the tumor evolution process, one needs at least two tumor samples to be able to rule out a linear phylogeny. However, given two samples and again assuming that SNV A is ancestral to both B and C, if the population frequency of B is larger than that of C in one sample, and vice versa in the other, than neither B nor C can be ancestral to the other, and the only phylogeny consistent with both sets of SNV frequencies is the branching one. We call this the “crossing rule” because the frequencies of B and C cross (see Figure 1C for an example). However, there is no guarantee that one can apply either rule to any set of SNV frequencies for all triplets of SNVs, although increasing the number of tumor samples does make it more likely that either the sum or crossing rule will be applicable for one or a pair of tumor samples, respectively. Furthermore, one needs to also consider the possibility of estimation error in the SNV population frequencies because these are inferred from discrete read counts. Note that these two rules also apply where SNV A is a mock SNV representing the wildtype state and having population frequency of ; as such these two rules also apply for multicentral tumors.
The PhyloSub algorithm
To explicitly model uncertainty in estimates of the SNV population frequencies and the precise tumor phylogeny, we have developed the PhyloSub model that we describe here. PhyloSub attempts to explain the observed read counts in terms of a latent phylogeny that associates SNVs with particular subclonal lineages. PhyloSub takes as input a set of read counts for a set of SNVs and the copy number status of each SNV, and estimates the number of major subclonal lineages, the mutational profile of each lineage, and the proportion of each lineage within the tumor cell population from which the read data was drawn. PhyloSub implements the parsimony assumptions detailed above using a nonparametric prior over tree structures. It is “generative” in that it attempts to explain the observed SNV frequencies in terms of an unobserved phylogeny; our model is also “Bayesian” in that it infers a posterior distribution over phylogenies and associated subclonal lineage frequencies. We introduce a new visualization, the partial order plot, to represent the posterior uncertainty in the phylogeny when the SNV frequencies alone do not provide sufficient information to uniquely reconstruct the phylogeny (Figure 2). The sum and crossing rule described above are implicitly incorporated into our generative model – our model assigns very low probability to any read counts that reflect deviations from either rule.
In the following, we provide a brief introduction to the PhyloSub model (see Methods for the full model) and we demonstrate its application to datasets where a single sample is profiled [17] and those where multiple samples are profiled [9]. We also report the application of the model on a simulated dataset to show that its prior parameterization allows it to represent a wide variety of phylogenies.
Results and Discussion
PhyloSub
PhyloSub represents the major subclonal lineages and their evolutionary relationships using a directed tree in which each node is associated with a subclonal lineage and the edges connect parental lineages to their direct child lineage. Each subclonal lineage is associated with a distinct subset of the SNVs input to the model, we call this subset the genotype of the lineage. Each node is also associated with (i) a set of SNVs that are present in this lineage but not its parent lineage and (ii) the population frequency of cells with the lineage genotype (and with no other SNVs from the input set). A subclonal lineage contains all of the SNVs associated with its parent, so its full genotype can be reconstructed by taking the union of the SNVs associated with its node and all of its ancestral nodes. Similarly, the population frequency of an SNV is the sum of the subclonal lineage frequencies of the lineage it appeared in and all of its descendent lineages. So, the subclonal lineage tree can be used to compute the population frequencies of each SNV and the genotype of each subclonal lineage. Associated with each SNV is a variable that indicates its zygosity and copy number in the cells that it appears (e.g., Aa indicating heterozygous and normal copy number), we assume all cells with the SNV have the same zygosity and copy number, and that all other cells have normal copy number at the SNV locus. The SNV genotype variable along with the population frequency is used to compute the allelic frequency of the SNV , . The data input to the model for each SNV is the number of reads mapping to the SNV locus, , and the number of these reads that do not contain the SNV, . We evaluate the likelihood of a given subclonal lineage tree (including the lineage population frequencies and the SNV genotype variables) by taking the product of the read count probabilities for each SNV, where the probability for the locus of SNV is computed using a binomial distribution whose parameter is derived from and an estimate of the error rate of the sequencer. PhyloSub also contains a vague prior over tree structures that is parameterized by three hyperparameters (see Methods) that govern how the prior scores trees with more or fewer nodes, and different average numbers of siblings. We use ranges for these hyperparameters that in simulations have a slight preference for trees with fewer nodes but a limited preference for sibling numbers (see below for details).
Simulations
PhyloSub’s Dirichlet process prior over tree structures depends on three hyperparameters: , , and . The hyperparameters and determine the number of nodes (subclones) in the tree, also affects the height of the tree and affects the number of siblings in the tree which in turn affects the width of the tree. In all the experiments, we sample these hyperparameters [19] as part of the MCMC sampling from a range whose upper and lower bounds we establish in this section.
To establish the ranges that we use for the hyperparameters in PhyloSub, we simulated read counts from clusters with an average of nine SNVs per cluster with SNV population frequencies , with a read depth of X which is a typical read depth for the targeting deep sequencing data that PhyloSub is designed for. We simulated heterozygous SNVs at loci with normal copy number and sample read counts for each SNV from a Binomial distribution (see Methods). The hyperparameter settings we used in the simulations are all possible combinations of , and . The SNV population frequencies are consistent with many different tree structures and Figure 3 shows that the tree structures with highest completedata likelihoods varies in the expected way for different settings of the tree prior hyperparameters. Although the preferred structure varies, the inferred SNV frequencies remain wellcorrelated with the baseline values (Pearson correlation ) for these hyperparameter ranges, so the prior is not overregularizing the SNV frequencies for these settings. To allow a range of tree structures, we integrate over these ranges by placing a uniform prior on the choice of these settings in our MCMC simulations (c.f., [19]).
Although we focused on high read depths in the above simulation, we found that PhyloSub works well for read depths X and was able to recover the clusters similar to the ones reported above and the SNV frequencies are wellcorrelated with the baseline values (Pearson correlation ). However, we found that the performance of the model degrades slightly at a read depth of X, due to merging of clusters whose nearby SNV frequencies could not be distinguished. Nonetheless, we note that the inferred SNV population frequency estimates remain wellcorrelated (Pearson correlation ) and that the majority of the clusters were recovered at read depth X.
The simulation as described above has no clear phylogeny by design. The SNV frequencies were consistent with multiple phylogenies and the main goal of this simulation was to establish the ranges for our hyperparameters that permit a wide variety of tree structures. We integrate over these parameter ranges on the real data in order to remove any prior bias towards particular structures. To determine whether PhyloSub can correctly recover the phylogenies from a single sample of SNV frequencies, we simulated read data from a chain phylogeny with SNV population frequencies . By the sum rule, these frequencies are only consistent with a chain phylogeny. PhyloSub was able to reliably recover this chain. The real datasets described in the later sections are representative of the types of problems that our methodology could be applied to as they contain single and multiple samples, some of which have clear phylogenies and some do not.
Results on AML datasets
To assess PhyloSub on single samples, we applied it to data from Jan et al. [17] who reported the coexistence of multiple subclonal lineages in hematopoietic stem cells (HSC) from acute myeloid leukemia (AML) patient samples. The deeptargeted sequencing of all SNV candidates identified by exome sequencing identified SNVs with differing allelic frequency, suggesting multiple clonal populations in the HSC cells. An independent singlecell assay confirmed the existence of multiple clones, and thus provides a ground truth tree that shows some of the major subclonal lineages within the populations. Here we apply PhyloSub to the two samples profiled by Jan et al. that had three or more SNVs profiled in a singlecell assay. These samples are SU048 and SU070 which have 6 and 10 SNVs in the singlecell assay, respectively. Although this assay confirmed the presence of some of the subclonal lineages, only 100200 cells were assayed, so lineages with low population frequency in the sample (e.g., ) may not be detected.
We applied PhyloSub providing it with the copy number and zygosity of each SNV (results were similar if we assume normal copy number and have a uniform prior on zygosity). For both SU048 and SU070, a number of different phylogenies were consistent with the SNV read counts, and we developed the “partial order plot” to represent the posterior uncertainty in the phylogeny (see Figure 2 and Methods).
Figure 4 shows that partial order plot for SU070. The ordering of the nodes in the partial order plot can also be used to infer ancestry via transitivity, for example, in Figure 4, the SNV CXorf66 has high probability of being in the subclonal lineage that is the direct parent of the one that DOCK9 is in, however, because the TET2T1884A SNV is sorted before CXorf66 (and has a small probability of being a direct parent of it), then in the PhyloSub posterior over lineages, TET2T1884A has a high probability of being in an ancestral lineage to the one CXorf66 is in.
Furthermore, one can interpret the partial order plot to indicate that both CXorf36 and CXorf66 are in the same lineage because they are both direct parents of DOCK9 (with high probability) and there are no edges between them. For reference, in Figure 4 we have included the results of the singlecell assay for SU070 in the partial order plot representation – Jan et al. report three subclonal lineages for SU070, as indicated by the SNV colorings[17]. We note that these plots are largely consistent. Indeed, we assign high posterior probability , to two of the three subclonal lineages detected by Jan et al. (see Table S3 in Additional file 2 for full lineage genotype probabilities). For reference, we also provide the list of the subclonal lineage trees along with their posterior probabilities in Table S1 (see Additional file 2).
The one major difference between PhyloSub’s estimates and the singlecell data from Jan et al. is that PhyloSub switches the order of the appearance of SNVs CXorf36 and TET2T1884A. In fact, there was not a single subclonal lineage that contained CXorf36 but not TET2T1884A in 5,000 subclonal lineage trees sampled from PhyloSub’s posterior. This switch is likely due to the observed SNV frequencies, indeed the 95% confidence intervals of the SNV frequencies of these two SNVs do not overlap (Table 1).
SNV  Variant allele read counts  Read depth  Allele frequency  Cluster ID 
CACNA1H 
12,085  24,860  0.486 (95 CI: 0.4810.491)  A 
TET2T1884A  4,220  8,772  0.481 (95 CI: 0.4720.490)  B 
TET2Y1649stop  7,792  16,211  0.481 (95 CI: 0.4740.487)  A 
CXorf66  3,684  8,150  0.452 (95 CI: 0.4430.461)  B 
CXorf36  3,523  8,060  0.437 (95 CI: 0.4280.446)  A 
DOCK9  3,391  8,676  0.391 (95 CI: 0.3820.400)  C 
NCRNA00200  9,201  25,413  0.362 (95 CI: 0.3570.367)  C 
CTCF  10,558  30,119  0.351 (95 CI: 0.3460.355)  C 
GABARAPL1  1,648  4,992  0.330 (95 CI: 0.3190341)  C 
SCN4B  5,113  16,386  0.312 (95 CI: 0.3060.318)  C 

One explanation for this difference is a systematic bias in the measurement of one or both of these SNVs; it is also possible that the labels of these two SNVs were switched in Jan et al.. We also note, however, that in Jan et al., the existence of the lineage that contains only CXorf36, TET2Y1649stop, and CACNA1H is only supported by 2 of the 189 clones that they profiled.
For the tumor sample SU048, both the partial order plot and the singlecell assay agree on TET2E1357stop event occurring early (at the root of the tree), and all other SNVs are secondary mutational events as shown in Figure 5B. Note that the partial order plot shows a large uncertainty in the structure for the rest of the SNVs and this is also reflected in the posterior over subclonal lineage trees and genotypes (see Tables S2 and S4 in Additional file 2, respectively). There is no strong evidence for either a linear or branching lineage or for particular clustering among these SNVs. Also, from Table 2, we see a lot of variation in the allele frequencies of these SNVs suggesting that they may not belong to the same subclonal lineages.
SNV  Variant allele read counts  Read depth  Allele frequency  Cluster ID 

TET2E1357stop  7,436  19,553  0.380 (95 CI: 0.3750.386)  A 
SMC1A  182,974  660,069  0.277 (95 CI: 0.2760.278)  B 
ACSM1  17,149  127,236  0.135 (95 CI: 0.1330.136)  B 
OLFM2  13,828  122,523  0.113 (95 CI: 0.1110.114)  B 
TET2D1384V  1,833  17,687  0.104 (95 CI: 0.1000.107)  B 
ZMYM3  18,536  307,346  0.060 (95 CI: 0.0600.061)  B 
The subclonal lineage inferred by Jan et al.’s singlecell assay is shown in Figure 5A and only contains two lineages, one with only TET2E1357stop and the other with the other five SNVs. The TET2E1357stop lineage genotype has probability in our posterior, however the second genotype has a relatively small probability () under the posterior although we note that the genotype that contains all SNVs but ZMYM3 has a posterior probability of (see Table S4 in Additional file 2). For reference, we also provide the list of the subclonal lineage trees along with their posterior probabilities in Table S2 (see Additional file 2).
In summary, the subclonal lineage trees inferred by PhyloSub on single samples of SNV frequencies are largely consistent with ground truth but there remains substantial uncertainty in SU048 about whether there was a linear or chain lineage. On the other hand, the SNV frequencies in SU070 were only consistent with a linear lineage and PhyloSub almost perfectly reconstructed the results of the singlecell assay with one misordering of the SNVs. This difference may be explained by unmodeled systematic biases in the deep sequencing data or experimental error. Nonetheless, we have shown that in some cases, it is possible to achieve a good estimate of the genotype of multiple subclonal lineages as well as their evolution from a single, targeted deep sequencing sample of SNV frequencies.
Results on CLL datasets
To evaluate PhyloSub on a multiple sample dataset, we used data from a study of chronic lymphocytic leukemia (CLL) by Schuh et al. [9] which quantified SNV frequencies of a set of SNVs during different time points spanning the patient therapy cycle. The candidate SNVs were identified by exome sequencing and then subjected to targeted resequencing. The tumor samples from the three patients in the study labeled CLL077, CLL006 and CLL003 have 11, 16 and 20 SNVs respectively with SNV frequencies for five different time points. Originally, Schuh et al. reconstructed the evolutionary histories of each cancer by a semimanual procedure in which they first automatically grouped SNVs into subclonal lineages using means clustering on the allele frequencies and the differences in allele frequencies between the time points for each patient and then reconstructed the evolutionary structure of those lineages using a procedure that they do not describe in the paper. In PhyloSub, we model multiple samples from the same cancer as sharing the same evolutionary history but we allow subclonal frequencies to change between samples.
We applied PhyloSub to the SNV read count data, providing the algorithm with the likely zygosity estimates – in most cases, SNVs appeared to be heterozygous with normal copy number but in a few cases, SNVs appeared to be hemizygous and were input to the model as such. For these data, because of the multiple samples per tumor, there is very little posterior uncertainty in the best fitting tree – as such, we only show the best single tree structure corresponding to the MCMC sample with the highest completedata likelihood [19].
For the tumor samples CLL077 and CLL003, the best tree structure estimated by PhyloSub and the tree structure from Schuh et al. [9] are in exact agreement and the population frequencies of the subclonal lineages are wellcorrelated. Figures 6 and 8 compare the PhyloSub estimates with those reported by Schuh et al. [9].
For the tumor sample CLL006, PhyloSub inferred a chain structure similar to the chain structure from Schuh et al., but the major difference in PhyloSub’s best estimate of the tree structure is the splitting of cluster A into two clusters as shown in Figure 7. However, we found that the completedata log likelihood of PhyloSub’s best estimate of the tree structure is higher than the one for the chain structure of Schuh et al. and therefore PhyloSub prefers the splitting of the cluster A into two clusters.
In the CLL dataset, there is no ground truth but to allow the reader to compare the two estimates of the evolutionary history, Figure 9 plots the frequency of each SNV in the three samples, and we have colored SNVs according to their subclonal lineage assignments by Schuh et al.. These SNV frequencies are not corrected for copy number, however, the hemizygous SNVs are clear from examination of the figure.
In summary, having multiple samples of SNV frequencies greatly reduces the posterior uncertainty in the evolutionary history of the tumor and PhyloSub is able to reconstruct histories produced by a semimanual procedure.
Conclusions
We presented a nonparametric Bayesian model called PhyloSub that uses a Dirichlet process prior over trees [19] to model the clonal evolutionary structure of tumors from next generation sequencing data. We also introduced a new visualization method, the partial order plot, to represent the posterior uncertainty in the phylogeny when the clonal frequencies alone do not provide sufficient information to uniquely reconstruct the phylogeny and mutational profiles of each subclonal lineage represented in the tumor. By enforcing a set of structural constraints on the SNV population frequencies using MCMC methods, we were able to infer the phylogenetic relationships between subclones from both single and multiple tumor samples.
We have demonstrated that it is possible, in some cases, to detect a linear lineage from a single, high cellularity sample of the tumor. We have also shown that multiple samples highly constrain the possible lineages that are consistent with the SNV frequency data. PhyloSub’s inferred subclonal lineage trees were in good agreement with single cell assays on single sample data and with an expertdriven, semimanual reconstruction procedure on multiple sample data.
PhyloSub’s ability to detect and characterize subclonal lineages depends on the frequency of the lineage in the population (compared to its descendant lineages), the number of SNVs that define the lineage, as well as the accuracy with which the SNV population frequencies are estimated which depends on both the sequencing depth as well as uncertainty about the copy number of the SNV. Simply put, for lineages defined by a single SNV, the read depth has to be high enough that the uncertainty in the estimated SNV frequency is less than the frequency of the subclonal population. Having more lineagedefining SNVs can relax this hard constraint. As such, the phylogenies of tumors with large numbers of subclonal lineages, each defined by a small number of SNVs (possibly due to a pronounced hypermutability phenotype), will be hard to reconstruct with PhyloSub, or any other method, unless the SNV frequencies are very accurately estimated. Indeed, it is not clear how ground truth could be uncovered in such a case: the gold standard of single cell sequencing would require an exceptionally large number of single cells to survey this highly heterogeneous population, and each of these cells would need to be sequenced deeply in order to ensure precise somatic variant calling.
One potential difficulty in scaling our approach to orders of magnitude more SNVs is that the Markov chain may not mix in a timely manner, in other words, may get stuck in local minima. We note that finding suboptimal solutions is an issue for any method based on these data. In our case, the mixing time of the chain would depend largely on the number of subclones represented in the population with less of a dependence on the number of SNVs. There are various techniques for determining whether or not a Markov chain is wellmixed and we refer the reader to a recent excellent review [25].
PhyloSub extends recent work on inferred cellularity and subclonal structure from somatic mutations. ABSOLUTE uses whole genome sequence data or array CGH data to identify regions of copy number change in the tumor and based on this infers cellularity and copy number changes associated with different subclones [11]. THetA [13] also attempts to infer both the copy number profiles and their relative proportions using the whole genome sequencing data based on an infinite sites assumptions. Neither of these algorithms explicitly reconstructs tumor phylogenies. Our work is closest to PyClone [10] which uses a flat Dirichlet process mixture model to group SNVs into subclonal lineages based on their frequencies; PhyloSub extends this work by reconstructing the phylogenetic relationships among these lineages and, in doing so, allows the full SNV genotype of each subclonal population to be reconstructed.
We designed PhyloSub to assume a single clonal origin for the cancerous cells in the sample. We made this decision to increase the applicability of the sum rule for low cellularity tumors (i.e., tumors with high normal contamination). However, removing this assumption would be a simple change to the model, which we have not evaluated.
Another area of future innovation would be in modeling sequencing biases and uncertainty in SNV allele frequencies resulting from them. We did not evaluate replacing our binomial model with a negative binomial one that would have allowed greater variability in the observed read counts for a given SNV allele frequency [26].
Methods
Dirichlet process mixture models
Consider the problem of clustering objects using a Bayesian finite mixture model of components (clusters) with the following generative process [27]:
(1) 
where are the mixing weights such that , is the concentration parameter of the symmetric Dirichlet prior placed on the mixing weights, is the cluster assignment variable, is the prior distribution from which the component parameters are drawn, is the component distribution parameterized by . The finite mixture model can be extended to a model with an infinite number of mixture components by replacing the Dirichlet prior with a Dirichlet process (DP) prior resulting in what is known as the DP mixture model (DPMM) [28]. Unlike finite mixture models, DPMMs automatically estimate the number of components from the data thereby circumventing the problem of fixing the number of components a priori. The stickbreaking construction [29] given below provides a precise recipe to draw samples from a Dirichlet process:
(2) 
where is a point mass centered at and , i.e., is a draw from a DP with base distribution and concentration parameter . The stickbreaking process can be viewed as recursively breaking sticks of length , starting with a stick of unit length. The beta variates determine the random location at which the stick is broken. The concentration parameter determines the number of clusters with high values resulting in large number of clusters. Let denote the stickbreaking process over . Replacing the Dirichlet prior in the finite mixture model (1) with the stickbreaking process prior results in the following generative process for infinite mixture models:
An alternative view of the above generative process produces component parameters by drawing samples from resulting in the following generative process:
(3) 
Note that in the above process every object is associated with a component parameter and that all objects assigned to the same cluster will have the same component parameter. In other words, multiple elements in the set will take on the same value from the set of unique parameters.
Treestructured stickbreaking process
The stickbreaking construction (2) described above can be used to produce a flat clustering of objects, where the clusters are independent of each other. Adams et al. [19] extended this construction for hierarchical clustering by interleaving two stickbreaking processes. This construction results in a relational clustering of objects where the clusters are connected to form a rooted tree structure. Unlike classical hierarchical clustering algorithms such as agglomerative clustering, this construction allows data to reside in the internal nodes of the tree; a feature we exploit to model the association of SNVs with subclonal lineages.
We borrow notation from Adams et al. [19]. Let denote a sequence of positive integers used to index the nodes of the tree. Let denote the zerolength string, i.e., the root of the tree. Let indicate the length of the sequence and therefore the depth of node . Let denote the sequence formed by appending to . The children of node is the set and let the ancestors of be denoted by the set . The interleaved, twolayered stickbreaking construction is as follows:
(4)  
The and determine the amount of mass allocated to and its descendants respectively, whereas determines the probability of a particular sequence of children. The construction ensures that the mixing weights sum to one. The parameters and control the height and the width of the tree respectively. Note that the concentration parameter is a function of the depth of the tree () and is defined to be with and [19].
PhyloSub model
We follow Shah et al. [30, 10] to model the allelic count data. For each genetic variant that is detected by highthroughput sequencing methods, cells containing the genetic variant is referred to as variant population and those without the variant as reference population. Let denote the set of nucleotides. Let and denote the number of reads matching the reference allele and the variant allele respectively at position , and let . The genotype would depend on the copy number at the variant location. Let denote the probability of sampling a reference allele from the reference population. This value depends on the error rate of the sequencer. Let denote a vector whose entries, , are the probabilities of sampling a reference allele from the variant population with genotype at position . Let denote the vector whose entries, , are the probabilities of the variant population at position to have the genotype . Let denote the pseudocount parameters of the Dirichlet prior over . Let denote the genotype of the variant population at position . Let denote the fraction of cells from the variant population, i.e., the SNV population frequency at position , and denote the fraction of cells from the reference population at position . The observation model for allelic counts has the following generative process [10]:
(5)  
The posterior distribution of is
Each of the terms appearing inside the summation over genotypes is the probability distribution of a Dirichlet compound multinomial (with a single draw) [31]. The posterior distribution can thus be rewritten as
(6) 
where is the Gamma function.
The Dirichlet process prior in the observation model of allelic counts (5) is useful to infer groups of mutations that occur at the same SNV population frequency [10]. Furthermore, being a nonparametric prior, it is useful to avoid the problem of selecting the number of groups of mutations a priori. However, it cannot be used to model the clonal evolutionary structure which takes the form of a rooted tree. In order to model this, we propose to use the treestructured stickbreaking process prior (4) described in the previous section. The probabilistic graphical model for allelic counts with the treestructured stickbreaking process prior is shown in Figure 10.
Inputs to the model including the hyperparameters are indicated in shaded nodes, whereas the latent variables including the set of SNV population frequencies are indicated in unshaded nodes. The prior/base distribution of the SNV population frequencies is the uniform distribution for the root node and for any other node in the tree, where denotes the parent node of and is the set of siblings of . This ensures that the clonal evolutionary constraints (discussed in the next section) are satisfied when adding a new node in the tree. The crucial difference between our model and the model of Shah et al. [10] is that we use a treestructured stickbreaking process instead of a Dirichlet process (cf. (5)) to generate the set of SNV population frequencies . Given this model and a set of observations/inputs , the tree structure and the SNV population frequencies are inferred using Markov chain Monte Carlo sampling. In particular, we use Gibbs sampling [19] to generate posterior samples of the SNV population frequencies (6). Each iteration of the Gibbs sampler involves multiple subsampling procedures: sampling cluster assignments , sampling stick lengths and , sampling stickbreaking hyperparameters , and , and sampling the SNV population frequencies . Our main algorithmic contribution, described below, is a method to sample SNV population frequencies in such a way that the tumor evolution proceeds according to the assumptions from the clonal evolutionary theory. The rest of the subsampling procedures follow directly from Adams et al. [19] and we refer the reader to it for further technical details.
Sampling SNV population frequencies
Given the current state of the tree structure, we sample SNV population frequencies in such a way that the SNV population frequency of every nonleaf node in the tree is greater than or equal to the sum of the SNV population frequencies of its children. To enforce this constraint, we introduce a set of auxiliary weights , one for each node, that satisfy , and rewrite the observation model for allelic counts (5) explicitly in terms of these weights resulting in the following posterior distribution:
(7) 
where we have used to denote the auxiliary weights for each SNV. The prior/base distribution of the auxiliary weights is defined such that it is 1 for the singleton root node and for any other node in the tree, where denotes the parent node of . When a new node is added to the tree, we sample and update . This ensures that .
This change is crucial as it allows us to design a Markov chain that converges to the stationary distribution of . The SNV population frequency for any node can then be computed via
(8) 
where and are the sets of all descendants and children of node respectively. This construction ensures that the SNV population frequencies of mutations appearing at the parent node is greater than or equal to the sum of the frequencies of all its children. The procedure to generate a random sample of SNV population frequencies is given in Algorithm 1 where we generate for every node by traversing the tree in a breadthfirst fashion. The input to this algorithm is the current state of the tree where is the set of vertices and is the set of edges, and the output is a multidimensional sample of SNV population frequencies (where ) and the corresponding auxiliary weights . A sample from this algorithm is shown in Figure 11.
We use MetropolisHastings algorithm to sample from the posterior distribution of the auxiliary weights (7) as shown in Algorithm 2 and derive the SNV population frequencies from these samples. We use an asymmetric Dirichlet distribution as the proposal distribution. This ensures that the Markov chain converges to the stationary distribution of . The inputs to the sampling algorithm are the current state of the tree , a scaling factor , and the number of iterations . The output is a sample from the posterior distribution of and its corresponding .
Extension to multiple tumor samples
PhyloSub (cf. Figure 10) can be easily extended for multiple tumor samples. We allow the treestructured stickbreaking process prior (4) to be shared across all the samples. Let and denote the number of reads matching the reference and the variant allele respectively at position for sample , and let . Let denote the fraction of cells from the variant population, i.e., the SNV population frequency at position for sample , and denote its corresponding auxiliary weight. The graphical model of PhyloSub for multiple tumor samples is shown in Figure 12.
The main technical difference between the single and the multiple sample models lies in the sampling procedure for SNV population frequencies. In the multiple sample model, we ensure that the clonal evolutionary constraints described in the previous section are satisfied separately for each tumor sample and then make a global MetropolisHastings move based on the distribution , where is the set of auxiliary weights for all the tumor samples.
Partial order plot
We construct a partial order plot to summarize and visualize the trees from all the posterior MCMC samples. It is important to note that the nodes of this partial order plot are the SNVs and not the SNV clusters. The thickness of a directed edge in the tree is proportional to the fraction of MCMC samples in which SNV first appears in a subclonal lineage that was the parent of the subclonal lineage that first appears in. The color of the border of the SNVs represents the subclonal lineage cluster that the SNV is placed into posthoc using an algorithm called correlation clustering [32]. Note that the main purpose of this clustering algorithm is only to color the nodes in the partial order plot by aggregating the clustering information from all the MCMC samples obtained from our model; this clustering is a summary but does not represent any (possibly quite large) uncertainty in the cluster assignments. The input to this algorithm is a symmetric coclustering matrix C, whose elements is the difference between the number of samples in which and were assigned to the same SNV cluster and the number of samples in which and were assigned to different SNV clusters. The algorithm estimates a symmetric cluster indicator matrix , whose elements if and are assigned to the same cluster and otherwise. This cluster indicator matrix has all the information about the number of SNV clusters as well as the SNVs assigned to each of them.
MCMC settings
In all the experiments, we fix the number of MCMC iterations to 5,000 with a burnin of 100 samples. We also fix the number of iterations in the MetropolisHastings algorithm to 5,000 and set the scaling factor for the Dirichlet proposal distribution to . We run the MCMC samplers multiple times with different random initializations and pick a single run based on the completedata likelihood trace and its autocorrelation function. We use all the 5,000 samples without thinning [33] to construct the partial order plots. We use the CODA R package [34] for MCMC diagnostics to monitor the convergence of the samplers. The completedata log likelihood traces and the corresponding autocorrelation function plots after the burnin period of 100 samples for all the experiments on AML and CLL datasets are shown in Figures S3 to S7 (Additional file 1).
Datasets and inputs to PhyloSub
All datasets used in the experiments, including details about the inputs to PhlyoSub, i.e., the set of observations , are provided in the Additional file 2 (Tables S5 – S10).
Acknowledgements
This work was funded by a National Science and Engineering Research Council (NSERC) operating grant and a Early Researcher Award from the Ontario Research Fund to QM. WJ and LS were supported by The Ministry of Research and Innovation, Province of Ontario.
We thank Andrew Roth and Sohrab Shah for sharing a prepublication version of the PyClone software with us; helping us to install it and to duplicate their published results; and for providing unpublished details of the PyClone observation model in their user manual.
References
 D. Hanahan and R. A. Weinberg. The hallmarks of cancer. Cell, 100(1):57–70, 2000.
 D. Hanahan and R. A. Weinberg. Hallmarks of cancer: The next generation. Cell, 144(5):646–674, 2011.
 T. A. Yap, M. Gerlinger, P. A. Futreal, L. Pusztai, and C. Swanton. Intratumor heterogeneity: Seeing the wood for the trees. Science Translational Medicine, 4(127):127ps10, 2012.
 J. E. Visvader. Cells of origin in cancer. Nature, 469(7330):314–322, 2011.
 N. E. Navin and J. Hicks. Tracing the tumor lineage. Molecular Oncology, 4(3):267–283, 2010.
 Charles G. Mullighan, Letha A. Phillips, Xiaoping Su, Jing Ma, Christopher B. Miller, Sheila A. Shurtleff, and James R. Downing. Genomic analysis of the clonal origins of relapsed acute lymphoblastic leukemia. Science, 322(5906):1377–1380, 2008.
 Marco Gerlinger, Andrew J. Rowan, Stuart Horswell, James Larkin, David Endesfelder, Eva Gronroos, Pierre Martinez, Nicholas Matthews, Aengus Stewart, Patrick Tarpey, Ignacio Varela, Benjamin Phillimore, Sharmin Begum, Neil Q. McDonald, Adam Butler, David Jones, Keiran Raine, Calli Latimer, Claudio R. Santos, Mahrokh Nohadani, Aron C. Eklund, Bradley SpencerDene, Graham Clark, Lisa Pickering, Gordon Stamp, Martin Gore, Zoltan Szallasi, Julian Downward, P. Andrew Futreal, and Charles Swanton. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. The New England Journal of Medicine, 366(10):883–892, 2012.
 A. Marusyk and K. Polyak. Tumor heterogeneity: Causes and consequences. Biochimica et Biophysica Acta, 1805(1):105–117, 2010.
 Anna Schuh, Jennifer Becq, Sean Humphray, Adrian Alexa, Adam Burns, Ruth Clifford, Stephan M. Feller, Russell Grocock, Shirley Henderson, Irina Khrebtukova, Zoya Kingsbury, Shujun Luo, David McBride, Lisa Murray, Toshi Menju, Adele Timbs, Mark Ross, Jenny Taylor, and David Bentley. Monitoring chronic lymphocytic leukemia progression by whole genome sequencing reveals heterogeneous clonal evolution patterns. Blood, 120(20):4191–4196, 2012.
 Sohrab P. Shah, Andrew Roth, Rodrigo Goya, Arusha Oloumi, Gavin Ha, Yongjun Zhao, Gulisa Turashvili, Jiarui Ding, Kane Tse, Gholamreza Haffari, Ali Bashashati, Leah M. Prentice, Jaswinder Khattra, Angela Burleigh, Damian Yap, Virginie Bernard, Andrew McPherson, Karey Shumansky, Anamaria Crisan, Ryan Giuliany, Alireza HeraviMoussavi, Jamie Rosner, Daniel Lai, Inanc Birol, Richard Varhol, Angela Tam, Noreen Dhalla, Thomas Zeng, Kevin Ma, Simon K. Chan, Malachi Griffith, Annie Moradian, S.W. Grace Cheng, Gregg B. Morin, Peter Watson, Karen Gelmon, Stephen Chia, SuetFeung Chin, Christina Curtis, Oscar M. Rueda, Paul D. Pharoah, Sambasivarao Damaraju, John Mackey, Kelly Hoon, Timothy Harkins, Vasisht Tadigotla, Mahvash Sigaroudinia, Philippe Gascard, Thea Tlsty, Joseph F. Costello, Irmtraud M. Meyer, Connie J. Eaves, Wyeth W. Wasserman, Steven Jones, David Huntsman, Martin Hirst, Carlos Caldas, Marco A. Marra, and Samuel Aparicio. The clonal and mutational evolution spectrum of primary triplenegative breast cancers. Nature, 486(7403):617–656, 2012.
 Scott L Carter, Kristian Cibulskis, Elena Helman, Aaron McKenna, Hui Shen, Travis Zack, Peter W. Laird, Robert C. Onofrio, Wendy Winckler, Barbara A. Weir, Rameen Beroukhim, David Pellman, Douglas A. Levine, Eric S. Lander, Matthew Meyerson, and Gad Getz. Absolute quantification of somatic DNA alterations in human cancer. Nature Biotechnology, 30(5):413–421, 2012.
 Dan A. Landau, Scott L. Carter, Petar Stojanov, Aaron McKenna, Kristen Stevenson, Michael S. Lawrence, Carrie Sougnez, Chip Stewart, Andrey Sivachenko, Lili Wang, Youzhong Wan, Wandi Zhang, Sachet A. Shukla, Alexander Vartanov, Stacey M. Fernandes, Gordon Saksena, Kristian Cibulskis, Bethany Tesar, Stacey Gabriel, Nir Hacohen, Matthew Meyerson, Eric S. Lander, Donna Neuberg, Jennifer R. Brown, Gad Getz, and Catherine J. Wu. Evolution and impact of subclonal mutations in chronic lymphocytic leukemia. Cell, 152(4):714–726, 2013.
 Layla Oesper, Ahmad Mahmoody, and Benjamin J. Raphael. THetA: Inferring intratumor heterogeneity from highthroughput DNA sequencing data. Genome Biology, 14(7):R80, 2013.
 A. Marusyk, V. Almendro, and K. Polyak. Intratumour heterogeneity: A looking glass for cancer? Nature Reviews Cancer, 12(5):323–334, 2012.
 M. Meyerson, S. Gabriel, and G. Getz. Advances in understanding cancer genomes through secondgeneration sequencing. Nature Reviews Genetics, 11(10):685–696, 2010.
 Serena NikZainal, Peter Van Loo, David C. Wedge, Ludmil B. Alexandrov, Christopher D. Greenman, King Wai Lau, Keiran Raine, David Jones, John Marshall, Manasa Ramakrishna, Adam Shlien, Susanna L. Cooke, Jonathan Hinton, Andrew Menzies, Lucy A. Stebbings, Catherine Leroy, Mingming Jia, Richard Rance, Laura J. Mudie, Stephen J. Gamble, Philip J. Stephens, Stuart McLaren, Patrick S. Tarpey, Elli Papaemmanuil, Helen R. Davies, Ignacio Varela, David J. McBride, Graham R. Bignell, Kenric Leung, Adam P. Butler, Jon W. Teague, Sancha Martin, Goran Jönsson, Odette Mariani, Sandrine Boyault, Penelope Miron, Aquila Fatima, Anita Langerød, Samuel A.J.R. Aparicio, Andrew Tutt, Anieta M. Sieuwerts, Åke Borg, Gilles Thomas, Anne Vincent Salomon, Andrea L. Richardson, AnneLise BørresenDale, P. Andrew Futreal, Michael R. Stratton, and Peter J. Campbell. The life history of 21 breast cancers. Cell, 149:994–1007, 2012.
 Max Jan, Thomas M. Snyder, M. Ryan CorcesZimmerman, Paresh Vyas, Irving L. Weissman, Stephen R. Quake, and Ravindra Majeti. Clonal evolution of preleukemic hematopoietic stem cells precedes human acute myeloid leukemia. Science Translational Medicine, 4(149):149ra118, 2012.
 Peter J. Campbell, Erin D. Pleasance, Philip J. Stephens, Ed Dicks, Richard Rance, Ian Goodhead, George A. Follows, Anthony R. Green, P. Andy Futreal, and Michael R. Stratton. Subclonal phylogenetic structures in cancer revealed by ultradeep sequencing. Proceedings of the National Academy of Sciences, 105(35):13081–13086, 2008.
 Ryan Prescott Adams, Zoubin Ghahramani, and Michael I. Jordan. Treestructured stick breaking for hierarchical data. In Advances in Neural Information Processing Systems 23, pages 19–27, 2010.
 J.A. Brosnan and C.A. IacobuzioDonahue. A new branch on the tree: Nextgeneration sequencing in the study of cancer evolution. Seminars in cell and developmental biology, 23(2):237–242, 2012.
 Motoo Kimura. The number of heterozygous nucleotide sites maintained in a finite population due to steady flux of mutations. Genetics, 61(4):893, 1969.
 Richard R Hudson. Properties of a neutral allele model with intragenic recombination. Theoretical population biology, 23(2):183–201, 1983.
 T. M. Shattuck, W. H. Westra, P. W. Ladenson, and A. Arnold. Independent clonal origins of distinct tumor foci in multifocal papillary thyroid carcinoma. New England Journal of Medicine, 352(23):2406–2412, 2005.
 Graphviz  Graph Visualization Software. http://graphviz.org.
 David A. Levin, Yuval Peres, and Elizabeth L. Wilmer. Markov chains and mixing times. American Mathematical Society, 2008.
 Mark D Robinson and Alicia Oshlack. A scaling normalization method for differential expression analysis of RNAseq data. Genome Biology, 11(3):R25, 2010.
 Y. W. Teh. Dirichlet processes. In Encyclopedia of Machine Learning. Springer, 2010.
 Charles E. Antoniak. Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. Annals of Statistics, 2(6):1152–1174, 1974.
 J. Sethuraman. A constructive definition of Dirichlet process. Statistica Sinica, 4:639–650, 1994.
 Pyclone. http://compbio.bccrc.ca/software/pyclone/.
 Thomas P. Minka. Estimating a Dirichlet distribution. Technical report, Microsoft Research Cambridge,, 2012.
 Nikhil Bansal, Avrim Blum, and Shuchi Chawla. Correlation clustering. Machine Learning, 56(13):89–113, 2004.
 W. A. Link and M. J. Eaton. On thinning of chains in MCMC. Methods in Ecology and Evolution, 3(1):112–115, 2012.
 Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. CODA: Convergence diagnosis and output analysis for MCMC. R News, 6(1):7–11, 2006.