TrAp: a Tree Approach for Fingerprinting Subclonal Tumor Composition
Francesco Strino, Fabio Parisi, Mariann Micsinai, Yuval Kluger
1 Department of Pathology, Yale University School of Medicine, New Haven, Connecticut, United States of America
Email: yuval.kluger@yale.edu
These authors contributed equally to this work.
Abstract
Revealing the clonal composition of a single tumor is essential for identifying cell subpopulations with metastatic potential in primary tumors or with resistance to therapies in metastatic tumors. Sequencing technologies provide an overview of an aggregate of numerous cells, rather than subclonalspecific quantification of aberrations such as single nucleotide variants (SNVs). Computational approaches to demix a single collective signal from the mixed cell population of a tumor sample into its individual components are currently not available.
Herein we propose a framework for deconvolving data from a single genomewide experiment to infer the composition, abundance and evolutionary paths of the underlying cell subpopulations of a tumor. The method is based on the plausible biological assumption that tumor progression is an evolutionary process where each individual aberration event stems from a unique subclone and is present in all its descendants subclones. We have developed an efficient algorithm (TrAp) for solving this mixture problem. In silico analyses show that TrAp correctly deconvolves mixed subpopulations when the number of subpopulations and the measurement errors are moderate. We demonstrate the applicability of the method using tumor karyotypes and somatic hypermutation datasets. We applied TrAp to SNV frequency profile from ExomeSeq experiment of a renal cell carcinoma tumor sample and compared the mutational profile of the inferred subpopulations to the mutational profiles of twenty single cells of the same tumor. Despite the large experimental noise, specific cooccurring mutations found in clones inferred by TrAp are also present in some of these single cells. Finally, we deconvolve ExomeSeq data from three distinct metastases from different body compartments of one melanoma patient and exhibit the evolutionary relationships of their subpopulations.
Introduction
The mechanisms of cancer evolution and metastatic onset are still largely unknown. The diversity, complexity and evasive nature of tumor biology are central reasons for the seemingly slow progress in the cure of most cancer types, particularly in controlling the ability of tumor populations to spread. Tumor populations are dynamic aggregates of constantly evolving subclones, each carrying a variety of genomic aberrations [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. This genetic heterogeneity is often associated with differences in the biological behavior of different cell subpopulations. Some of these subclones are likely to be the primary instigators of invasion, metastases or relapse following treatment [11, 17, 18, 19, 20]. An effective characterization of the aggressive potential of tumors at early stages has an enormous potential to guide new clinical interventions and translational research [21].
Recently, several efforts have been made to provide a complete genealogical perspective of cancer evolution [22, 23, 24, 25, 26]. Using fluorescent labeling techniques, or more recently, single cell sequencing, it is technically possible to separate single cells from tumor samples in order to investigate their evolutionary patterns [27, 22, 23, 24, 25, 26, 28, 29, 30, 31]. However, these approaches are limited to either a small number of fluorescent markers [23, 32], or to a relatively small number of single cells. On one hand, the identification and selection of uncharacterized subclones in highthroughput experiments is beyond the capabilities of current cellsorting technologies; on the other hand, isolation and profiling of enough single cells in order to obtain a statistically representative sample of a tumor composed of millions of cells has, currently, prohibitive costs. For this reason, genomics profiling of tumors still relies on pooling in order to provide global averaged signals over the subclonal population within a tumor sample [33, 34, 35, 36]. Computational methods for identifying subclones, quantifying their relative abundance and monitoring their emergence and dynamics could prove extremely useful for the analysis of the heterogeneity of these pooled samples. This problem has been often overlooked due to its mathematical complexity.
We herein present a mathematical approach to demix signals from heterogeneous cell populations into their subclonal components and subsequently unveil the underlying dynamic tumor heterogeneity. Our proposed method is related to the problem of blind source separation [37, 38, 39, 40, 41, 42, 43, 44], where both the underlying sources and their relative composition are unknown. In contrast to blind source separation methods we cannot assume that the underlying sources are statistically independent, we have no prior knowledge of the number of sources and we have at our disposal only one mixture of the unknown sources. This mathematical problem has a vast number of solutions and can be addressed only if additional constraints are imposed. Hence, we assume that tumor cell populations develop in a parsimonious evolutionary process. Furthermore, based on empirical observations, we introduce a sparsity constraint that limits the number of subpopulations. Distinctively from the standard problem of phylogeny [45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56], where each species is observed and measured separately, our methodology, which we term Tree Approach to Clonality (TrAp), is specifically designed to deconvolve a single aggregate signal into its different subclonal components. TrAp incorporates biologically motivated constraints which allow us to infer: 1) the composition of the subclones in a single aggregate sample 2) the abundance of each subclone and 3) the evolutionary path of the subclones.
The paper is organized as follows: we first define the subclonal deconvolution problem and present an efficient algorithm for finding all its solutions. Using in silico simulated data we verify that the algorithm is able to correctly deconvolve mixed subclonal populations and that the method is robust to realistic measurement errors. Further, the solution is often unique when the number of populated subclones is moderate. In addition, we also show that TrAp can correctly deconvolve random mixtures of karyotypes of several cells from the same tumor biopsy or from mixture of sequences generated in a study involving somatic hypermutations in B cells. We subsequently compare the mutation profiles of twenty ExomeSeq single cell experiments to those inferred using an aggregate signal generated by exome sequencing from the same renal cell carcinoma tumor. Finally, we apply TrAp to ExomeSeq data from three metastases from three distinct body compartments and compare their subclonal compositions and evolutionary histories.
Results
The results are divided in four parts. In the first part we describe our novel TrAp algorithm for subclonal deconvolution of aggregate genomic signals consisting of aberrations’ frequencies and we show that the TrAp algorithm always identifies at least one solution; we also introduce an extension to our work that is suited for signals where each locus can be affected by distinct consecutive aberrations (e.g. several consecutive point mutations affecting the same nucleotide). In the second part, we simulate noisy aggregate signals constructed by random linear combinations of simulated subclonal aberration profiles. We used these simulated data to show that TrAp can correctly deconvolve mixtures of evolutionarily related subclones even in presence of levels of noise that are typically found in current genomics experiments. In the third part, we generated realistic aggregate signals by mixing subclonal genomic profiles obtained from cytogenetic analyses using random coefficients. We generate these data separately for each patient and show that, for nearly all aggregate samples, TrAp recovered the subclonal components. Similarly, we generated realistic aggregate signals from somatic hypermutated (SHM) regions from B cells. As we show, somatic hypermutation is a particularly suitable system for the framework of our TrAp algorithm, which successfully recovered all components from the aggregate signals. In the fourth part, we apply our approach to exomesequencing experiments. We present an analysis of recent singlecell exome sequencing from a renal cell carcinoma study where, besides a collection of twenty single cells, an aggregate has also been measured. Despite the reported difference between the aggregate and mean aberration profile of the single cell experiments, TrAp could still identify subclones with cooccurring aberrations consistent with cooccurring aberrations found in direct single cell measurements. Finally, we analyze three metastases from separate body compartments of a melanoma patient and compare their inferred evolutionary patterns in a genomic region surrounding the DCC gene.
The TrAp algorithm
The subclonal deconvolution problem
We consider a population of cells where each cell can be described by a binary vector, which we call genotype. Each element of the genotype vector has an aberration state modeled as a binary variable, e.g. the presence/absence of a mutated nucleotide in a specific genomic position or the presence/absence of a specific copy number variation event in a specific locus. For each cell, the th element of the genotype vector is if the th aberration is present in the cell and if the aberration is absent. A subclone is a collection of all cells that have identical genomewide aberration profile. A subclone is populated in the sample if the fraction of cells sharing the subclone’s genome is greater than zero and can be detected by the experiment.
We define the subclonal deconvolution problem as the task of demixing an aggregate measurement into a linear combination of (unknown) subclonal genotypes. The only information that is required as input is the aggregate frequency vector , whose elements correspond to the frequency of the th aberration in the sample cell population. For efficiency, we remove from the genome all genotype entries that are homogenous within the population (i.e. or ), as they do not need to be deconvolved. Next, to ensure that the aberrationfree non cancerous cells (wildtype) are included in the solution of the problem, we add one dummy aberration to all the normal and cancerous cells in the sample. By construction, the aggregate frequency of this dummy aberration is equal to . Finally, without loss of generality, we sort the aggregate frequency vector in descending order such that , where is equal to the number of aberration events considered including the dummy aberration . The subclonal deconvolution problem can be written in matrix notation as
(1) 
where is a matrix of size whose elements are if aberration is present in subclone and otherwise; is the size of the vector ; is the number of subclones; and is a vector of size where each element corresponds to the frequency of subclone in the sample. We note that, without introducing the wildtype aberration, the wildtype subclone would correspond to a vector of zeros and we would not be able to infer the frequency of the wildtype component using the linear model of Equation (1). Furthermore, since the dummy aberration is present in the wildtype and all other subclones, it follows that (i) and (ii) . We note that since , and are all unknown in this problem there is an intractable number of possible solutions. As previously discussed [57], for the system is underdetermined and the aggregate signal cannot be uniquely deconvolved by solving the linear system and it is not even possible to find parsimonious unique solutions using sparse reconstruction methods. However, by introducing additional biologically motivated constraints to the model, we can dramatically reduce the number of possible solutions, such that the problem may have a tractable number of optimal solutions, ideally only one. We therefore seek the family of solutions (TrApsolutions) that sequentially satisfy the following constraints:

