Inferring clonal evolution of tumors from single nucleotide somatic mutations

Background High-throughput 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. But automated methods to do this reconstruction are not available and the conditions under which reconstruction is possible have not been described. Results 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. Conclusions 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 is available from https://github.com/morrislab/phylosub/.


Abstract
High-throughput 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. Phy-loSub 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 tumor-associated single nucleotide variants (SNVs) to achieve higher accuracy in estimated SNV frequencies [9,10,16,17]. 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 non-negligible 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 100%; 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 non-parametric 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. Nodes correspond to SNVs. Edge thickness is proportional to the posterior probability that parent SNVs is in a subclonal lineage that is the parent of the one containing the child SNV. SNV nodes are ordered based on a topological sort implemented as the "layered graph drawing method" in Graphviz [24]. In this example, SNV A is always in a subclonal lineage that is the parent of B but the same is true with probability 0.5 for C.

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 i, p i . The data input to the model for each SNV is the number of reads mapping to the SNV locus, d i , and the number of these reads that do not contain the SNV, a i . 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 i is computed using a binomial distribution whose parameter is derived from p i 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 (α 0 , γ, λ) (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: α 0 , γ, and λ. The hyperparameters α 0 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 {1.0, 0.85, 0.6, 0.35, 0.2, 0.08}, with a read depth of ≈ 10, 000X 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 α 0 ∈ {1, 2, 4, 10, 20, 50}, γ ∈ {1, 2, 4, 6, 8} and λ ∈ {0.25, 0.5, 1}. The SNV population frequencies are consistent with many different tree structures and Figure 3 shows that the tree structures with highest complete-data likelihoods varies in the expected way for different settings of the tree prior hyperparameters. Although the preferred structure varies, the inferred SNV frequencies remain well-correlated with the baseline values (Pearson correlation > 0.99) for these hyperparameter ranges, so the prior is not over-regularizing 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 ≈ 1, 000X and was able to recover the clusters similar to the ones reported above and the SNV frequencies are well-correlated with the baseline values (Pearson correlation > 0.99). However, we found that the performance of the model degrades slightly at a read depth of ≈ 200X, due to merging of clusters whose nearby SNV frequencies could not be distinguished. Nonetheless, we note that the inferred SNV population frequency estimates remain well-correlated (Pearson correlation > 0.96) and that the majority of the clusters were recovered at read depth ≈ 200X.
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 0.9 → 0.75 → 0.55 → 0.4 → 0.25. 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 deep-targeted 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 single-cell 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 single-cell assay. These samples are SU048 and SU070 which have 6 and 10 SNVs in the single-cell assay, respectively. Although this assay confirmed the presence of some of the subclonal lineages, only 100-200 cells were assayed, so lineages with low population frequency in the sample (e.g., < 1%) 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 TET2-T1884A SNV is sorted before CXorf66 (and has a small probability of being a direct parent of it), then in the PhyloSub posterior over lineages, TET2-T1884A 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 single-cell 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 > 0.96, to two of the three subclonal lineages detected by Jan et al. (see Table S3 in Additional file 2 for full lineage genotype NCRNA00200 CTCF GABARAPL1 SCN4B Figure 4: Clonal evolutionary structures of tumor sample SU070. (A) Ground truth from Jan et al. [17]. (B) PhyloSub's output summarized using a partial order plot. For clarity, we removed edges with probability less than 0.1 before laying out the nodes and we only show SNVs for which single-cell data is available. However, all the SNVs whose frequencies were reported in this study were used in the inference. The color of the border of the SNVs represents the subclonal lineage cluster that the SNV is placed into by a graph-based clustering algorithm that takes as input the co-clustering frequencies from the MCMC samples (see Methods). Note that unlike the thickness of the edges, this is simply a visualization aid, and does not fully represent the model's posterior uncertainty in the SNV clusterings.
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 single-cell data from Jan et al. is that PhyloSub switches the order of the appearance of SNVs CXorf36 and TET2-T1884A. In fact, there was not a single subclonal lineage that contained CXorf36 but not TET2-T1884A 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). 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, TET2-Y1649stop, 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 single-cell assay agree on TET2-E1357stop 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. The subclonal lineage inferred by Jan et  al.'s single-cell assay is shown in Figure 5A and only contains two lineages, one with only TET2-E1357stop and the other with the other five SNVs. The TET2-E1357stop lineage genotype has probability 0.81 in our posterior, however the second genotype has a relatively small probability (0.06) under the posterior although we note that the genotype that contains all SNVs but ZMYM3 has a posterior probability of 0.32 (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 single-cell 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 semi-manual procedure in which they first automatically grouped SNVs into subclonal lineages using k-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 complete-data 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 well-correlated. 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 semi-manual 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 expert-driven, semi-manual 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 lineage-defining 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 well-mixed 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].

Dirichlet process mixture models
Consider the problem of clustering N objects {x i } N i=1 using a Bayesian finite mixture model of K components (clusters) with the following generative process [27]: where ω are the mixing weights such that K k=1 ω k = 1, α is the concentration parameter of the symmetric Dirichlet prior placed on the mixing weights, z i ∈ {1, . . . , K} is the cluster assignment variable, H is the prior distribution from which the component parameters {φ k } are drawn, F (φ) 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: where δ φ is a point mass centered at φ and G ∼ DP(α, H), i.e., G is a draw from a DP with base distribution H and concentration parameter α. The stick-breaking process can be viewed as recursively breaking sticks of length k−1 =1 (1 − β ), starting with a stick of unit length. The beta variates {β k } 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 GEM(α) denote the stick-breaking process over ω. Replacing the Dirichlet prior in the finite mixture model (1) with the stick-breaking process prior results in the following generative process for infinite mixture models: An alternative view of the above generative process produces component parameters {φ i } by drawing samples from G resulting in the following generative process: Note that in the above process every object {x i } N i=1 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 {φ i } N i=1 will take on the same value from the set {φ k } K k=1 of unique parameters.

Tree-structured stick-breaking process
The stick-breaking 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 stick-breaking 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 = ( 1 , . . . , p ) denote a sequence of positive integers used to index the nodes of the tree. Let = κ denote the zero-length string, i.e., the root of the tree. Let | | indicate the length of the sequence and therefore the depth of node . Let i denote the sequence formed by appending i to . The children of node is the set { i : i ∈ 1, 2, . . .} and let the ancestors of be denoted by the set { : ≺ }. The interleaved, two-layered stick-breaking construction is as follows: The ν and (1 − ν ) 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 (α(·) : N → R + ) and is defined to be α(j) = λ j α 0 with α 0 > 0 and λ ∈ (0, 1] [19].

PhyloSub model
We follow Shah et al. [10,30] . .} denote the probability of sampling a reference allele from the reference population. This value depends on the error rate of the sequencer. Let µ v i denote a vector whose entries, . .}, are the probabilities of sampling a reference allele from the variant population with genotype g at position i. Let π i denote the vector whose entries, π g i ∈ {π B i , π AB i , π BB i , . . .}, are the probabilities of the variant population at position i to have the genotype g. Let δ i denote the pseudo-count parameters of the Dirichlet prior over π i . Let G i denote the genotype of the variant population at position i. Letφ i denote the fraction of cells from the variant population, i.e., the SNV population frequency at position i, and 1 −φ i denote the fraction of cells from the reference population at position i. The observation model for allelic counts has the following generative process [10]: The posterior distribution ofφ i 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 where Γ(·) is the Gamma function. The Dirichlet process prior DP(α, H) 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 tree-structured stick-breaking process prior (4) described in the previous section. The probabilistic graphical model for allelic counts with the tree-structured stick-breaking 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 {φ i } are indicated in unshaded nodes. The prior/base distribution H of the SNV population frequencies is the uniform distribution Uniform(0, 1) for the root node and Uniform(0, φ par(v) − w∈S(v) φ w ) for any other node v in the tree, where par(v) denotes the parent node of v and S(v) is the set of siblings of v. 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 tree-structured stick-breaking process instead of a Dirichlet process (cf. (5)) to generate the set of SNV population frequencies {φ i }. Given this model and a set of , the tree structure and the SNV population frequencies {φ i } 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 {z i }, sampling stick lengths ν and ψ i , sampling stick-breaking hyperparameters α 0 , γ and λ, and sampling the SNV population frequencies {φ i }. 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 φ v of every non-leaf node v in the tree is greater than or Figure 10: PhyloSub graphical model for single sample. Probabilistic graphical model for allelic counts with tree-structured stick-breaking process prior. Observed variables and hyperparameters (inputs to the model) are indicated in shaded nodes.
equal to the sum of the SNV population frequencies of its children. To enforce this constraint, we introduce a set of auxiliary weights {η v }, one for each node, that satisfy v η v = 1, and rewrite the observation model for allelic counts (5) explicitly in terms of these weights resulting in the following posterior distribution: where we have used {η i } 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 Uniform(0, η par(v) ) for any other node v in the tree, where par(v) denotes the parent node of v. When a new node w is added to the tree, we sample η w ∼ Uniform(0, η par(w) ) and update η par(w) ← η par(w) − η w . This ensures that v η v = 1. This change is crucial as it allows us to design a Markov chain that converges to the stationary distribution of {η v }. The SNV population frequency for any node v can then be computed via where D(v) and C(v) are the sets of all descendants and children of node v 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 (η v , φ v ) for every node v by traversing the tree in a breadth-first fashion. The input to this algorithm is the current state of the tree T = (V, E) where V is the set of vertices and E is the set of edges, and the output is a multi-dimensional sample of SNV population frequencies φ = {φ 1 , φ 2 , . . . , φ |V | } η ∼ Dirichlet(ση (t−1) + 1)

