## Abstract

Cancer is a clonal evolutionary process, caused by successive accumulation of genetic alterations providing milestones of tumor initiation, progression, dissemination, and/or resistance to certain therapeutic regimes. To unravel these milestones we propose a framework, tumor evolutionary directed graphs (TEDG), which is able to characterize the history of genetic alterations by integrating longitudinal and cross-sectional genomic data. We applied TEDG to a chronic lymphocytic leukemia (CLL) cohort of 70 patients spanning 12 years and show that: (a) the evolution of CLL follows a time-ordered process represented as a global flow in TEDG that proceeds from initiating events to late events; (b) there are two distinct and mutually exclusive evolutionary paths of CLL evolution; (c) higher fitness clones are present in later stages of the disease, indicating a progressive clonal replacement with more aggressive clones. Our results suggest that TEDG may constitute an effective framework to recapitulate the evolutionary history of tumors.

## eLife digest

A historical event is often the culmination of the preceding circumstances. The same can be said of cancer as a disease. Cancer results from genetic mutations that disrupt the normal biological processes within a cell, removing the fail-safes that prevent it from growing and reproducing uncontrollably.

Cancer is not caused by just one mutation, and once one gene is malfunctioning, other genes become much more likely to mutate. Although modern sequencing methods have revealed many of the genes that mutate in several different kinds of cancer, uncovering when each of these mutations occurs has been more difficult. Knowing when each mutation occurs could make it easier to predict how the cancer will progress and could also help target cancer treatments more effectively.

Wang, Khiabanian, Rossi et al. have devised a new method of studying the history of genetic mutations of cancer patients. This combines a ‘longitudinal’ method that looks at how mutations develop in a single tumor by taking samples from it at different times and ‘cross-sectional’ methods that make predictions based on data collected from a large number of patients. Wang, Khiabanian, Rossi et al. call this method ‘tumor evolutionary directed graphs’ (TEDG), as it produces a graph that shows how different gene mutations are related to each other. Initial tests showed that the TEDG method could accurately decipher the main chain of events in cancer evolution when used on data collected from at least 30 patients.

Wang, Khiabanian, Rossi et al. then used TEDG on data from 164 tumor samples collected over 12 years from 70 patients with chronic lymphocytic leukemia, the type of leukemia that is most widespread amongst adults in Western countries. This uncovered two separate ways that this cancer may develop, one of which has a higher risk of life-threatening complications.

Knowing which of the two ways chronic lymphocytic leukemia is progressing in a patient could help treat the disease, as each pathway responds differently to different treatments. In addition, understanding the paths that cancer progression follows could also provide early warning signals of the mutations that will occur next. This could help to develop alternative, targeted cancer treatments.

## Main text

### Introduction

Cancer is a complex, Darwinian, adaptive clonal evolutionary process, driven by the accumulation of genetic alterations that confer high proliferative and survival advantage (Merlo et al., 2006; Greaves and Maley, 2012). Recent advances in sequencing technologies have allowed uncovering the most common genetic alterations of many tumors, but the temporal order of most of these alterations is still unknown (Futreal et al., 2004; Greenman et al., 2007; Santarius et al., 2010; Lawrence et al., 2014). Temporal patterns of genetic alterations may indicate the fate of tumor progression, allowing early diagnosis of tumor subtypes and improving the choice of therapeutic strategies.

To understand the evolutionary history of tumors, several experimental and computational strategies have been used. Longitudinal strategies require samples at multiple time points spanning the clonal tumor evolution process (Fearon and Vogelstein, 1990; Campbell et al., 2010; Ding et al., 2010; Notta et al., 2011; Gerlinger et al., 2012; Turajlic et al., 2012; Landau et al., 2013). Landau et al. sampled leukemia cells from 18 CLL patients at two time points, revealing that *SF3B1* and *TP53* mutations are late events in subclonal tumor cells (Landau et al., 2013). The study of different stages of colorectal carcinogenesis showed the sequence of genetic events to be *APC*, *KRAS*, and then *TP53* (Fearon and Vogelstein, 1990). Another alternative approach is a cross-sectional strategy, which makes use of a large cohort of patients to computationally predict the preferred orders. RESIC is a stochastic process model to identify the order of mutations (Attolini et al., 2010), which successfully confirmed the results in colorectal cancer, suggesting that cross-sectional data is informative for the prediction of mutation order. However, RESIC does not consider a critical aspect of carcinogenesis that most tumors are heterogeneous (Parsons, 2011).

Following the assumption that different tumors proceed through related temporally ordered alterations, we propose to summarize tumor histories using a newly developed analytical approach that integrates the genomic information from different longitudinally characterized patients. Our method, termed tumor evolutionary directed graphs (TEDG), proceeds in two steps to ensemble in a simplified way cancer clonal evolutionary histories of large number of patients: first, by merging the evolutionary history of each patient, and second, by removing indirect relationships using spectral techniques for network deconvolution (Feizi et al., 2013). The resulting TEDG is a directed graph with nodes representing driver genes and arrows representing temporal order of gene lesions. A non-randomly distributed TEDG shows that cancer proceeds in an orchestrated fashion and indicates the main paths and the alternative routes of cancer evolution.

In this study, we have applied TEDG to study the dynamics of the acquisition of alterations in chronic lymphocytic leukemia (CLL), which represents the most common adult leukemia in Western countries (Hallek et al., 2008; Müller-Hermelink et al., 2008). CLL is an ideal model for studying clonal dynamics because it is possible to collect highly purified sequential samples over time, and its clinical course is well characterized by serial cycles of response, remissions, and relapse ending in some instances with the development of lethal complications such as chemoresistant progression or transformation into an aggressive lymphoma (Richter syndrome) (Pasqualucci et al., 2011; Zenz et al., 2012; Fabbri et al., 2013). No systematic approach has been followed to disentangle and characterize the ensemble of evolutionary histories of this disease. For this purpose, we envision a dual cross-sectional and longitudinal strategy by collecting genomic information from the most common alterations in a cohort of 70 CLL patients spanning over a period of 12 years (2001–2012).

### Results

#### Tumor Evolutionary Directed Graphs

To recapitulate and compare the history of genetic alterations in many patients, we propose a framework to infer TEDG by integrating longitudinal and cross-sectional genomic data of cancer patients. First, we reconstruct the sequential network of genetic alterations in each patient by analyzing genomic data from different time points. Specifically, the techniques of high-depth next generation sequencing (NGS) and fluorescence in situ hybridization (FISH) are separately carried out to assess the mutation allele frequency (MAF) and copy number abnormalities (CNA) of selected driver genes. To unify both types of data, and to adjust the MAF of mutations in genes with CNA, we introduce mutation cell frequency (MCF, defined as the fraction of tumor cells with a particular alteration) for quantification of genetic lesions (‘Materials and methods’, Figure 1—figure supplement 1). Based on MCF, we investigate alterations represented in at least 5% of leukemic cells (see examples of CLL patients in Figure 1—figure supplement 2). First, if a given genetic lesion is observed to be temporally earlier than another lesion, we connect them with a directed edge to represent their sequential order of development (Figure 1A). Second, we pool many sequential networks from different patients to construct an Integrated Sequential Network (ISN). Third, we infer TEDG from ISN by removing indirect associations with spectral techniques and minimal spanning tree algorithm. TEDG is the backbone of ISN, representing an optimal explanation of the mutation order across many patients (Figure 1B).