Evolutionarity. The subclones are generated in an evolutionary process starting from a subclone with no aberrations. Every subclone is generated from an existing subclone by adding to it a single new aberration event.

Parsimony. The number of subclones that are generated during the evolution process is minimal.

Sparsity. The number of populated subclones (i.e. the number of subclones for which ) is minimal.

Shallowness. The depth of the evolutionary tree (i.e. the number of generations) is minimal.
A schema of a TrAp solution is shown in Figure 1.
The evolutionarity constraint is used in many biological systems, in particular when studying development of cell populations [45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55]. The parsimony constraint is typically satisfied because the expected probability of a specific aberration event in a nucleotide or a locus is low and it is therefore unlikely that an event occurs more than once and independently in distant subclones. This constraint is the main criterion used to determine the optimality of maximum parsimony methods commonly used in phylogeny [58, 59, 48, 56]. The parsimony constraint dramatically reduces the number of possible solutions of Equation (1) because it limits the number of subclones . The sparsity constraint is justified by the fact that some subclones may die out or may be too rare to be detected. Also, it has been shown in several studies that few subclones acquire an evolutionary advantage and outgrow the other subclones [60, 61, 5, 13], thus reducing the number of populated subclones. The shallowness constraint is justified as excessive genomic instability may not be viable, thus evolutionary trees tend to be shallow and wide rather than deep and narrow.
Properties of TrAp solutions
In this subsection we show that in the subclonal deconvolution problem the evolutionarity and parsimony constraint can always be satisfied by a naïve model in which each aberration event occurs exactly once during evolution (i.e. ). We call any solution with an Nsolution and its evolutionary tree an tree.
The solution optimally satisfies the evolutionary and parsimony constraints. Since all detected aberrations need to occur at least once in the evolutionary tree, the number of subclones must be greater than or equal to the the number of aberration events considered, i.e. . We note that it is always possible to construct a valid solution for Equation (1) of exactly subclones for every aggregate frequency vector y. Specifically, for a cascadelike evolutionary process with no branching, where the wildtype subclone is the root of the tree and every other subclone is a direct descendant of the subclone , the solution of Equation (1) is given by:
(2) 
While this cascadelike solution satisfies both evolutionarity and parsimony constraints, this solution is not necessarily optimal with respect to the sparsity constraint and is the least optimal in terms of the shallowness constraint. Furthermore, the existence of this solution guarantees that there is always at least one solution to the subclonal deconvolution problem.
Since we imposed that the four constraints to the subclonal deconvolution problem must be satisfied sequentially, any solution with subclones will always be less optimal than the solution given in Equation (2), regardless of its sparseness and shallowness. We can thus limit the search space of the TrAp algorithm to solutions. For this subset of solutions, the vector is of size and is a square matrix of size . Importantly, the index of each subclone indicates the subclone for which the th aberration occurs for the first time. Hence, for any TrAp solution there is a onetoone correspondence between the th aberration and the subclone whose index indicates that it evolved from its parent subclone by acquiring the th aberration. Moreover, each aberration event of an tree occurs exactly once and cannot be reverted during evolution. Each aberration is thus present only in subclone and its subclonal descendants:
(3) 
where is the ancestor indicator variable that is equal to if subclone is an ancestor of subclone and otherwise. We note that , which means that the wildtype clone is always the root of the evolutionary tree, as required by the evolutionarity constraint.
In Equation (3), we expressed the aggregate frequency as the sum of the frequencies of all subclones descending from . We now express the aggregate frequency as a function of the direct descendants of . We define the parent indicator variable , which is if is the parent of (i.e. if subclone is the result of a single aberration event in subclone ) and otherwise. Finally, using the parent indicator variables we express in terms of the aggregate frequencies of the direct descendants of
(4) 
Equation (4) can be rearranged to express the vector in terms of and the parent indicator matrix :
(5) 
where is the matrix whose elements are given by the parent indicator variables . Because of the evolutionary constraint, the matrix has only nonzero elements, reflecting the fact that each subclone except the wildtype has exactly one parent subclone, i.e. and . Furthermore, an important corollary of Equation (4) is that the subclone is not populated if and only if
(6) 
In other words, the clone is not populated when the aggregate frequency of aberration is equal to the sum of the aggregate frequencies of all the direct descendants of the subclone . Therefore, the number of nonpopulated subclones of the tree encoded by is given by the number of aberrations that satisfy Equation (6).
Finally, we summarize the relationships between and the indicator variables and . Using Equation (3), we can express the subclonal deconvolution problem as , where is the identity matrix and is the matrix of whose elements are given by the ancestor indicator variable . Furthermore, Equation (5) allows us to write the subclonal deconvolution problem as . We can therefore express the relationships between the matrices , and as:
(7) 
A bruteforce algorithm for solving the subclonal deconvolution problem
We now observe that the relationship between two subclones and such that (which also implies as the vector is sorted), must be one of the following: i) is an ancestor of , i.e. , and ; or ii) and are on separate branches, i.e. , and . This property implies that all evolutionary trees can be generated iteratively by starting with the wildtype clone and adding the clone at step to all trees generated at step . In detail, for any tree that can be generated using subclones we generate a new tree by adding the subclone as direct descendant of subclone for all for which the resulting (calculated using Equation (5)) remains nonnegative after adding .
For completeness, when the subclones and can be either on separate branches or on the same branch. If they are on the same branch, then there is an ambiguity regarding the order in which the two aberrations occur (Figure 2). However, in the case that these two subclones are on the same branch, the aberration profile of the ancestor subclone (shown in green in Figure 2) is not informative because this subclone is not populated () and aberrations and are both present in the descendant subclone regardless of the order in which they occur. Since these two aberrations cannot be observed separately (i.e. the coefficient associated with the ancestor aberration is zero, whereas the associated with both aberrations could be nonzero), we only output the solution for which is an ancestor of (left solution in Figure 2). This choice ensures that for every pair of aberrations , cannot be an ancestor of . A stepbystep example of the TrAp solution obtained using the bruteforce algorithm is given in Figure 3.
The upper bound on the number of possible evolutionary trees is thus , as every subclone can only be the direct descendant of subclones. This bound is significantly smaller than the number of all possible trees with labeled vertices, which is according to Cailey’s formula [62, 63]. Furthermore, we note that the parent indicator matrix and are upper triangular and that both and are of rank and invertible, which guarantees that Equation (7) can always be used to switch between the representation with the parent indicator variable (Equation (5)) and the representation with the subclone matrix (Equation (1)). We also note that, given a matrix (or ), the vector is uniquely determined. In particular, if the matrix is known, the vector can be efficiently found by solving the linear system using backsubstitution (i.e. by solving Equation (3) first for , then using to solve for and repeating through ).
TrAp: a fast algorithm for solving the subclonal deconvolution problem
As we have shown, the number of nonpopulated subclones of a given tree is equal to the number of aberrations that satisfy Equation (6). Moreover, in order to satisfy the sparsity constraint of a solution, we do not need to know the topology of the whole evolutionary tree, but only the subset of rows of the matrix that satisfy Equation (6). We now leverage on this property to efficiently generate sparse solutions to the subclonal deconvolution problem.
First, we group each subset of subclones that satisfy Equation (6) into a first generation tree , which we define as the subset of subclones such that the subclone is not populated (i.e. ) and that subclones are its direct descendants. Each first generation tree is represented by a row of the matrix. For example, there are three first generation trees for the aggregate signal , namely , and (Figure 4). We note that the optimal TrAp solution for this example contains the first generation trees and (Figure 1). Furthermore, a matrix associated with a first generation tree must follow the evolutionary constraints previously described () and thus the first generation tree also gives complete information about the columns of corresponding to the direct descendant subclones . For example, the first generation tree implies that and for every (Figure 4).
Next, we define a partial tree as a collection of first generation trees that can jointly be contained in a full evolutionary tree. Since each first generation tree can be represented by a row of the matrix, a partial tree that is comprised of first generation trees can be represented by rows of the matrix. In the example above, the partial tree that contains the first generation trees and is represented by rows and of the matrix in the bottom panel of Figure 4. Similarly to first generation trees, the matrix associated with a partial tree must follow the evolutionary constraint, which implies that not all combinations of first generation trees give rise to partial trees. In the example above, the partial trees and can be combined to generate the partial tree (Figure 4 bottom), whereas the pairs and cannot be combined to generate partial trees. Therefore, in the example above the possible partial trees are , the empty partial tree and the partial trees , and . Moreover, we note that all trees that contain a partial tree comprising of first generation trees have populated subclones. This implies that TrAp solutions, which must satisfy the sparsity constraint, must also contain one of the partial trees comprising the maximum number of first generation trees. In the example above, the optimal TrAp solution (Figure 1) contains the partial tree , which is the only partial tree comprising two first generation trees.
Since all TrAp solutions contain the maximum number of first generation trees, the TrAp algorithm dramatically reduces the search space by identifying the optimal partial trees and later utilizing them as starting points for rapidly building all the sparsest solutions to the subclonal deconvolution problem. In details, the TrAp algorithm solves the subclonal deconvolution problem in the following steps (Figure 5):