8:
if log(r) < a then 9: η (t) ← η 10: We use Metropolis-Hastings algorithm to sample from the posterior distribution of the auxiliary weights {η i } (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 {η i }. The inputs to the sampling algorithm are the current state of the tree T = (V, E), a scaling factor σ, and the number of iterations T . The output is a sample from the posterior distribution of η = {η 1 , η 2 , . . . , η |V | } and its corresponding φ = {φ 1 , φ 2 , . . . , φ |V | }.

Extension to multiple tumor samples
PhyloSub (cf. Figure 10) can be easily extended for multiple tumor samples. We allow the treestructured stick-breaking process prior (4) to be shared across all the samples. Let a t i and b t i denote the number of reads matching the reference and the variant allele respectively at position i for sample t ∈ {1, . . . , S}, and let d t i = a t i + b t i . Letφ t i denote the fraction of cells from the variant population, i.e., the SNV population frequency at position i for sample t, andη t i 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 Metropolis-Hastings move based on the distribution S t=1 p(η t | ·), where {η 1 , η 2 , . . . , η S } 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 P → Q in the tree is proportional to the fraction of MCMC samples in which SNV P first appears in a subclonal lineage that was the parent of the subclonal lineage that Q 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 N × N co-clustering matrix C, whose elements C ij is the difference between the number of samples in which i and j were assigned to the same SNV cluster and the number of samples in which i and j were assigned to different SNV clusters. The algorithm estimates a symmetric N × N cluster indicator matrix Y , whose elements Y ij = 1 if i and j are assigned to the same cluster and Y ij = 0 otherwise. This cluster indicator matrix Y 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 burn-in of 100 samples. We also fix the number of iterations in the Metropolis-Hastings algorithm to 5,000 and set the scaling factor for the Dirichlet proposal distribution to σ = 100. We run the MCMC samplers multiple times with different random initializations and pick a single run based on the complete-data likelihood trace and its auto-correlation 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 complete-data log likelihood traces and the corresponding autocorrelation function plots after the burn-in 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 {(a i , d i , µ r i , µ v i , δ i )} N i=1 , are provided in the Additional file 2 (Tables S5 -S10).