To test TEDG method and also to show how many patients are required to approximate the ground truth, we employ artificial examples by both simulating linear evolution and branching evolution of cancer, where the longitudinal data are generated by one-step Markov process and Nordling's multi-mutation model (Nordling, 1953) (Figure 2A,B, ‘Materials and methods’). For example, in a cohort of 15 patients with linear evolution, we start the Markov process of each case from no mutations at time zero. Mutation status at three time points 10, 20, and 30 are then successively updated based on the Markov chain transition probability described in ‘Materials and methods’. We finally pool all temporal information of those 15 patients to generate ISN (left panel of Figure 2D). The TEDG is further deduced by deconvoluting ISN. Suppose *G*_{dir} is the adjacent matrix of all direct interactions/orders, the observed ISN should be a summary of direct and indirect orders. The deconvolution is formulated by *G*_{dir} = *G*_{obs}(*I* + *βG*_{obs})^{−1}, where, *G*_{obs} is the observed weighted adjacent matrix, *I* is the identify matrix, and *β* is the scaling factor between zero and one, indicating the degree of deconvolution. To evaluate TEDG, we define the accuracy by how frequently its results match the input model. To calculate the accuracy, we generate 10 simulated datasets of *N* patients (with three sequential samples per patient). In each simulation, we apply the TEDG analysis to reconstruct the sequential order and compare the result with the input model. We then apply the concept of accuracy to find a reasonable *β* for our simulation, by calculating the accuracy of TEDG with different *β* ∈ (0, 1) when *N* ∈ {1, 2, …, 100}. Figure 2C indicates that the optimal value of *β* is related to the number of samples, but the wide range of high accuracy region suggests that the TEDG results are robust to the parameter selection. With 30 patients, TEDG's accuracy for a linear model is 80%, and for a branching model, it reaches 90% (Figure 2E–G, ‘Materials and methods’).

#### Tumor Evolutionary Directed Graph of CLL

In order to investigate the evolutionary history of CLL, we apply the TEDG framework to the driver genetic lesions of this leukemia. We study the most common alterations of 164 temporally sequential samples from 70 patients by high-depth next generation sequencing (NGS) and fluorescence in situ hybridization (FISH) analysis (‘Materials and methods’). Half (35 out of 70) of the patients have at least one subclonal genetic lesion with mutation frequency less than 20% at diagnosis (Figure 3—figure supplement 1). We firstly build sequential networks for each patient based on this longitudinal genetic data. Then, we pool all sequential networks to construct the ISN of CLL. We ask whether the genetic lesions in CLL are temporally ordered or randomly accumulated and reason that if the genetic alterations driving CLL progression follow a preferential order, there exists a well ordered directed flow in ISN. Figure 3A represents a hierarchical layout of ISN, which depicts an ordered structure of lesions in genes represented by sources (nodes with more outgoing arrows) and sinks (nodes with more incoming arrows) (p-value < 0.0001, chi square distribution).

To understand the evolutionary pattern of genetic alterations in CLL, we infer TEDG by removing indirect orders in ISN. The TEDG of CLL (Figure 3B) shows a clear pattern of the flow of alterations, revealing that genetic alterations of CLL are sequentially ordered in a branching mode. To statistically assess the temporal pattern of the genetic lesions, we use the binomial test to assess significance by assuming that the number of in-degree and out-degree of each alteration are randomly distributed. Consistent with TEDG, this analysis suggests the following temporal order of the lesion: mutations of *MYD88*, deletion 13q14, +12, mutations of *NOTCH1*, *RB1* deletion, 17p13 deletion, mutations of *SF3B1*, 11q22-q23 deletion, *BIRC3* deletion, mutations of *TP53*, and mutations of *BIRC3* (Figure 3C). Deletion of 13q14, +12 and mutations of *NOTCH1* are significant early events in CLL development, while mutations of *TP53* and *BIRC3* are significant late events. Though the sample size prevents statistical considerations, *MYD88* mutations might occur even before 13q14 deletion (one of six patients), indicating that activation of the Toll-like receptor pathway could be important in the initiation of a fraction of CLL tumors (Arvaniti et al., 2011). The sequential order of the genetic lesions remains consistent (Pearson's correlation = 0.9, p-value < 1e-3) when evolutionary networks are constructed with patients who received chemotherapy (Figure 3—figure supplement 2), suggesting that the order of development of the genetic lesions may be affected by therapy (the Pearson's correlation of the order in patients without chemotherapy is 0.4 with p-value > 0.1).

By considering each single type of mutation affecting the same gene as a distinct and independent node, we construct the comprehensive ISN of CLL mutations (Figure 4A). It is very difficult to capture useful information directly from ISN, while TEDG simplified the topology by capturing the backbone of tumor evolution (Figure 4B). We observe that the monoallelic 13q14 deletion, RB1 deletion, and +12 are significant early events (p-value < 0.01), while the *BIRC3* E537fs and the *TP53* R248Q are significant late lesions (p-value < 0.05) (Figure 4—figure supplement 1). Also, the analysis of ISN and TEDG shows that different lesions affecting the same gene may occur in distinct branches and stages. For example, mutations K700E and K666E of *SF3B1* are late events in cases harboring 13q14 deletion, while mutations R273C of *TP53* are late events in cases with +12. Though the sample size prevents statistical considerations, TEDG reveals the H179R and Y234C missense substitutions in *TP53* may be early events, while the R248Q, R273C, P152fs, and N239T substitutions are late events.

Based on pivotal NGS studies, two different evolutionary models have been proposed in CLL, namely gradual linear and branching evolution (Knight et al., 2012; Schuh et al., 2012). The analysis of our cohort (excluding patients with Richter's transformation, refer to ‘Materials and methods’ for details) shows that a minority of patients (*n* = 3/60, 5%; FDR = 0.1) are characterized by a significantly decreased or undetectable representation of the founding clone, coupled with a significant increase of a second subclone that represented a small subpopulation at an earlier time point, consistent with a branching evolution model (highlighted by yellow circle in Figure 5A–B, and Table 1). Interestingly, in all three cases clonal replacement events involve *SF3B1* mutations and occur after treatment (Figure 5C), suggesting that the branching evolution model is closely connected to the combination of treatments and the emergence of *SF3B1* mutations (p-value = 0.0016 by Fisher's exact test). This observation implies that, at the time of treatment requirement, limiting the knowledge of disease genetics to the dominant clone will likely be uninformative for accurate therapeutic decisions. Of particular interest in this scenario is the development of therapeutic strategies to prevent the branching evolution of the tumor, with the goal of eradicating dominant as well as minor clones (Anderson et al., 2011; Notta et al., 2011; Ding et al., 2012; Egan et al., 2012; Keats et al., 2012; Walker et al., 2012; Rossi et al., 2014).

Patient ID | Sampling time | Treatment | Increased alterations | Decreased alterations | Detectable alterations |
---|---|---|---|---|---|

#1 | 2001 | None | – | – | SF3B1 (K666E), SF3B1 (K700E), del13q |

2005 | FC | SF3B1 (K666E) | None | SF3B1 (K666E), SF3B1 (K700E), del13q | |

2008 | FCR/CAM/RBEN | SF3B1 (K700E), del11q | None | SF3B1 (K666E), SF3B1 (K700E), del13q, del11q | |

2011 | BENCAM | SF3B1 (K700E), del11q, del13q, delBIRC3 | SF3B1 (K666E) | SF3B1 (K700E), del13q, del11q, delBIRC3 | |

#14 | 2004 | RF | – | – | +12, del13q |

2007 | FC/CAM | None | None | +12, del13q | |

2010 | BENCAM | TP53 (R248Q), del17p, SF3B1 (K666E) | None | +12, del13q, TP53 (R248Q), del17p, SF3B1 (K666E) | |