Identify all first generation trees from the aggregate signal vector .

Combine all first generation trees to generate all partial trees.

Discard partial trees that do not have the minimum number of populated subclones.

Generate all evolutionary trees consistent with the partial trees comprising the maximum number of first generation trees. This step is done iteratively for each partial tree, similarly to the way described for the bruteforce algorithm. The only difference is that, when the parent clone of a first generation tree is added to the tree, the subclones are automatically added as direct descendants of .

Optimize the shallowness constraint by sorting the generated solutions by the depth of the generated tree.
The performance of the TrAp algorithm is equivalent to the bruteforce approach when there are no first generation trees (i.e. when all subclones are populated) and becomes superior to a bruteforce approach when . Furthermore, the algorithm can be modified by imposing additional constraints that need to be satisfied at each step of the iterative tree reconstruction procedure (an example of such modification is presented in the subsection on extension of TrAp to nonbinary aberrations). The TrAp algorithm by default returns only the solution(s) that optimize all the constraints. In addition, the user can specify parameters to relax the sparsity and shallowness constraints.
Generalization of TrAp to nonbinary aberrations
In the previous sections we represented the genome of each subclone by a vector of binary values whose entries represent if a genomic position is in a normal state (0) or in an aberrated state (1). In general the number of states in a given genomic position could be larger than two and hence subclones cannot be represented by vectors of binary values without loss of information. For example, a nucleotide found in the reference genome or in the germline at a specific position may undergo multiple distinct point mutation events into more than one specific nucleotide. In this subsection, we describe an extension of the TrAp algorithm to deal with such cases.
We consider the generalized subclonal deconvolution problem in which the genome consists of positions each of which can have different states. We also assume that the genome of the wildtype subclone is known. The only information that we utilize as input is the aggregate frequency matrix whose elements correspond to the observed frequency of the aberrated state at position . We note that, by construction, and that . To utilize our framework for solving the subclonal deconvolution problem, we convert the information encoded in as a vector whose elements represent frequencies of binary events. We perform this transformation in several steps. First, we vectorize the matrix by concatenating its rows to construct the vector , which has elements. Then, we remove from the vector the entries for which as they are not informative. As a result, for each position there are entries, where ranges from (only the unmutated state is observed) to (all aberrated states are observed). For illustration of these first steps, we consider a toy example of a genome of length three whose wildtype sequence is ”CAT” and we analyze an aggregate sample made of three subclones with sequences ”TCT”, ”TAC” and ”CGT” mixed with frequencies , and , respectively (Figure 6). In this example the vector consists of the elements , , , , , and (Figure 6).
Next, we wish to design a binarization matrix to encode the information contained in as a vector whose elements represent the frequency of binary aberrations and can thus be used as input to the subclonal deconvolution problem for the whole genome (Equation (1)). For every position , we assume that each state () is reached by a sequence of aberration events. We denote by an aberration event to state at position . In the example above, there are two states at position . The unmutated state is reached by a dummy aberration (in analogy to the dummy aberration of the wildtype clone in the subclonal deconvolution problem) and the aberrated state is reached by the sequence of the dummy aberration followed by the aberration . Since the dummy aberration is present when we observe states or at position , the frequency of the unmutated aberration is . However, the aberration is present only when we observe the state , therefore the frequency of the aberration is (Figure 6). Since measurements at different genomic positions do not affect one another, we construct as a blockdiagonal matrix:
(8) 
where and . Combining Equation (1) and Equation (8) gives
(9) 
It is important to note that for any pair of different states and at position where and , the value of defines the ancestral relationship between the aberration events and . Using the ancestor indicator variable we can express this relationship as . To preserve these ancestral relationships in both sides of Equation (9) (we recall that ), the element of must be equal to the element of for every pair of states and at position , where , and .
In order to find the solutions to Equation (8), we solve independently each subproblem and we require that the matrix encodes a tree which satisfies the evolutionary and parsimony constraints and whose root is the dummy unmutated aberration at position . The number of solutions for each subproblem is equal to the number of trees with labeled vertices and, as discussed earlier, this number is equal to . We then denote by the th solution () and by the aggregate frequency vector associated to it. The solutions to the problem for are:
 .

There is only one solution. Since the input consists only of the unmutated state, it follows that . The solution is and the corresponding aggregate frequency vector is .
 .

There is only one solution. Assuming the input is and the unmutated state is , the solution is and the corresponding aggregate frequency vector is .
 .

There are solutions. Assuming the input is and unmutated state is , the solutions are , and and their corresponding vectors are given by , and . The values and define the ancestral relationship between aberrations and . In the first solution and must be on separate branches, in the second solution must be an ancestor of , and in the third solution must be an ancestor of .
The number of solutions grows dramatically with ( solutions for , solutions for , solutions for ). Therefore, herein we explicitly show the solutions for and implement the method to adress practical scenarios, such as nucleotide point mutations ().
The set of all possible solutions to Equation (8) is given by the Cartesian product of the solution sets of each subproblem. This implies that the number of solutions of Equation (8) is . We denote by the block matrix solution associated with the blocks representing the solutions of each subproblem. Similarly, we denote by the aggregate signal vector given by . It is possible to reduce the number of rows and columns in Equation (9) by substituting all the rows corresponding to the unmutated states (shown in green in Figure 6) with a single dummy wildtype aberration (shown in blue in Figure 6). In our toy example, one aberrated state is observed at positions and (), while two aberrated states ( and ) are observed at position (). Therefore, there are three solutions to Equation (8): , and . Figure 7 illustrates the best solution for each of these binarization matrices in a graphical and matrix form. The solution associated with is the only TrApsolution of the generalized subclonal deconvolution problem since it has the minimum number of populated subclones (sparsity constraint).
In summary the generalized subclonal deconvolution problem can be solved as follows:

Vectorize the aggregate frequency matrix and identify all binarization matrices (Equation (8)) compatible with the vector .

For each binarization matrix , identify all first generation trees from the aggregate signal vector and combine the first generation trees to generate all partial trees compatible with .

Discard partial trees that do not have the minimum number of populated subclones.

Generate all evolutionary trees consistent with the partial trees comprising a maximum number of first generation trees. This step is performed as described above, but with the additional constraint that for any pair of states and at position , where , and .

Optimize the shallowness constraint by sorting the generated solutions by the depth of the generated tree.
Benchmarks using simulated data
The models presented above show that TrAp is an efficient algorithm for inferring subclonal components from the aggregate measure. In particular we have shown that in the absence of noise TrAp returns the exact solution when the underlying subclonal population satisfies reasonable constraints and that the algorithm is always able to find at least one solution. However, experimental measurements are often noisy and can only have finite precision. In this section, we first discuss two approaches to treat noisy input and then we apply our algorithm to simulated aggregate signals from random in silico trees showing that TrAp is robust to typical noise levels found in genomic experiments.
Handling Measurement Errors
In order to accommodate different types of noise that may arise in genomic data, we incorporated two error models in the TrAp algorithm. In both error models, the input to the TrAp algorithm requires an additional vector of size whose elements are related to the precision of the measure . The error related to the dummy variable is denoted by and is set to as is a constraint of the model and thus must vanish.
First, we examine the bound error model in which we assign a threshold to the error of every underlying aggregate signal such that each measured signal will be in the range . Equation (6) is then modified accordingly and we can state that the subclone defined by aberration is not populated if and only if:
(10) 
Next, we implement a normal error model assuming normally distributed measurements errors using a confidence interval approach to determine whether a subclone is populated or not. Specifically, we assume that the underlying aggregate signal is normally distributed around the observed signal, i.e. . We substitute each term of the lefthand side of Equation (6) by its normal distribution in order to derive the distribution with mean and variance . Using the distribution of and a desired confidence level (default ), we can define that clone is not populated if falls within the confidence interval , where is the quantile of the distribution of .
Once the error model is chosen, the algorithm generates all optimal TrAptrees in a similar fashion to the noisefree case. The main difference is that in the first step of the TrAp algorithm, Equation (10) (or a confidence interval on the distribution ) is used instead of Equation (6) to identify first generation trees. Moreover, instead of using backsubstitution for finding the vector , we solve the nonnegative linear least square problem with the constraint for all nonpopulated subclones associated with the parents of the first generation trees. This fitting allows us to obtain a value of exactly zero for all nonpopulated subclones and to distribute measurement error more evenly in the vector .
Deconvolution of simulated noisy aggregates
To confirm that TrAp can correctly infer the subclonal composition from aggregate noisy signals with typical noise levels found in genomic experiments, we performed simulations starting with random in silico evolutionary trees (see Methods) with different numbers of aberrations and different numbers of populated subclones . For each tree, we also studied the effect of different magnitudes of measurement errors and we investigated the operating conditions for which TrAp would assign the best score to the true solution.
For each aggregate signal from a random tree we examined: (i) whether the true solution (i.e. the solution associated with the simulated tree), which is by construction one of the possible solutions to the subclonal deconvolution problem, had the minimum number of populated subclones among all solutions (sparsity constraint), (ii) whether the true solution had the minimum number of populated subclones and minimum number of levels of the evolutionary tree among all solutions generated by TrAp (sparsity and shallowness constraints), and (iii) whether the true solution was the only TrAp solution (sparsity and shallowness constraints and uniqueness of the optimal solution).
The results of the simulations show that aggregate signals from sparse trees are deconvolved correctly even in presence of typical noise levels of sequencing experiments (Figure 8). We note that for simulations of nonsparse trees, TrAp generates a large number of possible solutions of which only one is the true solution. Furthermore, in the presence of high levels of noise, the TrAp algorithm identifies a large number of first generation trees that satisfy Equation (10) and generates solution trees whose number of populated subclones is smaller than the number of populated subclones of the true solution.
Analysis of simulated mixtures of biological data
Deconvolution of a mathematical mixtures of karyotyping data from single tumor biopsies
After showing that our approach can correctly deconvolve aggregate signals of subclones with a treelike genealogy, we sought to investigate whether actual subclonal populations can be charted on evolutionary trees. For this purpose we analyzed the Mitelman database, consisting of cytogenetic analyses of more than biopsies (see Methods). For each tumor type we counted how frequently the relationships between cancer clones from the same biopsy could be explained by an evolutionary tree that follows the evolutionarity and parsimony constraints (but not necessarily the sparsity and shallowness constraints). We found that almost all biopsies in the Mitelman database can be represented by evolutionary trees (Table 1), with the notable exception of astrocytoma of grades III and IV (Figure 9).
1  2  3  4  4  