2012 | BENDOFA | TP53 (R248Q), del17p, SF3B1 (K666E), T12 | del13q | +12, del13q, TP53 (R248Q), del17p, SF3B1 (K666E) | |

#33 | 2002 | None | – | – | NOTCH1 (P2514-) |

2004 | FCR | SF3B1 (K700E) | None | NOTCH1 (P2514-), SF3B1 (K700E) | |

2009 | FCR | SF3B1 (K700E) | NOTCH1 (P2514-) | SF3B1 (K700E) |

FC, fludarabine, cyclophosphamide; FCR, fludarabine, cyclophosphamide, rituximab; CAM, Campath; RBEN, rituximab, bendamustine; BENCAM, bendamustine, Campath; RF, rituximab, fludarabine; BENDOFA, bendamustine, ofatumumab.

#### Statistical association analysis of TEDG

To further investigate the relationship between driver genetic lesions in TEDG, we assess their associations and anti-associations in a cross-sectional cohort of 1403 CLL patients, of which 1054 cases are informative in at least one lesion (Figure 6A, Table 2). Most of the co-mutations (connected in red in Figure 6B) are experimentally confirmed by previous studies, including: co-occurrence of 17p13 deletion and *TP53* mutations, reflecting a typical two-hit model for tumor suppressor genes; co-occurrence of *BIRC3* deletion, 11q22-q23 deletion, and *BIRC3* mutations; and the relationship between *NOTCH1* mutations and +12 (Balatti et al., 2012; Del Giudice et al., 2012). In addition to previously reported associated lesions, this large cohort of patients allows the statistical power to reveal other significant and previously unreported co-occurrence interactions, including the co-occurrence of *BIRC3* abnormalities with +12 and *NOTCH1* mutations and the co-occurrence of 13q14 deletion and *BIRC3* deletion. This large cohort also reveals mutually exclusive relationships between +12 and 13q14 deletion and between *NOTCH1* mutations and 13q14 deletion (connected in blue in Figure 6B). The statistical association analysis is consistent with the prediction of TEDG and indirectly validates the ability of TEDG in successfully capturing the major information in tumor evolution.

TEDG analysis | SA analysis | |||||
---|---|---|---|---|---|---|

N | Total | % | N | Total | % | |

IGHV homology >98% | 37 | 69 | 53.6 | 532 | 1380 | 38.6 |

del13q14 | 35 | 70 | 50.0 | 682 | 1403 | 48.6 |

+12 | 19 | 70 | 27.1 | 205 | 1403 | 14.6 |

del11q22-q23 | 8 | 70 | 11.4 | 136 | 1403 | 9.7 |

del17p13 | 20 | 70 | 28.6 | 114 | 1403 | 8.1 |

TP53 mutation | 22 | 70 | 31.4 | 114 | 1401 | 8.1 |

NOTCH1 mutation | 22 | 70 | 31.4 | 150 | 1403 | 10.7 |

SF3B1 mutation | 17 | 70 | 24.3 | 100 | 1403 | 7.1 |

MYD88 mutation | 6 | 70 | 8.6 | 52 | 1403 | 3.7 |

BIRC3 mutation | 10 | 70 | 14.3 | 41 | 1403 | 2.9 |

BIRC3 deletion | 6 | 70 | 8.6 | 79 | 1403 | 5.6 |

TEDG analysis, Tumor Evolutionary Directed Graphs analysis; SA analysis, Statistical Association analysis;

*IGHV*, immunoglobulin heavy variable gene; FISH, fluorescence in situ hybridization.

By integrating the topology of TEDG to the association analysis, two distinct and mutually exclusive evolutionary paths of CLL evolution are disclosed. The first evolutionary path involves CLL harboring +12 and *NOTCH1* mutations as early driver events. In this path, clonal evolution proceeds toward the development of *TP53* and *BIRC3* abnormalities. The second evolutionary path involves 13q14 deletion as early driver lesion and proceeds toward the development of *SF3B1* mutations and *BIRC3* abnormalities. Overall, our results are consistent with the different clinico-biological phenotype of +12 CLL and 13q14 deleted CLL and further support the hypothesis that at least two distinct genetic subtypes of CLL exist.

#### Rate of allele frequency change

We reason that if a subclone replaces the major clone during tumor evolution, it is because of its higher fitness. Although fitness of a clone is not directly measurable, the growth rate of mutation frequency of drivers in this clone can serve as an indicator. Given this, mutations related to high fitness clones can be identified by extracting the alterations that rapidly change their allele/cell frequency within the tumor. We define the growth rate of each genetic lesion as the average increasing speed of mutation cell frequency per year (‘Materials and methods’). Interestingly, the growth rates of late events (17p13 deletion, mutations of *SF3B1*, 11q22-q23 deletion, *BIRC3* deletion, mutations of *TP53*, and mutations of *BIRC3*) are found higher than those of early events (mutations of *MYD88*, 13q14 deletion, +12, mutations of *NOTCH1*, and *RB1* deletion) with p-value = 0.005 by Wilcoxon Rank-Sum Test (Figure 6C), revealing that late events may drive higher fitness or more aggressive clones. Note that the initiating lesions are usually clonal, presenting in most of the tumor cells, and therefore do not increase in frequency at the same magnitude as subclonal mutations. However, after eliminating clonal genetic lesions with MCF > 20% in the above analysis, we still find the growth rates of late and early events to be significantly different with p-value = 0.0114 by Wilcoxon Rank-Sum Test (Figure 6—figure supplement 3).

We define *maximal mutation frequency slope* (MMFS) as the rate of allele frequency change of the fastest-growing clone in a patient (Table 3, ‘Materials and methods’). MMFS is a function that aims at characterizing the relative fitness of a particular clone carrying a particular mutation. In our cohort, only genetic lesions affecting the *TP53* gene show a statistically significant association with MMFS (p-value = 0.0164, by Wilcoxon Rank-Sum Test). Clinically, the fastest-growing clones are strongly correlated with poor survival and Richter syndrome transformation, consistent with the fact that CLL transformed to Richter syndrome presents significantly higher MMFS values (p-value = 0.002, by Wilcoxon Rank-Sum Test) (Figure 6—figure supplement 1). Indeed, by survival analysis, having a clone with high MMFS associates with an approximately threefold significant increase in the hazard of death (HR: 3.17; 95% CI: 1.38–7.40; p-value = 0.005) and a significant shortening of overall survival (47% at 5 years) (Figure 6—figure supplement 2). Most deaths (77%) in patients carrying a clone with high MMFS are due to Richter syndrome transformation. Consistently, patients having a clone with high MMFS show an approximately fourfold increased risk of Richter syndrome development (HR: 4.61; 95% CI: 1.54–13.80; p-value = 0.003) and an ∼50% of them are projected to develop Richter syndrome at 5 years (Figure 6—figure supplement 2).

Patient | MMFS* | Fast growing mutations | MCF1 (%) | MCF2 (%) | dT (moths) | Treatment | Richter syndrome transformation |
---|---|---|---|---|---|---|---|

#52 | 117.7 | TP53 (P152+) | 2.1 | 11.4 | <1 | None | Yes |

#68 | 9.1 | TP53 (R273C) | 0.0 | 100.0 | 11 | None | No |

#51 | 6.7 | SF3B1 (I704F) | 12.2 | 68.9 | 5 | RCVP | Yes |

#37 | 6.4 | TP53 (R248Q) | 3.8 | 90.0 | 10 | RDHAOX | Yes |

#57 | 5.9 | NOTCH1 (P2514-) | 0.8 | 94.2 | 12 | None | Yes |

#47 | 4.6 | del13q | 37.6 | 100.0 | 13 | None | No |

#14 | 3.9 | del17p | 18.9 | 90.0 | 16 | A | No |

#42 | 3.4 | TP53 (N239T) | 0.0 | 100.0 | 16 | RDHAOX | Yes |

#4 | 3.3 | del17p | 53.6 | 100.0 | 15 | CLB-O | No |

#54 | 3.0 | NOTCH1 (P2415-) | 100.0† | 100.0 | 22 | FCO | No |

#22 | 2.9 | BIRC3 deletion | 0.0 | 93.9 | 29 | FCR | No |

#20 | 2.8 | del13q | 0.0 | 100.0 | 36 | CLB | No |

#63 | 2.5 | BIRC3 (M388V) | 0.0 | 86.0 | 24 | FCR | Yes |

#6 | 2.3 | del17p | 0.0 | 92.1 | 34 | CLB | No |

#38 | 2.3 | NOTCH1 (P2514-) | 41.7 | 42.9 | 4 | CVP | Yes |

#13 | 2.2 | TP53 (G136H) | 0.0 | 100.0 | 45 | RDHAOX | No |

#1 | 2.1 | del11q | 11.8 | 100.0 | 41 | FCR/A/BR | No |

MMFS, maximal mutation frequency slope (in standard deviation per year); MCF1, mutation cell frequency of selected mutation at the first time point; MCF2, mutation cell frequency at the second time point; dT, the elapsed time between two samples; RCVP, rituximab, cyclophosphamide, vincristine, prednisone; RDHAOX, rituximab, dexamethasone, high dose cytarabine, oxaliplatin; A, alemtuzumab; CLB-O, chlorambucil, ofatumumab; FCO, fludarabine, cyclophosphamide, ofatumumab; FCR, fludarabine, cyclophosphamide, rituximab; CLB, chlorambucil; CVP, rituximab, cyclophosphamide, vincristine, prednisone; BR, bendamustine, rituximab.

Total number of the cancer cells with

*NOTCH1*alteration does not change, but the allele frequency of the mutation increases because of the deletion of the wild-type allele.

### Discussion

It is not known whether there is a preferred order of mutations in the development of cancer and how the order of mutations may impact clinical outcomes. We propose a TEDG framework, which is able to integrate longitudinal and cross-sectional genomic data into a directed graph of tumor evolution. The flow in this graph reveals underlying paths of tumor progression. Starting from time-series genomic data in one patient, a sequential network is reconstructed to capture the possible historical events of the evolution in the particular tumor. Integrating sequential data from many patients, we collect the ensemble of tumor histories by ISN, which presents a comprehensive topology of evolutionary landscape. A recent technology of network deconvolution is able to distinguish direct and indirect interactions using spectral methods, which reduce the weights of indirect connections and remove low weighted edges by selecting appropriate cut-offs (Feizi et al., 2013). To adapt this method to our particular problem depending on the number of samples, we introduce a degree parameter, β. To gain insights of the selection of beta, we propose a strategy to simulate tumor evolution by a one-step Markov process, with transition probability derived from Nordling's multi-mutation model. A linear and a branching evolution model are separately simulated in a study of four alterations. The proof-of-concept simulation shows that the TEDG strategy can obtain excellent performance in capturing the evolutionary history when the number of cases is beyond 30.

We apply TEDG to CLL and: (i) reconstruct an evolutionary network representing the sequential order of genetic lesions occurring during the course of this disease; (ii) investigate statistical associations and competitiveness of driver genetic lesions to uncover evolutionary paths; and (iii) correlate the order of alterations with the kinetics of changes in their allelic abundance, to identify the molecular alterations associated with the fastest growth of a subclone and their clinical impact on outcome.

Phylogenetic trees are often employed to infer temporal order of gene mutations by assuming that common ancestors are early events, but phylogenetic methods usually require higher number of alterations than the ones available in our cohort. Standard statistical techniques to assess the significance of different branches rely on bootstrapping segregating sites, and the statistical power of these techniques using longitudinal data from cancer patients in few selected driver genes is extremely limited. None of the standard phylogenetic methods such as distance matrix methods, Bayesian phylogenetic methods, or parsimony methods can produce branches with enough bootstrap support (>80%). On the other hand, TEDG works with aggregate data from several patients allowing statistical power for robust estimates.

According to TEDG analysis, the molecular lesions of CLL are temporally ordered in a specific fashion rather than being randomly accumulated. Among recurrent lesions, 13q14 deletion and +12 are initiating events, while mutations of *TP53* and *BIRC3* are late. This observation is consistent with the notion that del13q and +12 occur at a similar prevalence in all CLL phases, including monoclonal B-cell lymphocytosis, a condition that often anticipates overt CLL, thus suggesting that they are early events (Rawstron et al., 2008; Nieto et al., 2009; Rossi et al., 2009b; Fazi et al., 2011; Pasqualucci et al., 2011; Kern et al., 2012). Also, 13q14 deletion has been directly implicated in CLL initiation (Klein et al., 2010). On the other hand, late-onset abnormalities are known to accumulate in more advanced phases of the disease, thus suggesting that they are second-hit lesions (Zenz et al., 2009; Rossi et al., 2012). In our model, mutations of *SF3B1* and 11q22-q23 deletions appear to be acquired at an intermediate time point. This observation is consistent with the clinical observation that mutations of *SF3B1* and 11q22-q23 deletion co-segregate with an intermediate prognostic group of CLL (Döhner et al., 2000; Oscier et al., 2013; Rossi et al., 2014).

By integrating the cross-sectional data of association and anti-association between genetic lesions at CLL diagnosis with longitudinal data of clonal evolution from TEDG, two distinct and mutually exclusive evolutionary paths of evolution emerge. The first evolutionary path involves CLL initially harboring +12 and *NOTCH1* mutations. In this path, clonal evolution proceeds toward the development of *TP53* and *BIRC3* abnormalities. The second evolutionary path involves CLL initially harboring 13q14 deletion and proceeds toward the development of *SF3B1* mutations and *BIRC3* abnormalities. Our data, as well as previous analyses (Rossi et al., 2013), show that deletion of 13q14 and +12 are mutually exclusive in CLL. From a clinical standpoint, +12 CLL is known to stand out of typical 13q14 deleted CLL because of their atypical cytomorphology and phenotype, the more intense expression of CD20, the preferential nodal presentation and the higher risk of transformation to Richter syndrome. These notions along with our novel observation that +12 and 13q14 deleted CLL proceed through distinct paths of clonal evolutions further support the hypothesis that at least two distinct genetic subtypes of CLL exist.

The presence at diagnosis of fully clonal genetic lesions that are considered late genetic events in CLL evolution, such as *TP53* abnormalities, is already known to have an adverse impact on disease outcome (Zenz et al., 2008, 2010; Dicker et al., 2009; Malcikova et al., 2009; Rossi et al., 2009a; Gonzalez et al., 2011). We observe that a fraction of *TP53* abnormalities, though subclonal at presentation, lead to expansion of fitter subclones that progressively predominate with time in the tumor architecture. The clinical impact of this observation is supported by the evidence that harboring higher fitness alterations correlates with poor survival and increased risk of CLL transformation into an aggressive lymphoma, an often-lethal complication known as Richter syndrome.

In this study, we have considered some of the most important and clinically relevant drivers of CLL (Rossi et al., 2013). The *ATM* mutations are also commonly found in CLL; however, the *ATM* gene is large and highly polymorphic without well-known hotspots. Therefore, the distinction of its driver mutations from constitutional variants is challenging. Also, from a clinical standpoint, the prognostic relevance of ATM mutations in CLL is still controversial, so we have excluded this gene from the current study (Lozanski et al., 2012; Ouillette et al., 2012; Skowronska et al., 2012).