1  100% (19078)  n/a  n/a  n/a  n/a 
2  100% (5150)  100% (923)  0% (3)  n/a  n/a 
3  100% (1830)  100% (367)  94% (83)  0% (2)  n/a 
4  100% (991)  100% (182)  89% (27)  89% (18)  n/a 
5  100% (656)  100% (120)  88% (33)  100% (8)  100% (5) 
6  100% (445)  100% (66)  92% (13)  100% (6)  50% (4) 
7  100% (333)  100% (58)  89% (9)  25% (4)  100% (2) 
8  100% (241)  100% (37)  86% (7)  100% (3)  50% (2) 
9  100% (228)  100% (26)  60% (5)  0% (1)  100% (1) 
10  100% (174)  100% (14)  100% (2)  n/a  50% (2) 
11  100% (196)  100% (25)  67% (3)  67% (3)  n/a 
12  100% (156)  100% (16)  100% (3)  0% (1)  50% (2) 
13  100% (137)  100% (21)  50% (2)  n/a  100% (1) 
14  100% (94)  100% (12)  n/a  100% (1)  100% (1) 
14  100% (152)  100% (22)  57% (7)  25% (4)  25% (4) 
We mathematically mix all karyotypes of each single patient from the Mitelman database and apply the TrAp algorithm for each of these mixtures. The ability of the TrAp algorithm to extract the correct underlying clonal or subclonal components depends on the number of actual components (columns) and the multiplicity of aberrations studied in each mixture (rows). The frequency in which TrAp is able to recover the correct underlying components is shown in percentages. The number of mixtures for a given size of aberration multiplex (row) and given number of actual underlying components (column) is shown in parentheses. Note that when the column index is greater than the row index (entries shown in bold) the parsimony constraint cannot be satisfied.
Next, we investigated whether the TrAp algorithm could uniquely deconvolve mixtures of the cancer subclones within a biopsy. As these aggregate signals are obtained by mixing actual subclonal profiles, we consider these signals to be more realistic than our previous in silico simulations. For each biopsy, we generated in silico mixtures by combining the cytogenetic profiles of each subclone using random nonnegative coefficients (see Methods). We then applied our TrAp method to deconvolve in silico mixtures of these cancer clones and found that TrAp could correctly deconvolve of these realistic aggregate signals. However, we note that this task was trivial for a significant fraction of the biopsies whose number of aberrations and/or subclones is small. Moreover, the TrAp algorithm inferred also intermediate nodes in the evolutionary tree that did not correspond to any of the cytogenetic profiles for the biopsy, providing a plausible picture of the evolutionary order in which the aberrations occurred. Figure 10 shows the result of two deconvolution simulations, one from a melanoma sample with subclonal populations [64] and one from an adenocarcinoma sample with subclonal populations [65]. An interesting, yet rare (only examples, shown in bold in Table 1), case of biopsies showed more clones than aberrations. Albeit a treelike genealogical relationship can be constructed, these biopsies do not satisfy the parsimony constraint since the number of subclones is greater than the number of mutations . For this reasons, their genealogy cannot be reconstructed by the TrAp algorithm or by any other method that makes use of a similar parsimony constraint [58, 59, 48, 56].
Deconvolution of a mathematical mixture of somatic hypermutation data with polyallelic mutations in a single nucleotide
Somatic hypermutation (SHM) introduces mutations that target the variable regions associated with immune adaptivity, i.e. Ig loci. In particular, SHM involves a programmed process of mutations that affects the variable regions of immunoglobulin genes and starts from an initial dividing single cell (a naïve B cell in this case). All descendants of the founder cell accumulate mutations and, at the same time, are subjected to a strong selective pressure. For this reason, SHM is a particularly good biological model system to test our deconvolution method, which imposes treelike evolutionary constraints.
We considered a dataset where a total of mutated nucleotides in the V(D)J region of the Ig locus were measured in sequences extracted from the same germinal center (see Methods) [66]. This dataset was particularly interesting because of the high number of mutations found and because of the presence of polyallelic mutations.
We mathematically mixed the multisubclonal data and applied our TrAp algorithm taking into account that the SHM scenario consists of nonbinary aberrations. In all simulations, TrAp was able to recover the original sequences and the solution was unique in of the simulations. The TrApsolution of one simulation is shown in Figure 11. However, even if the solution was not always unique, in more than of the simulations there were only five or less candidate solutions satisfying the evolutionary, parsimony and sparsity constraints, all of which correctly identified six out of seven subclones.
In addition to the identification of the underlying subclones, the TrAp algorithm generates evolutionary trees, which represent the B cell lineage during the somatic hypermutation process. The reconstruction of B cell lineage can provides important insights into the mechanisms that regulate adaptive immunity. B cell lineage reconstruction is generally performed using maximum parsimony constraints [56] using the sequences of several Ig loci as input. However, in contrast to these approaches, the TrAp algorithm is able to generate maximum parsimony trees when only the relative frequency of mutations at each nucleotide is available. Therefore, the TrAp algorithm can be used to generate parsimonious evolutionary trees when only partial sequence information is available, e.g. when only short read sequences from a single aggregate sample are available, or when the loci analyzed span a region that is too large to be fully sequenced.
Analysis of tumor biopsies
Comparison between subclonal aberration profiles inferred from heterogeneous cell populations and singlcell aberration profiles
We analyzed data from a recent study on renal cell carcinoma where two aggregate samples and twenty single cells were isolated from a clear cell renal cell carcinoma (ccRCC) and subjected to exome sequencing. Interestingly, the original study only showed partial similarity between the single cells and the aggregate [24]. However, since the single cells and the aggregate used in the experiments are from the same tumor, we sought to investigate whether any subclones inferred by TrAp would share a similar combination of mutations found in any of the single cells.
We applied our TrAp algorithm to the aggregate sample and obtained an evolutionary tree consisting of three main subclones. Due to the lack of extensive validations, we limited ourselves to investigate whether mutations that cooccur in the TrApsolution also cooccur in single cell samples. We considered the mutations that were validated by bioinformatics analysis (Table S3A from Xu et al. [24]) and by PCR validation (TableS3B). The fraction of correctly estimated cooccurrences was for mutations validated by bioinformatics analysis and for mutations validated by PCR.
Melanoma data
Finally, we applied our algorithm to investigate evolutionary mechanisms in tumor metastases using exome sequencing data from three tumor metastases (TM1: left lateral chest wall, TM4: pleural cavity, and TM3: right axilla) a matched normal sample (N: left lateral chest wall) of one melanoma patient [67]. TrAp can efficiently handle aggregate signal vectors of about unique frequencies and therefore we perform deconvolution analysis only on one chromosome. We selected chromosome as it contains the tumor suppressor Deleted in Colorectal Cancer (DCC) gene, which is known to exhibit a high load of mutations only in melanoma [68], in contrast to low expression, loss of heterozygosity or epigenetic silencing in other tumors.
To apply the TrAp algorithm, we first preprocessed the data and selected mutations on chromosome (see Methods). We labeled each mutation according to the gene affected and the amino acid change caused by the mutation (e.g. the label DCC.L1099H indicates a mutation in the DCC gene that causes a mutation from a Leucine to a Histidine at position in the DCC protein). There are six mutations with frequency in all samples (including the normal): ADNP2.G54G, ALPK2.I2157V, CD226.S307G, DCC.F23L, NETO1.S481N, and SLC39A6.E119D. The only other mutation found in the normal sample was TCEB3CL.S301C, which occurs with frequency in all samples. Moreover, the mutations ALPK2.R136S, CHST9.S122N, FAM38B.V2463, LAMA1.S1577A, LAMA1.K2002E, MYOM1.T215M, SERPINB10.R246C and SLC14A2.A880T were found in all three tumor samples and shared a similar frequency profile. The mutations DSC3.A28D, DSG1.M11V and IMPACT.D125E were found only in the metastases samples TM3 and TM4 and shared a similar frequency profile. Finally, the mutation DCC.L1099H was found only in the sample TM3.
Independent runs on the three metastatic samples gave optimal solutions for TM1, for TM3 and only TrApsolution for TM4. These high number of solutions is due to the substantial noise of the experiment (in the range ) and the fact that in samples TM1 and TM3, two of the subclones have similar frequencies and are thus difficult to separate from one another. However, TrAp identified an unique solution in the sample TM4, where the three populated subclones are distributed with significantly different frequencies (Figure 12 middle). Next, we reasoned that the metastatic TM1, TM3 and TM4 samples may share common ancestors and that their evolutionary profiles may be related. We then applied our TrAp algorithm while also imposing that all evolutionary trees must be a subset of one global evolutionary tree. This approach returned an unique solution for each sample, all of which were among the solution sets identified in the previous analyses. We observe that this approach can be very powerful since, in principle, it allows the reconstruction of very large trees by combining several snapshots of the related subclonal populations.
The results of the deconvolution are shown in Figure 12. We observe the presence of two main branches. Mutations in the left branch of TM4 (77%) are more abundant than in TM1 and TM3 (48%). We note that LAMA1, a protein that is involved in cell adhesion, is present in the right branch. 44% (19%) of the subclones of TM3 (TM4) have mutations in DSC3, DSG1 and IMPACT. The TM3 subclone also acquires a second mutation in the DCC gene (DCC.L1099HL) in addition to the mutation DCC.F23L, which was hereditary. The novel mutation in the DCC gene occurs close to the boundary between the extracellular domain and the transmembrane domain of the protein product. The resulting Histidine aminoacid is positively charged, opposed to the Leucine aminoacid of wildtype, which is neutrally charged. Because this change is next to the cell membrane, it may have repercussions on the functionality of the DCC protein product, perhaps causing inactivation, similar to the inactivation caused by loss of heterozygosity and transcript suppression observed in other cancer types.
Discussion
In the present study we presented the TrAp algorithm, a tool for inferring subclonal composition and abundance from a single aggregate measurement experiment. As we have shown, TrAp is robust to noise and it can deconvolve mixtures where multiple mutations occurs at the same locus. TrAp efficiently enumerates all possible solutions that are compatible with the model constraints, thus always identifying the sparsest and most parsimonious solution(s). However, TrAp will also generate trees (cf. Equation 2) in cases where no tree structure can be inferred. As we have shown, such structures, while deviating from the true underlying population structure, can still capture relevant cooccurrences of mutations that are specific to certain subclones. Further, in contrast to parsimonious neighborjoining approaches, which rely on sampling single subclones from the population (e.g. single cell experiments), TrAp utilizes one single aggregate experiment as input, thus overcoming the issue of small sampling size, which may be insufficient to cover the whole spectrum of subclones in a sample. We successfully deconvolved systems of up to aberrations. Although this number may not always be large enough to consider all somatic mutations found in a tumor sample, this problem can be circumvented by clustering aberrations with similar frequencies before running the TrAp algorithm.
The level of characterization achieved by subclonal deconvolution holds high potential for personalized therapies. Possible applications include the classification of subclones in primary tumors, the identification of the seeds of metastases, tracing of resistant subclones especially after drug treatments and developing treatment strategies to eliminate resistant subclones. Furthermore, our proposed model can be applied to other medical problems, such as tracing bacterial or viral paths of adaptation within the infected host, detailed genomewide reconstruction of the epigenetic differentiation program, or class specification in the hematopoietic system or in other systems.
Methods
In this section, we describe the procedures used to preprocess input data for the TrAp algorithm.
Data processing and generation
Random in silico evolutionary trees
We performed simulations by sampling genotypes whose size ranged from to and with underlying number of populated subclones ranging from to . The simulations were repeated for measurement error values equal to , and . For each combination of these quantities, we performed runs using in silico generated data as follows: during each run, a random evolutionary tree with aberrations was generated. The set of populated subclones was then selected by first including all leaves of the generated tree and then adding the remaining subclones randomly. The frequency of the populated subclones was randomly assigned and the frequency of the nonpopulated clones was set to zero. Next, the aggregate frequency vector was calculated from the generated tree. Finally, we perturbed each element of by adding an error drawn from a uniform distribution . The elements and the error are used as input for the bound error model option of the TrAp algorithm.
Cytogenetic data
The cytogenetic data was obtained from the Mitelman database, which contains biopsies as of August 15, 2012. We accessed the database through the CGAP Website[69] and we filtered out biopsies with uncertain calls (indicated by a ”?” in the karyotype data). We focused our search only on aberrations that are binary by nature, namely chromosome deletions and translocations. For each biopsy, we performed in silico simulations in which we mixed the subclones using random nonnegative coefficients.
Somatic hypermutation data
SHM sequencing data was derived from B cells undergoing somatic hypermutation (SHM), a process that leads to highaffinity antibody molecules [70]. In details, we analyzed sequences of the V(D)J region extracted from the same germinal center from the sample 11930d16_4print.2, which was sequenced by Anderson et al. [66, 71]. The sequences were aligned using the IMGT alignment tool [72, 73] in order to properly align the V, D and J regions. We selected sequences that were aligned to the same V(D)J sequence (). Since these sequences are from the same germinal center and align to the same V(D)J sequence, they are expected to stem from a single naïve B cell and have evolved through the somatic hypermutation process. From the sequencing experiment, mutated nucleotides were identified. Furthermore, a polyallelic mutations were found at position , where both and mutations were observed. Next, we considered the unique sequences (sequences and were identical) as representatives of the genome of different subclones. Similarly to the previous application, we mixed these subclones using random nonnegative coefficients and performed simulations using random nonnegative coefficients.
Exome capture sequencing data
Exomecapture data [74] was obtained from a recent clear cell renal cell carcinoma (ccRCC) study [24] and from the melanoma patient YUHUY of the Yale cohort, for which DNA from normal circulating lymphocytes and three metastases (TM1, TM3 and TM4) were subjected to exomecapture Illumina HiSeq sequencing [67].
ExomeSeq reads from the aggregate samples of the ccRCC patient were combined and aligned to the human reference genome (assembly hg19) using Bowtie [75] with parameters ”n3 k1 m20 l32 –chunkmbs 1024 –best –strata”. The frequency and coverage of each point mutation was computed using VarScan [76]. We further selected the mutations that were validated by Xu et al. [24] by PCR validation (TableS3B) and by bioinformatics analysis (Table S3A ), whose genomic coordinated were realigned to the assembly hg19 using the LiftOver tool of Galaxy [77].
For the melanoma patient YUHUY [67], we selected mutations that were homozygous in the normal sample (i.e. or ), had high sequence coverage (i.e. more than reads covering the specific nucleotide) and were localized on chromosome .
Since none of the genomic positions analyzed contained polyallelic mutations, we assign a binary state (normal/mutated) to each selected genomic position and we estimated the aggregate signal and measurement error for each mutation event using a normal approximation as
(11) 
where is total the number of reads covering position and is the number of reads covering position where the nucleotide is mutated. Finally, the and vectors were used as input for the TrAp algorithm, which was set to use the normal error model.
Implementation of the TrAp algorithm
TrAp was programmed in Java. TrAp makes use of the Java Matrix package [78] for linear regression and code by Josh Vermaas to solve the nonnegative least square problem using JAMA. The Java Universal Network/Graph Framework (JUNG) [79] is used for creating pictorial representations of evolutionary trees. TrAp is released under the GNU Lesser General Public License 3.0 and can be downloaded on the SourceForge repository of Kluger’s lab at the URL http://sourceforge.net/projects/klugerlab/files/TrAp/.
Acknowledgments
We are thankful to Ruth Halaban, Michael Krauthammer and Perry Evans for providing the melanoma sequencing data. We are grateful to Steven Kleinstein, Mohamed Uduman and Gur Yaari for providing the somatic hypermutation data and for helpful discussions. Finally, we thank Jiaqi Jin and Kelly Stanton for critical readings of the manuscript.
References
 P. Nowell, “The clonal evolution of tumor cell populations,” Science, vol. 194, no. 4260, pp. 23–28, 1976.
 M. Greaves and C. C. Maley, “Clonal evolution in cancer,” Nature, vol. 481, pp. 306–313, 01 2012.
 K. Anderson, C. Lutz, F. W. van Delft, C. M. Bateman, Y. Guo, S. M. Colman, H. Kempski, A. V. Moorman, I. Titley, J. Swansbury, L. Kearney, T. Enver, and M. Greaves, “Genetic variegation of clonal architecture and propagating cells in leukaemia,” Nature, vol. 469, pp. 356–361, 01 2011.
 J. Cairns, “Mutation selection and the natural history of cancer,” Sci Aging Knowledge Environ, vol. 2006, no. 10, p. cp1, 2006.
 C. A. Klein, “Parallel progression of primary tumours and metastases,” Nat Rev Cancer, vol. 9, pp. 302–312, 04 2009.
 M. Gerlinger, A. J. Rowan, S. Horswell, J. Larkin, D. Endesfelder, E. Gronroos, P. Martinez, N. Matthews, A. Stewart, P. Tarpey, I. Varela, B. Phillimore, S. Begum, N. Q. McDonald, A. Butler, D. Jones, K. Raine, C. Latimer, C. R. Santos, M. Nohadani, A. C. Eklund, B. SpencerDene, G. Clark, L. Pickering, G. Stamp, M. Gore, Z. Szallasi, J. Downward, P. A. Futreal, and C. Swanton, “Intratumor heterogeneity and branched evolution revealed by multiregion sequencing,” N Engl J Med, vol. 366, pp. 883–892, Mar 2012.
 P. J. Campbell, E. D. Pleasance, P. J. Stephens, E. Dicks, R. Rance, I. Goodhead, G. A. Follows, A. R. Green, P. A. Futreal, and M. R. Stratton, “Subclonal phylogenetic structures in cancer revealed by ultradeep sequencing,” Proc Natl Acad Sci U S A, vol. 105, no. 35, pp. 13081–13086, 2008.
 D. Hanahan and R. A. Weinberg, “Hallmarks of cancer: The next generation,” Cell, vol. 144, no. 5, pp. 646–674, 2011.
 O. Podlaha, M. Riester, S. De, and F. Michor, “Evolution of the cancer genome,” Trends Genet, vol. 28, pp. 155–163, 04 2012.
 S. P. Shah, A. Roth, R. Goya, A. Oloumi, G. Ha, Y. Zhao, G. Turashvili, J. Ding, K. Tse, G. Haffari, A. Bashashati, L. M. Prentice, J. Khattra, A. Burleigh, D. Yap, V. Bernard, A. McPherson, K. Shumansky, A. Crisan, R. Giuliany, A. HeraviMoussavi, J. Rosner, D. Lai, I. Birol, R. Varhol, A. Tam, N. Dhalla, T. Zeng, K. Ma, S. K. Chan, M. Griffith, A. Moradian, S. W. G. Cheng, G. B. Morin, P. Watson, K. Gelmon, S. Chia, S.F. Chin, C. Curtis, O. M. Rueda, P. D. Pharoah, S. Damaraju, J. Mackey, K. Hoon, T. Harkins, V. Tadigotla, M. Sigaroudinia, P. Gascard, T. Tlsty, J. F. Costello, I. M. Meyer, C. J. Eaves, W. W. Wasserman, S. Jones, D. Huntsman, M. Hirst, C. Caldas, M. A. Marra, and S. Aparicio, “The clonal and mutational evolution spectrum of primary triplenegative breast cancers,” Nature, vol. 486, pp. 395–399, 06 2012.
 C. G. Mullighan, L. A. Phillips, X. Su, J. Ma, C. B. Miller, S. A. Shurtleff, and J. R. Downing, “Genomic analysis of the clonal origins of relapsed acute lymphoblastic leukemia,” Science, vol. 322, no. 5906, pp. 1377–1380, 2008.
 F. Notta, C. G. Mullighan, J. C. Y. Wang, A. Poeppl, S. Doulatov, L. A. Phillips, J. Ma, M. D. Minden, J. R. Downing, and J. E. Dick, “Evolution of human BCR—ABL1 lymphoblastic leukaemiainitiating cells,” Nature, vol. 469, pp. 362–367, 01 2011.
 G. Driessens, B. Beck, A. Caauwe, B. D. Simons, and C. Blanpain, “Defining the mode of tumour growth by clonal analysis,” Nature, vol. 488, pp. 527—530, 08 2012.
 J. Chen, Y. Li, T.S. Yu, R. M. McKay, D. K. Burns, S. G. Kernie, and L. F. Parada, “A restricted cell population propagates glioblastoma growth after chemotherapy,” Nature, vol. 488, pp. 522–526, 08 2012.
 A. G. Schepers, H. J. Snippert, D. E. Stange, M. van den Born, J. H. van Es, M. van de Wetering, and H. Clevers, “Lineage tracing reveals Lgr5 stem cell activity in mouse intestinal adenomas,” Science, vol. 337, no. 6095, pp. 730–735, 2012.
 L. A. Loeb, “Human cancers express mutator phenotypes: origin, consequences and targeting,” Nat Rev Cancer, vol. 11, pp. 450–457, 06 2011.
 I. J. Fidler and M. L. Kripke, “Metastasis results from preexisting variant cells within a malignant tumor,” Science, vol. 197, no. 4306, pp. 893–895, 1977.
 R. Bouabdallah, P. Abéna, B. Chetaille, T. AurranSchleinitz, D. Sainty, P. Dubus, C. Arnoulet, D. Coso, L. Xerri, and J.A. Gastaut, “True histiocytic lymphoma following Bacute lymphoblastic leukaemia: case report with evidence for a common clonal origin in both neoplasms,” Br J Haematol, vol. 113, no. 4, pp. 1047–1050, 2001.
 A. L. Feldman, F. Berthold, R. J. Arceci, C. Abramowsky, B. M. Shehata, K. P. Mann, S. J. Lauer, J. Pritchard, M. Raffeld, and E. S. Jaffe, “Clonal relationship between precursor Tlymphoblastic leukaemia/lymphoma and Langerhanscell histiocytosis,” Lancet Oncol, vol. 6, no. 6, pp. 435–437, 2005.
 S. P. Shah, R. D. Morin, J. Khattra, L. Prentice, T. Pugh, A. Burleigh, A. Delaney, K. Gelmon, R. Guliany, J. Senz, C. Steidl, R. A. Holt, S. Jones, M. Sun, G. Leung, R. Moore, T. Severson, G. A. Taylor, A. E. Teschendorff, K. Tse, G. Turashvili, R. Varhol, R. L. Warren, P. Watson, Y. Zhao, C. Caldas, D. Huntsman, M. Hirst, M. A. Marra, and S. Aparicio, “Mutational evolution in a lobular breast tumour profiled at single nucleotide resolution,” Nature, vol. 461, pp. 809–813, 10 2009.
 I. J. Fidler and M. L. Kripke, “Genomic analysis of primary tumors does not address the prevalence of metastatic cells in the population,” Nat Genet, vol. 34, pp. 23–23, 05 2003.
 N. Navin, J. Kendall, J. Troge, P. Andrews, L. Rodgers, J. McIndoo, K. Cook, A. Stepansky, D. Levy, D. Esposito, L. Muthuswamy, A. Krasnitz, W. R. McCombie, J. Hicks, and M. Wigler, “Tumour evolution inferred by singlecell sequencing,” Nature, vol. 472, pp. 90–94, 04 2011.
 J. M. Irish, R. Hovland, P. O. Krutzik, O. D. Perez, Ø. Bruserud, B. T. Gjertsen, and G. P. Nolan, “Single cell profiling of potentiated phosphoprotein networks in cancer cells,” Cell, vol. 118, no. 2, pp. 217–228, 2004.
 X. Xu, Y. Hou, X. Yin, L. Bao, A. Tang, L. Song, F. Li, S. Tsang, K. Wu, H. Wu, W. He, L. Zeng, M. Xing, R. Wu, H. Jiang, X. Liu, D. Cao, G. Guo, X. Hu, Y. Gui, Z. Li, W. Xie, X. Sun, M. Shi, Z. Cai, B. Wang, M. Zhong, J. Li, Z. Lu, N. Gu, X. Zhang, L. Goodman, L. Bolund, J. Wang, H. Yang, K. Kristiansen, M. Dean, Y. Li, and J. Wang, “Singlecell exome sequencing reveals singlenucleotide mutation characteristics of a kidney tumor,” Cell, vol. 148, no. 5, pp. 886–895, 2012.
 Y. Hou, L. Song, P. Zhu, B. Zhang, Y. Tao, X. Xu, F. Li, K. Wu, J. Liang, D. Shao, H. Wu, X. Ye, C. Ye, R. Wu, M. Jian, Y. Chen, W. Xie, R. Zhang, L. Chen, X. Liu, X. Yao, H. Zheng, C. Yu, Q. Li, Z. Gong, M. Mao, X. Yang, L. Yang, J. Li, W. Wang, Z. Lu, N. Gu, G. Laurie, L. Bolund, K. Kristiansen, J. Wang, H. Yang, Y. Li, X. Zhang, and J. Wang, “Singlecell exome sequencing and monoclonal evolution of a JAK2negative myeloproliferative neoplasm,” Cell, vol. 148, no. 5, pp. 873–885, 2012.
 S. NikZainal, P. V. Loo, D. C. Wedge, L. B. Alexandrov, C. D. Greenman, K. W. Lau, K. Raine, D. Jones, J. Marshall, M. Ramakrishna, A. Shlien, S. L. Cooke, J. Hinton, A. Menzies, L. A. Stebbings, C. Leroy, M. Jia, R. Rance, L. J. Mudie, S. J. Gamble, P. J. Stephens, S. McLaren, P. S. Tarpey, E. Papaemmanuil, H. R. Davies, I. Varela, D. J. McBride, G. R. Bignell, K. Leung, A. P. Butler, J. W. Teague, S. Martin, G. Jönsson, O. Mariani, S. Boyault, P. Miron, A. Fatima, A. Langerød, S. A. Aparicio, A. Tutt, A. M. Sieuwerts, Å. Borg, G. Thomas, A. V. Salomon, A. L. Richardson, A.L. BørresenDale, P. A. Futreal, M. R. Stratton, and P. J. Campbell, “The life history of 21 breast cancers,” Cell, vol. 149, no. 5, pp. 994–1007, 2012.
 N. Navin and J. Hicks, “Future medical applications of singlecell sequencing in cancer,” Genome Med, vol. 3, no. 5, p. 31, 2011.
 L. D. Sjö, C. B. Poulsen, M. Hansen, M. B. Møller, and E. Ralfkiaer, “Profiling of diffuse large Bcell lymphoma by immunohistochemistry: identification of prognostic subgroups,” Eur J Haematol, vol. 79, no. 6, pp. 501–507, 2007.
 M. Varma and B. Jasani, “Diagnostic utility of immunohistochemistry in morphologically difficult prostate cancer: review of current literature,” Histopathology, vol. 47, no. 1, pp. 1–16, 2005.
 N. Yamamoto, M. Yang, P. Jiang, M. Xu, H. Tsuchiya, K. Tomita, A. R. Moossa, and R. M. Hoffman, “Determination of clonality of metastasis by cellspecific colorcoded fluorescentprotein imaging,” Cancer Res, vol. 63, no. 22, pp. 7785–7790, 2003.
 C. S.O. Attolini and F. Michor, “Evolutionary theory of cancer,” Ann N Y Acad Sci, vol. 1168, no. 1, pp. 23–51, 2009.
 N. Navin, A. Krasnitz, L. Rodgers, K. Cook, J. Meth, J. Kendall, M. Riggs, Y. Eberling, J. Troge, V. Grubor, D. Levy, P. Lundin, S. Månér, A. Zetterberg, J. Hicks, and M. Wigler, “Inferring tumor progression from genomic heterogeneity,” Genome Res, vol. 20, no. 1, pp. 68–80, 2010.
 R. Halaban, W. Zhang, A. Bacchiocchi, E. Cheng, F. Parisi, S. Ariyan, M. Krauthammer, J. P. McCusker, Y. Kluger, and M. Sznol, “PLX4032, a selective kinase inhibitor, activates the erk pathway and enhances cell migration and proliferation of melanoma cells,” Pigment Cell Melanoma Res, vol. 23, no. 2, pp. 190–200, 2010.
 S. Banerji, K. Cibulskis, C. RangelEscareno, K. K. Brown, S. L. Carter, A. M. Frederick, M. S. Lawrence, A. Y. Sivachenko, C. Sougnez, L. Zou, M. L. Cortes, J. C. FernandezLopez, S. Peng, K. G. Ardlie, D. Auclair, V. BautistaPina, F. Duke, J. Francis, J. Jung, A. MaffuzAziz, R. C. Onofrio, M. Parkin, N. H. Pho, V. QuintanarJurado, A. H. Ramos, R. RebollarVega, S. RodriguezCuevas, S. L. RomeroCordoba, S. E. Schumacher, N. Stransky, K. M. Thompson, L. UribeFigueroa, J. Baselga, R. Beroukhim, K. Polyak, D. C. Sgroi, A. L. Richardson, G. JimenezSanchez, E. S. Lander, S. B. Gabriel, L. A. Garraway, T. R. Golub, J. MelendezZajgla, A. Toker, G. Getz, A. HidalgoMiranda, and M. Meyerson, “Sequence analysis of mutations and translocations across breast cancer subtypes,” Nature, vol. 486, pp. 405–409, 06 2012.
 H. Matsushita, M. D. Vesely, D. C. Koboldt, C. G. Rickert, R. Uppaluri, V. J. Magrini, C. D. Arthur, J. M. White, Y.S. Chen, L. K. Shea, J. Hundal, M. C. Wendl, R. Demeter, T. Wylie, J. P. Allison, M. J. Smyth, L. J. Old, E. R. Mardis, and R. D. Schreiber, “Cancer exome analysis reveals a Tcelldependent mechanism of cancer immunoediting,” Nature, vol. 482, pp. 400–404, 02 2012.
 I. Varela, P. Tarpey, K. Raine, D. Huang, C. K. Ong, P. Stephens, H. Davies, D. Jones, M.L. Lin, J. Teague, G. Bignell, A. Butler, J. Cho, G. L. Dalgliesh, D. Galappaththige, C. Greenman, C. Hardy, M. Jia, C. Latimer, K. W. Lau, J. Marshall, S. McLaren, A. Menzies, L. Mudie, L. Stebbings, D. A. Largaespada, L. F. A. Wessels, S. Richard, R. J. Kahnoski, J. Anema, D. A. Tuveson, P. A. PerezMancera, V. Mustonen, A. Fischer, D. J. Adams, A. Rust, W. Chanon, C. Subimerb, K. Dykema, K. Furge, P. J. Campbell, B. T. Teh, M. R. Stratton, and P. A. Futreal, “Exome sequencing identifies frequent mutation of the SWI/SNF complex gene PBRM1 in renal carcinoma,” Nature, vol. 469, pp. 539–542, 01 2011.
 A. Hyvärinen, “Fast and robust fixedpoint algorithms for independent component analysis,” IEEE Trans Neural Netw, vol. 10, no. 3, pp. 626–634, 1999.
 H. Attias, “Independent factor analysis,” Neural Comput, vol. 11, no. 4, pp. 803–851, 1999.
 J. J. Li, C.R. Jiang, J. B. Brown, H. Huang, and P. J. Bickel, “Sparse linear modeling of nextgeneration mRNA sequencing (RNASeq) data for isoform discovery and abundance estimation,” Proc Natl Acad Sci U S A, vol. 108, no. 50, pp. 19867–19872, 2011.
 T. Gong, N. Hartmann, I. S. Kohane, V. Brinkmann, F. Staedtler, M. Letzkus, S. Bongiovanni, and J. D. Szustakowski, “Optimal deconvolution of transcriptional profiling data using quadratic programming with application to complex clinical blood samples,” PLoS One, vol. 6, p. e27156, 11 2011.
 D. Repsilber, S. Kern, A. Telaar, G. Walzl, G. F. Black, J. Selbig, S. K. Parida, S. H. E. Kaufmann, and M. Jacobsen, “Biomarker discovery in heterogeneous tissue samples taking the insilico deconfounding approach,” BMC Bioinformatics, vol. 11, no. 1, p. 27, 2010.
 T. Erkkilä, S. Lehmusvaara, P. Ruusuvuori, T. Visakorpi, I. Shmulevich, and H. Lähdesmäki, “Probabilistic analysis of gene expression measurements from heterogeneous tissues,” Bioinformatics, vol. 26, no. 20, pp. 2571–2577, 2010.
 G. Quon and Q. Morris, “ISOLATE: a computational strategy for identifying the primary origin of cancers using highthroughput sequencing,” Bioinformatics, vol. 25, no. 21, pp. 2882–2889, 2009.
 S. S. ShenOrr, R. Tibshirani, P. Khatri, D. L. Bodian, F. Staedtler, N. M. Perry, T. Hastie, M. M. Sarwal, M. M. Davis, and A. J. Butte, “Cell typespecific gene expression differences in complex tissues,” Nat Methods, vol. 7, no. 4, pp. 287–289, 2010.
 J. P. Huelsenbeck, F. Ronquist, R. Nielsen, and J. P. Bollback, “Bayesian inference of phylogeny and its impact on evolutionary biology,” Science, vol. 294, no. 5550, pp. 2310–2314, 2001.
 J.P. Doyon, V. Ranwez, V. Daubin, and V. Berry, “Models, algorithms and programs for phylogeny reconciliation,” Brief Bioinform, vol. 12, no. 5, pp. 392–400, 2011.
 S.C. Chen and B. G. Lindsay, “Building mixture trees from binary sequence data,” Biometrika, vol. 93, no. 4, pp. 843–860, 2006.
 L. Kannan and W. C. Wheeler, “Maximum parsimony on phylogenetic networks,” Algorithms Mol Biol, vol. 7, no. 1, p. 9, 2012.
 A. von Heydebreck, B. Gunawan, and L. Füzesi, “Maximum likelihood estimation of oncogenetic tree models,” Biostatistics, vol. 5, no. 4, pp. 545–556, 2004.
 R. Desper, F. Jiang, O. P. Kallioniemi, H. Moch, C. H. Papadimitriou, and A. A. Schäffer, “Inferring tree models for oncogenesis from comparative genome hybridization data,” J Comput Biol, vol. 6, no. 1, pp. 37–51, 1999.
 R. Desper, F. Jiang, O. P. Kallioniemi, H. Moch, C. H. Papadimitriou, and A. A. Schaffer, “Distancebased reconstruction of tree models for oncogenesis,” J Comput Biol, vol. 7, no. 6, pp. 789–803, 2000.
 M. Hjelm, M. Höglund, and J. Lagergren, “New probabilistic network models and algorithms for oncogenesis,” J Comput Biol, vol. 13, no. 4, pp. 853–865, 2006.
 M. D. Radmacher, R. Simon, R. Desper, R. Taetle, A. A. Schäffer, and M. A. Nelson, “Graph models of oncogenesis with an application to melanoma,” J Theor Biol, vol. 212, no. 4, pp. 535–548, 2001.
 K. Sprouffske, J. Pepper, and C. C. Maley, “Accurate reconstruction of the temporal order of mutations in neoplastic progression,” Cancer Prev Res, vol. 4, no. 7, pp. 1135–1144, 2011.
 C. D. Greenman, E. D. Pleasance, S. Newman, F. Yang, B. Fu, S. NikZainal, D. Jones, K. W. Lau, N. Carter, P. A. Edwards, P. A. Futreal, M. R. Stratton, and P. J. Campbell, “Estimation of rearrangement phylogeny for cancer genomes,” Genome Res, vol. 22, no. 2, pp. 346–361, 2012.
 M. Clement, D. Posada, and K. A. Crandall, “TCS: a computer program to estimate gene genealogies,” Mol Ecol, vol. 9, no. 10, pp. 1657–1659, 2000.
 F. Parisi, S. Ariyan, D. Narayan, A. Bacchiocchi, K. Hoyt, E. Cheng, F. Xu, P. Li, R. Halaban, and Y. Kluger, “Detecting copy number status and uncovering subclonal markers in heterogeneous tumor biopsies,” BMC Genomics, vol. 12, p. 230, 2011.
 J. S. Farris, “Estimation of conservatism of characters by constancy within biological populations,” Evolution, vol. 20, no. 4, pp. 587–591, 1966.
 A. G. Kluge and J. S. Farris, “Quantitative phyletics and the evolution of anurans,” Syst Biol, vol. 18, no. 1, pp. 1–32, 1969.
 A. Yao and H. Rubin, “Automatic enumeration and characterization of heterogeneous clonal progression in cell transformation,” Proc Natl Acad Sci U S A, vol. 90, no. 22, pp. 10524–10528, 1993.
 E. P. Malaise, N. Chavaudra, and M. Tubiana, “The relationship between growth rate, labelling index and histological type of human solid tumours,” Eur J Cancer, vol. 9, no. 4, pp. 305–312, 1973.
 C. W. Borchardt, “Über eine interpolationsformel für eine art symmetrischer functionen und über deren anwendung,” Math. Abh. der Akademie der Wissenschaften zu Berlin, pp. 1–20, 1860.
 A. Cayley, “A theorem on trees,” Quart J Math, vol. 23, pp. 376–378, 1889.
 F. H. Thompson, J. Emerson, S. Olson, R. Weinstein, S. A. Leavitt, S. P. Leong, S. Emerson, J. M. Trent, M. A. Nelson, S. E. Salmon, and R. Taetle, “Cytogenetics of 158 patients with regional or disseminated melanoma subset analysis of neardiploid and simple karyotypes,” Cancer Genet Cytogenet, vol. 83, no. 2, pp. 93–104, 1995.
 K. Junker, A. Schlichter, U. Junker, B. Knöfel, H. Kosmehl, J. Schubert, and U. Claussen, “Cytogenetic, histopathologic, and immunologic studies of multifocal renal cell carcinoma,” Cancer, vol. 79, no. 5, pp. 975–981, 1997.
 S. M. Anderson, A. Khalil, M. Uduman, U. Hershberg, Y. Louzoun, A. M. Haberman, S. H. Kleinstein, and M. J. Shlomchik, “Taking advantage: Highaffinity B cells in the germinal center have lower death rates, but similar rates of division, compared to lowaffinity cells,” J Immunol, vol. 183, no. 11, pp. 7314–7325, 2009.
 M. Krauthammer, Y. Kong, B. H. Ha, P. Evans, A. Bacchiocchi, J. P. McCusker, E. Cheng, M. J. Davis, G. Goh, M. Choi, S. Ariyan, D. Narayan, K. DuttonRegester, A. Capatana, E. C. Holman, M. Bosenberg, M. Sznol, H. M. Kluger, D. E. Brash, D. F. Stern, M. A. Materin, R. S. Lo, S. Mane, S. Ma, K. K. Kidd, N. K. Hayward, R. P. Lifton, J. Schlessinger, T. J. Boggon, and R. Halaban, “Exome sequencing identifies recurrent somatic RAC1 mutations in melanoma,” Nat Genet, vol. 44, pp. 1006—1014, 07 2012.
 A. Bernet and J. Fitamant, “Netrin1 and its receptors in tumour growth promotion,” Expert Opin Ther Targets, vol. 12, no. 8, pp. 995–1007, 2008.
 F. Mitelman, B. Johansson, and F. Mertens, “Mitelman database of chromosome aberrations and gene fusions in cancer,” 2012.
 M. Liu, J. L. Duke, D. J. Richter, C. G. Vinuesa, C. C. Goodnow, S. H. Kleinstein, and D. G. Schatz, “Two levels of protection for the B cell genome during somatic hypermutation,” Nature, vol. 451, pp. 841–845, 02 2008.
 S. H. Kleinstein, “Computational laboratory for immunology & pathology.”
 J. Robinson, A. Malik, P. Parham, J. G. Bodmer, and S. G. E. Marsh, “IMGT/HLA database – a sequence database for the human major histocompatibility complex,” Tissue Antigens, vol. 55, no. 3, pp. 280–287, 2000.
 M. P. Lefranc, V. Giudicelli, Q. Kaas, E. Duprat, J. JabadoMichaloud, D. Scaviner, C. Ginestoux, O. Clément, D. Chaume, and G. Lefranc, “IMGT, the international ImMunoGeneTics information system,” Nucleic Acids Res, vol. 33, no. suppl 1, pp. D593–D597, 2005.
 J. S. Parla, I. Iossifov, I. Grabill, M. S. Spector, M. Kramer, and W. R. McCombie, “A comparative analysis of exome capture,” Genome Biol, vol. 12, no. 9, p. R97, 2011.
 B. Langmead, C. Trapnell, M. Pop, and S. L. Salzberg, “Ultrafast and memoryefficient alignment of short DNA sequences to the human genome,” Genome Biol, vol. 10, no. 3, p. R25, 2009.
 D. C. Koboldt, K. Chen, T. Wylie, D. E. Larson, M. D. McLellan, E. R. Mardis, G. M. Weinstock, R. K. Wilson, and L. Ding, “VarScan: variant detection in massively parallel sequencing of individual and pooled samples,” Bioinformatics, vol. 25, no. 17, pp. 2283–2285, 2009.
 B. Giardine, C. Riemer, R. C. Hardison, R. Burhans, L. Elnitski, P. Shah, Y. Zhang, D. Blankenberg, I. Albert, J. Taylor, W. Miller, W. J. Kent, and A. Nekrutenko, “Galaxy: A platform for interactive largescale genome analysis,” Genome Res, vol. 15, no. 10, pp. 1451–1455, 2005.
 J. Hicklin, C. Moler, P. Webb, R. Boisvert, B. Miller, R. Pozo, and K. Remington, “JAMA: A java matrix package,” 2000.
 D. Fisher, J. O’Madadhain, P. Smyth, S. White, and Y.B. Boey, “Analysis and visualization of network data using JUNG,” J Stat Softw, 2005.