In conclusion, the application of TEDG to CLL provides the proof-of-principle that this method is able to: (i) improve cancer classification and dissection into genetic subgroups following different paths of clonal evolution; and (ii) anticipate the genetic composition of the progressive/relapsed disease according to the genetic composition of the tumor clone at the time of diagnosis, including the development of genetic lesions associated with chemorefractoriness. TEDG provides a general framework that could be used to study and compare the evolutionary histories of other tumors.

### Materials and methods

#### Samples

We collected 202 CLL patients provided with at least 2 sequential samples and followed for at least 2 years after presentation (median interval between baseline and last sequential sample: 62.8 months, range 24–150 months). 70 out of the 202 patients were informative for *TP53*, *NOTCH1*, *SF3B1*, *MYD88*, or *BIRC3* mutations (which are the most common mutations in this leukemia), overall accounting for 164 paired sequential samples collected at diagnosis, progression, and last follow-up. To analyze the dynamics of those mutations, we carried out high-depth next generation sequencing (NGS) to quantify the mutation allele frequencies at each time point of the disease course and to establish their modifications during leukemia progression. Additionally, copy number abnormalities at 13q14, chromosome 12, 11q22-q23, 17p13, as well as at the *RB1*, and *BIRC3* loci were investigated by FISH. Inclusion criteria for the longitudinal analysis of clonal evolution were: (i) having at least two years of follow-up after diagnosis; and (ii) availability of >2 sequential samples collected at: (a) diagnosis; (b) each progression requiring treatment; (c) last follow-up. CLL diagnosis was according to 2008 IWCLL-NCI criteria and confirmed by a flow cytometry score >3 in all cases. Monoclonal B-cell lymphocytosis (MBL) was excluded. The study was approved by the institutional ethical committee of the Azienda Ospedaliero-Universitaria Maggiore della Carità di Novara affiliated with the Amedeo Avogadro University of Eastern Piedmont, Novara, Italy (Protocol Code 59/CE; Study Number CE 8/11). Patients provided informed consent in accordance with local IRB requirements and Declaration of Helsinki.

CLL samples were extracted from fresh or frozen peripheral blood mononuclear cells (PBMC) isolated by Ficoll-Paque gradient centrifugation. In all cases, the fraction of tumor cells corresponded to 70–98% as assessed by flow cytometry. Matched normal DNA from the same patient was obtained from saliva or from purified granulocytes and confirmed to be tumor-free by PCR of tumor-specific IGHV-D-J rearrangements. High-molecular-weight (HMW) genomic DNA was extracted from tumor and normal samples according to standard procedures. DNA was quantified by the NanoDrop 2000C spectrophotometer (Thermo Scientific, Wilmington, DE). To validate TEDG, 1403 newly diagnosed and previously untreated CLL were enrolled in the study, of whom 931 (66%) were provided with clinical data and regular follow-up (Table 2). The frequency of alterations was higher in TEDG analysis group, because only informative patients with known mutated driver genes were included. Cross-sectional investigation of the associations and anti-associations between the most recurrent genetic lesions at diagnosis was based on the entire CLL cohort of 1403 patients. Survival analysis was based on CLL cases provided with clinical data (n = 931). Mutation screening of the *IGHV*, *TP53*, *NOTCH1*, *SF3B1*, *MYD88*, and *BIRC3* genes were performed by Sanger sequencing.

#### Estimation of the variant allele frequency

Among cases that were informative for *TP53*, *NOTCH1*, *SF3B1*, *MYD88*, or *BIRC3* mutations at any time point in the disease course, we carried out high-depth NGS to quantify the variant allele frequencies at each stage of their disease. Positions known to harbor *TP53*, *NOTCH1*, *SF3B1*, *MYD88*, or *BIRC3* mutations by Sanger sequencing were amplified from genomic DNA by oligonucleotides containing the gene-specific sequences, along with 10-bp MID tag for multiplexing and amplicon library A and B sequencing adapters. The obtained amplicon library was subjected to deep sequencing on the Genome Sequencer Junior instrument (454 Life Sciences). In order to obtain at least 700-fold coverage per amplicon, no more than 100 amplicons/run were analyzed. The obtained sequencing reads were mapped to reference sequences and analyzed by the Amplicon Variant Analyzer software (Roche) to establish the mutant allele frequency. The sequencing depth in this study is on average 1200×, which is sufficient for a highly sensitive detection of mutations with allele frequency >1% out of the background error noise (Rossi et al., 2014).

#### Fluorescence in situ hybridization (FISH)

Probes used for FISH analysis were: (i) LSID13S319 (13q14 deletion), CEP12 (trisomy 12), LSIp53 (17p13/*TP53* deletion), and LSIATM (11q2-q23/ATM deletion) (Abbott, Rome, Italy); and (ii) the RP11-177O8 (*BIRC3*) BAC clone. The labeled *BIRC3* BAC probe was tested against normal control metaphases to verify the specificity of the hybridization. For each probe, at least 400 interphase cells with well-delineated fluorescent spots were examined. Nuclei were counterstained with 4′,6′-diamidino-2-phenylindole (DAPI) and antifade reagent, and signals were visualized using an Olympus BX51 microscope (Olympus Italia, Milan, Italy).

#### Fluorescence-activated cell sorting (FACS) analysis

The count of CD19^{+}CD5^{+} cells is a standard assay to define the representation of CLL cells in a diagnostic or research sample as measured by Fluorescence-activated cell sorting (FACS) analysis. A FACSCalibur flow cytometer (Becton–Dickinson) was utilized for the analysis. Expression of CD5 and CD19 was analyzed by combining Peridinin-Chlorophyll-Protein–Cyanine-5.5 (PerCP–Cy5.5)-conjugated anti-CD19 mAbs and fluorescein isothiocyanate-conjugated anti-CD5 mAbs. In order to estimate the proportion of cells co-expressing CD19 and CD5, for each sample events were acquired by gating on low forward and side scatter (FSC/SSC) CD19^{+} cells, which were further divided into CD5^{−} and CD5^{+} subsets. Irrelevant isotype-matched antibodies (Becton–Dickinson) were used to determine background fluorescence. FACS data of all samples are listed in Supplementary file 2.

#### Adjustment of mutation frequency

To consider tumor content, we performed CD19/CD5 FACS analysis to quantify the fraction of tumor cells. To unify FISH and NGS data, and also to adjust MAF of mutations in genes with copy number abnormalities, we introduced mutation cell frequency (MCF). As shown in Figure 1—figure supplement 1A, MCF represents the fraction of cancer cells containing particular alterations. Unlike MAF, copy numbers or tumor content do not affect MCF of a point mutation. For instance, if tumor purity is 70% and MAF of a heterogeneous variance in a diploid region is 0.2, then MCF is approximated by 0.2 × 2/0.7 ≈ 0.57. To infer MCF of different types of lesions, we applied the following strategy:

i. MCF of a copy number abnormality (i.e., T12, del11q, del13q_x1, del13q_x2, del17p, delBIRC3, delRB1_x1, delRB1_x2) was calculated by FISH analysis divided by the fraction of CD19/CD5 cells based on FACS, correcting for the tumor purity.

ii. Genes whose copy numbers were not frequently changed, such as

*MYD88*and*SF3B1*, were considered as diploid. So the MCF should be close to twofolds of MAF. However, MAF of some mutations might exceed 0.5, because of the random noise introduced in PCR. Here, MCF could be simply defined as a Piecewise function:

But this function was not smooth and it was not able to distinguish MAF = 0.5 and MAF = 0.9. To smooth this function, a Hill function was introduced, assuming$MCF=\raisebox{1ex}{$K{\left(MAF\right)}^{n}$}\!\left/ \!\raisebox{-1ex}{$\left[1+K{\left(MAF\right)}^{n}\right]$}\right.,$where *K* and *n* are parameters. To find optimal parameters, such that *MCF* ≈ 2 × *MAF*, we optimize the following objective function:${\text{min}}_{K,n}F\left(K,n\right)=\left[\underset{0}{\overset{0.5}{{\displaystyle \int}}}{\left(\raisebox{1ex}{$K{\left(MAF\right)}^{n}$}\!\left/ \!\raisebox{-1ex}{${\left[1+K{\left(MAF\right)}^{n}\right]}^{-2\times MAF}$}\right.\right)}^{2}dMAF\right].$

A grid-search method was applied to exhaust potential values of *K* and *n*. Particularly, for *K* ∈ {2^{0},⋯, 2^{9}} and *n* ∈ {1,⋯, 10}, *F* was calculated by the numerical integration of the above equation. Figure 1—figure supplement 1B shows that the optimal solution is $\widehat{K}=16$ and $\widehat{n}=2$. So for mutations of *MYD88* and *SF3B1,*$MCF=\raisebox{1ex}{$16{\left(MAF\right)}^{2}$}\!\left/ \!\raisebox{-1ex}{$\left[1+16{\left(MAF\right)}^{2}\right]$}\right..$

This Hill function, which was used to smooth the piecewise function, showed a more smooth correlation between MAF and MCF (Figure 1—figure supplement 1C). Additionally, by comparing FACS measurement of tumor purity, this Hill function form of MCF is more robust than either MAF without adjustment or simplistic piecewise function in assessing the fraction of cancer nuclei (Figure 1—figure supplement 1E).

iii. Genes that are frequently deleted in CLL, such as

*TP53*and*BIRC3*, were analyzed in the model without considering copy-neutral LOH or homozygous deletions (Figure 1—figure supplement 1D). We assume four different cell types: wild type (upper left), heterozygous deletion (upper right), mutation (lower left), and both mutation and deletion (lower right). The cell fractions were represented by*x*_{i}(*i*= 1,⋯, 4), respectively. So

and then MCF of mutation depends on MAF of mutation, as well as MCF of the deletion:${MCF}_{mut}={MAF}_{mut}(2-{MCF}_{del}).$

To smooth it, we optimized${\text{min}}_{K,n}F\left(K,n;{MCF}_{del}\right)=\left[\underset{0}{\overset{0.5}{{\displaystyle \int}}}{\left(\raisebox{1ex}{$K{\left(MAF\right)}^{n}$}\!\left/ \!\raisebox{-1ex}{${\left[1+K{\left(MAF\right)}^{n}\right]}^{-2\times MAF+MAF\times {MCF}_{del}}$}\right.\right)}^{2}dMAF\right].$

Here if we use $\widehat{K}$ and $\widehat{n}$ to represent the optimal solution given *MCF*_{del}, for mutations of *TP53* and *BIRC3*,$MCF=\raisebox{1ex}{$\widehat{K}{\left(MAF\right)}^{\widehat{n}}$}\!\left/ \!\raisebox{-1ex}{$\left[1+\widehat{K}{\left(MAF\right)}^{\widehat{n}}\right]$}\right..$

The MCF values were then divided by the fraction of CD19^{+}CD5^{+} cells based on FACS.

#### Simulation of linear and branching cancer models

To generate the artificial data for TEDG, we simulate cancer clonal evolution as a one-step Markov process, that is, the transition probability of mutation profile $P\left({\psi}_{{t}_{k}}|{\psi}_{{t}_{k-1}},\cdots ,{\psi}_{{t}_{1}}\right)=P\left({\psi}_{{t}_{k}}|{\psi}_{{t}_{k-1}}\right)$, where *ψ*_{tk} is the observed mutation profile at time *t*_{k}, which depended on the mutation profile at time *t*_{k−1}. To define the transition probability, we focused on two simple models of four mutations *x*_{1},*x*_{2},*x*_{3},*x*_{4}.

In the linear evolution model, we assumed that all mutations were mutated in a linear order shown in the left panel of Figure 2A. All possibilities of mutation status are$\left\{{\pi}_{0}=\varphi ,\hspace{0.17em}{\pi}_{1}=\left({x}_{1}\right),\hspace{0.17em}{\pi}_{2}=\left({x}_{1},{x}_{2}\right),\hspace{0.17em}{\pi}_{3}=\left({x}_{1},{x}_{2},{x}_{3}\right),{\pi}_{4}=\left({x}_{1},{x}_{2},{x}_{3},{x}_{4}\right)\hspace{0.17em}\right\}.$

According to Nordling's multi-mutation model (Nordling, 1953), the transition probability was defined as$P\left({\psi}_{{t}_{k}}={\pi}_{j}\text{|}{\psi}_{{t}_{k-1}}={\pi}_{i}\right)=\{\begin{array}{ll}{\left[1-{e}^{-f\xb7\left({t}_{k}-{t}_{k-1}\right)}\right]}^{j-i}-{\left[1-{e}^{-f\xb7\left({t}_{k}-{t}_{k-1}\right)}\right]}^{j-i+1}& ,j>i;\\ {e}^{-f\xb7\left({t}_{k}-{t}_{k-1}\right)}& ,j=i;\\ 0& ,j<i.\end{array},$where *f* represents fitness of the new mutation. The transition probability is shown in Figure 2B. The longitudinal data were then generated for each patient following the above model. To simplify the model, we fixed the number of time points (three) and the length of interval between time points (ten). We asked whether or not TEDG framework can reconstruct the order of mutations and how many patients are required. To answer this question, we applied a grid-search strategy to optimize parameter *β* and number of patients. Specifically, all possible combinations of parameter *β* from 0 to 1 (step by 0.1) and patient number from 1 to 100 (step by 1) are exhausted. For each combination, we randomly simulated the cancer patients and applied TEDG framework to reconstruct the order for 10 times. The frequency of reconstructing the exact order was defined as the accuracy of TEDG. Figure 2D shows an example with *β* equals 0.2 and the number of patients equals 15. Edge weights of ISN present the number of simulated patients with one mutation happening before another. The techniques of deconvolution and minimal spanning tree reconstruct the real structure by removing indirect interactions. Figure 2E summarizes the correlation between number of patients and optimal accuracy of TEDG.

Different from linear model, the branching model assumes that *x*_{3} and *x*_{4} are independently following the mutation of *x*_{2}. All possibilities of mutation status are$\left\{{\pi}_{0}=\varphi ,\hspace{0.17em}{\pi}_{1}=\left({x}_{1}\right),\hspace{0.17em}{\pi}_{2}=\left({x}_{1},{x}_{2}\right),\hspace{0.17em}{\pi}_{3}=\left({x}_{1},{x}_{2},{x}_{3}\right),{\pi}_{4}=\left({x}_{1},{x}_{2},{x}_{4}\right)\hspace{0.17em}\right\}.$

According to Nordling's multi-mutation model (Nordling, 1953), when *j* ≤ 2, the transition probability is the same as linear model:$P\left({\psi}_{{t}_{k}}={\pi}_{j}\text{|}{\psi}_{{t}_{k-1}}={\pi}_{i}\right)=\{\begin{array}{ll}{\left[1-{e}^{-f\xb7\left({t}_{k}-{t}_{k-1}\right)}\right]}^{j-i}-{\left[1-{e}^{-f\xb7\left({t}_{k}-{t}_{k-1}\right)}\right]}^{j-i+1}& ,j>i;\\ {e}^{-f\xb7\left({t}_{k}-{t}_{k-1}\right)}& ,j=i;\\ 0& ,j<i.\end{array}$

When *j* > 2, the transition probability is$P\left({\psi}_{{t}_{k}}={\pi}_{j}\text{|}{\psi}_{{t}_{k-1}}={\pi}_{i}\right)=\{\begin{array}{ll}{\left[1-{e}^{-f\xb7\left({t}_{k}-{t}_{k-1}\right)}\right]}^{3-i}& ,i<3;\\ 1& ,i=j;\\ 0& ,otherwise.\end{array}$

Figure 2F shows an example of branching model with *β* equals 0.2 and number of patient equals 15. Figure 2G summarized the correlation between number of patients and optimal accuracy of TEDG in branching model.

#### Network construction and analysis

To construct the sequential network of alterations for each tumor, we monitored the presence or absence of each genetic lesion in each sample. By collecting the clonal representation of each lesion from NGS and FISH analysis, we defined the status of each lesion with a cut-off of 5%. If the frequency of a genetic lesion is larger than 5%, we will name it as present; otherwise absent. If event A predates event B, we added a directed link between A and B. The ISN was constructed by pooling all sequential networks from different patients. To simplify ISN, we removed self-loops by subtracting the weight of weaker direction. To hierarchically layout the simplified ISN, we use yFiles Hierarchical Layout algorithm in Cytoscape 2.8.2, which well represents main direction or ‘flow’ in a directed network. With this method, nodes were placed in hierarchically arranged layers and the nodes within each layer are ordered in such a way that minimizes the number of edge crossings (Smoot et al., 2011). To statistically test whether the alterations are temporally ordered or randomly accumulated, we assume that in-degree (the number of incoming arrows) ${d}_{in}^{i}$ is equal to out-degree (the number of outgoing arrows) ${d}_{out}^{i}$ for each node *i*, and then $\sum}_{i}\frac{{\left({d}_{in}^{i}-{d}_{out}^{i}\right)}^{2}}{{d}_{out}^{i}$ follows a chi square distribution with degree of freedom *n − 1,* where *n* is the number nodes in the network. To test the statistical association of lesions in TEDG, we counted the number of samples carrying each pair of lesions and then calculated the p-value based on hypergeometric distribution to test whether the two genetic lesions are independent or not. To consider the effect of multi-hypothesizes, we corrected p-values with Bonferroni method and made a cut-off of 0.05. If two lesions were significantly co-mutated, a red link was added, while if they were significantly mutually exclusive, a blue edge was added.

Suppose *G*_{dir} is the adjacent matrix of all direct interactions/orders, the simplified ISN should be a summary of direct and indirect orders in the deconvolution formula *G*_{dir} = *G*_{obs}(*I* + *βG*_{obs})^{−1}, where, *G*_{obs} is the observed weighted adjacent matrix, *I* is the identify matrix, and *β* is a scaling factor between zero and one indicating the degree of deconvolution. The resulting matrix of deconvolution formula indicates the score of an edge to be a direct interaction (edge weights in the middle panel of Figure 2D,F). We then introduced a minimal spanning tree-based method to determine the final TEDG, in which we transformed the deconvoluted weights with a negative exponential function, and then calculated the minimal spanning tree with Prim's algorithm.

#### Fitting the evolutionary models

We developed Fit the Evolutionary Model (FEM) to properly fit the evolutionary model by systematically identifying clonal replacement in all sample pairs. FEM defines *Z*-scores for variant allele frequency by NGS and percentage of nuclei harboring abnormalities by FISH to represent normalized changes of the frequency of the lesions. For *i*th mutation in *j*th patient at the *k*th time point *t*_{k}${Z}_{i}^{j}\left({t}_{k}\right)=\frac{{f}_{i}^{j}\left({t}_{k}\right)-{f}_{i}^{j}\left({t}_{k-1}\right)}{{\sigma}_{seq}\sqrt{\frac{{f}_{i}^{j}\left({t}_{k}\right)+{f}_{i}^{j}\left({t}_{k-1}\right)}{2}\left(1-\frac{{f}_{i}^{j}\left({t}_{k}\right)+{f}_{i}^{j}\left({t}_{k-1}\right)}{2}\right)}},$where *f* is the mutation frequency of some particular gene mutation and *σ*_{seq} is its standard variation. Similarly, for the *i*th copy number change in the *j*th patient at the *k*th time point *t*_{k}${Z}_{i}^{j}\left({t}_{k}\right)=\frac{{f}_{i}^{j}\left({t}_{k}\right)-{f}_{i}^{j}\left({t}_{k-1}\right)}{{\sigma}_{FISH}\sqrt{\frac{{f}_{i}^{j}\left({t}_{k}\right)+{f}_{i}^{j}\left({t}_{k-1}\right)}{2}\left(1-\frac{{f}_{i}^{j}\left({t}_{k}\right)+{f}_{i}^{j}\left({t}_{k-1}\right)}{2}\right)}},$where *f* is the frequency nuclei harboring the abnormality of some particular gene mutation and *σ*_{FISH} is its standard variation.

The change of genetic lesion frequency is a synergic effect of treatment, tumor progression, and experimental noises, such as sequencing error or change of tumor purity. To obtain the variance caused by background noises, FEM robustly fit *σ*_{seq} and *σ*_{FISH} by eliminating data obviously affected by treatments or tumor progression. Based on this, we could calculate p-values of each genetic lesion in all sample pairs. To assess whether there was a significantly increased (or decreased) subclone in a given sample pair, we use Fisher's combinational test to combine p-values of all increased (or decreased) genetic lesions. The resulting p-values were separately defined as P_increase and P_decrease. Sample pairs that were significant in both P_increase and P_decrease were defined as sample pairs containing replacement events, which were further fit to the branching evolution model. All the others were compatible with the gradual linear model.

#### Growth rate and maximal mutation frequency slope (MMFS)

To investigate the fitness of genetic lesions, we defined growth rate and the maximal mutation frequency slope (MMFS). For the *i*th mutation in the *j*th patient at the *k*th time point *t*_{k}, mutation frequency slope is defined as ${s}_{i}^{j}\left({t}_{k}\right)=\frac{{z}_{i}^{j}\left({t}_{k}\right)}{{t}_{k}-{t}_{k-1}}$, where *z* is the *Z*-scores defined above. The growth rate of *i*th mutation in the *j*th patient at the *k*th time point *t*_{k} is defined as ${g}_{i}^{j}\left({t}_{k}\right)=max\left\{{s}_{i}^{j}\left({t}_{k}\right),0\right\}$. Furthermore, for the *j*th patient, MMFS is defined as ${s}^{j}={max}_{i}\left\{{max}_{k}\left[{s}_{i}^{j}\left({t}_{k}\right)\right]\right\}$.

#### Survival analysis

Overall survival (OS) was measured from date of initial presentation to date of death from any cause (event) or last follow-up (censoring). The cumulative probability of Richter syndrome transformation was measured from date of initial presentation to date of the biopsy documenting Richter syndrome transformation (event), death or last follow-up (censoring). Survival analysis was performed by the Kaplan–Meier method. The crude association between time-fixed exposure variables at diagnosis and survival was estimated by Cox proportional hazard regression. The analysis was performed with the Statistical Package for the Social Sciences (SPSS) software v.20.0 (Chicago, IL).

#### Candidate gene mutation screening

The mutation hotspots of *TP53* (exons 4–9, including splicing sites; RefSeq NM_000546.5), *NOTCH1* (exon 34, including splicing sites; RefSeq NM_017617.2), *SF3B1* (exons 14, 15, 16, 18, including splice sites; RefSeq NM_012433.2), *MYD88* (exons 3, 5, including splicing sites; RefSeq NM_002468.4), and *BIRC3* (exons 6–9, including splicing sites; RefSeq NM_001165.4) genes were analyzed by PCR amplification and DNA direct sequencing of high-molecular weight genomic DNA. Sequences for all annotated exons and flanking splice sites were retrieved from the UCSC Human Genome database using the corresponding mRNA accession number as a reference. PCR primers, located ∼50 bp upstream or downstream to target exon boundaries, were either derived from previously published studies or designed in the Primer 3 program (http://frodo.wi.mit.edu/primer3/) and filtered using UCSC in silico PCR to exclude pairs yielding more than a single product. All PCR primers and conditions are listed in Supplementary file 1. Purified amplicons were subjected to conventional DNA Sanger sequencing using the ABI PRISM 3100 Genetic Analyzer (Applied Biosystems) and compared to the corresponding germline sequences using the Mutation Surveyor Version 4.0.5 software package (SoftGenetics) after automated and/or manual curation. Of the evaluated sequences, 99% had a Phred score of 20 or more and 97% had a score of 30 or more. Candidate variants were confirmed from both strands on independent PCR products. The following databases were used to exclude known germline variants: Human dbSNP Database at NCBI (Build 136) (http://www.ncbi.nlm.nih.gov/snp); Ensembl Database (http://www.ensembl.org/index.html); The 1000 Genomes Project (http://www.1000genomes.org/); five single-genome projects available at the UCSC Genome Bioinformatics resource (http://genome.ucsc.edu/). Synonymous variants, previously reported germline polymorphisms and changes present in the matched normal DNA were removed from the analysis.

*IGHV*-*IGHD*-*IGHJ* rearrangement analysis

PCR amplification of *IGHV*-*IGHD*-*IGHJ* rearrangements was performed on high molecular weight genomic DNA using *IGHV* leader primers or consensus primers for the *IGHV* FR1 along with appropriate IGHJ genes, as previously described. PCR products were directly sequenced with the ABI PRISM BigDye Terminator v1.1 Ready Reaction Cycle Sequencing kit (Applied Biosystems) using the ABI PRISM 3100 Genetic Analyzer (Applied Biosystems). Sequences were analyzed using the IMGT databases and the IMGT/V-QUEST tool (version 3.2.17, Université Montpellier 2, CNRS, LIGM, Montpellier, France).

## References

- A mathematical framework to determine the temporal sequence of somatic genetic events in cancerProceedings of the National Academy of Sciences of USA, 107, 17604-17609, 2010
- Genomic aberrations and survival in chronic lymphocytic leukemiaThe New England Journal of Medicine, 343, 1910-1916, 2000
- Chronic lymphocytic leukemia/small lymphocytic lymphomaWorld health organization classification of tumours, pathology and genetics of tumours of haematopoietic and lymphoid tissues, IARC, Lyon, France, 180-182, 2008
- Monoclonal tumor origin is an underlying misconception of the RESIC approachProceedings of the National Academy of Sciences of USA, 108, E15, 2011

## Acknowledgements

We thank Dr Riccardo Dalla-Favera, Alexandra Jacunski, Dr Joseph Chan, and Dr Yong Wang for their helpful comments. LP is on leave from the University of Perugia Medical School.

## Decision letter

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for sending your work entitled “Tumor Evolutionary Directed Graphs and the History of Chronic Lymphocytic Leukemia” for consideration at *eLife*. Your article has been favorably evaluated by Chris Ponting (Senior editor), a Reviewing editor, and 2 reviewers. The Reviewing editor has assembled the following comments to help you prepare a revised submission.

Overall, the reviewers found the work interesting but raised a number of concerns that need to be addressed regarding the analysis of the data as well as clarifications on the design and presentation of the work.

Major comments:

1) Relative timing of mutations: It is not entirely clear in the Materials and methods or main text exactly how the temporal order of mutations is inferred. At one point, mutations are described as 'present' or 'absent' on the basis of allele fraction >5% and compared across time points. Really, the only true method for inferring temporal ordering is through reconstruction of phylogenetic trees; this can be difficult from bulk sequencing data, but there are methods emerging to do this.

2) Copy-number adjusted allele fractions: On a similar note, it is not clear that the fraction of cells carrying a point mutation has been adjusted for the copy number at that locus. This is especially important for the *TP53* point mutations when the other allele has undergone LOH. Also, how are double hits in a given gene dealt with in the graph method?

3) Rate of increase: The argument that later mutations show more rapid increase than earlier mutations is difficult to sustain on the basis of these data. The initiating lesion, by definition, will be present in all the tumor cells, and therefore cannot increase at the same magnitude as later mutations. Given that it is a prerequisite for entry into the study that the patient be diagnosed with CLL, one would imagine that, at that diagnostic time point, all the earliest mutations might already be clonal, and therefore will not change over time.

4) A question for another type of inference approach is: How do we know ground truth? The only way to know is if there were a set of samples from which there was frequent sampling, and if there were within clinical defined groups. Robustness of the network is dependent on how good the data is going in. While they do mention that they selected 70 samples for TEDG for which there were at least two sequential samples, in all likelihood, 2 samples are inadequate for defining the hierarchical relationships. There needs to be some table that provides information on the distribution of sequential samples per patient, and what types of time intervals between sequential samples. While a figure about the samples is provided as Figure 3–figure supplement 1, this is very hard to digest. This is a very very heterogeneous dataset, with 10 of 70 proceeding to Richters, more than half with treatment. A question is whether or not there are enough samples, by the time that one breaks down the patients into discrete clinical groups, to understand the hierarchical relationships.

Related to this, the numbers of samples are very confusing. This is what the reviewer understood: it seems that they start out with 202 patients, but only 70 had known highly recurrent drivers. 10 of the 70 had Richters Transformation, and how many have been excluded? By the time we get to Figure 4, there are 32 samples, but there is no discussion how this number came about. In general, in what percentage of sample was a polyclonal tumor population or a subclonal driver identified at the time of diagnosis?

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled “Tumor Evolutionary Directed Graphs and the History of Chronic Lymphocytic Leukemia” for further consideration at *eLife*. Your revised article has been favorably evaluated by Chris Ponting (Senior editor), a member of the Board of Reviewing Editors, and the original two reviewers. The manuscript has been improved but there are two remaining major issues that will need to be addressed before acceptance, as outlined below. We have copied the exact wording of the reviewer to facilitate the response.

Major issues:

1) A major revision has been the addition of the concept of MCF. This was done to address the concern of adjusting for local copy number changes that could affect estimation of the clone size with a gene alteration. Unfortunately, the approach taken by the authors to calculate MCF seems not well defined. This reviewer is concerned that sampling error in VAF values might fluctuate above or below 0.5 in ways that have a large influence on the inferred MCF. No convincing data is presented by the authors to justify their novel and elaborate seeming approach. Several approaches for correcting mutation VAF values for tumor purity and somatic copy number changes have been described (e.g. Carter et al., 2012, Nat Biotech; Landau et al., 2013, Cell; Fischer et al., 2014, Cell Reports). Notably, these methods rely on analysis of germline heterozygous SNPs at the mutant locus to infer CN-LOH, which is a much more robust approach. Could the authors simply use one of these established (and principled) approaches? This would leverage their exome data more fully.

2) They rely on FISH analysis for estimations of trisomy 12 and deletions of various chromosomal regions. This is not accurate, as it is not possible to correct for purity when using FISH data. In addition, copy neutral LOH cannot be inferred from FISH.