1. Genomics and Evolutionary Biology
  2. Human Biology and Medicine
Download icon

Tumor evolutionary directed graphs and the history of chronic lymphocytic leukemia

  1. Jiguang Wang Is a corresponding author
  2. Hossein Khiabanian
  3. Davide Rossi
  4. Giulia Fabbri
  5. Valter Gattei
  6. Francesco Forconi
  7. Luca Laurenti
  8. Roberto Marasca
  9. Giovanni Del Poeta
  10. Robin Foà
  11. Laura Pasqualucci
  12. Gianluca Gaidano
  13. Raul Rabadan Is a corresponding author
  1. Columbia University, United States
  2. Amedeo Avogadro University of Eastern Piedmont, Italy
  3. Centro di Riferimento Oncologico, Italy
  4. University of Southampton, United Kingdom
  5. Southampton University Hospital Trust, United Kingdom
  6. Catholic University of the Sacred Heart, Italy
  7. University of Modena and Reggio Emilia, Italy
  8. Tor Vergata University, Italy
  9. Sapienza University, Italy
Research Article
Cited
6
Views
2,635
Comments
0
Cite as: eLife 2014;3:e02869 doi: 10.7554/eLife.02869

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.

https://doi.org/10.7554/eLife.02869.001

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.

https://doi.org/10.7554/eLife.02869.002

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).

Figure 1 with 2 supplements see all
Tumor Evolutionary Directed Graph (TEDG) framework.

(A) A toy example of clonal evolution of one patient. The evolutionary history of four alterations A, B, C, and D is shown in the left panel. We then sample different time points and analyze genomic data. Specifically, for each patient, tagged-amplicon library next generation sequencing (NGS) and fluorescence in situ hybridization (FISH) analyses are carried out at different time points to evaluate the presence and quantify the clonal abundance of possible driver genetic lesions. Then we use mutation cell frequency (MCF) to adjust and unify the data (middle panel). Based on this longitudinal data, we build sequential network of one patient (right panel). CNA: copy number abnormalities. (B) Sequential networks derived from different patients (left panel) are further pooled to generate Integrated Sequential Network (ISN), which is a cross-sectional integration of longitudinal data (middle panel). We then infer TEDG by removing the indirect interactions with network deconvolution and simplification algorithms (‘Materials and methods’). To construct TEDG, we calculate a minimal spanning tree-based on the deconvolution scores (right panel).

https://doi.org/10.7554/eLife.02869.003

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 Gdir 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 Gdir = Gobs(I + βGobs)−1, where, Gobs 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’).

The calibration of TEDG framework on the simulation of two basic evolutionary models.

(A) The representation of linear model and branching model showing the sequential orders of four alterations. (B) The probability of observing different mutation profiles by Nordling's multi-mutation model. Specifically, if patient i harbored k0 mutations at time point t1, the probability to observe k more mutations at time point t2 is (1ef·(t2t1))k, where f represents the fitness of the new mutation. (C) The selection of parameter β. The color of heat map represents the accuracy of TEDG analysis, which is defined by the probability of TEDG reconstructing the input model. (D) An example showing the deconvolution algorithm on linear model. (F) An example showing the deconvolution algorithm on branching model. The left panel represents the ISNs of simulations of 15 patients. The weight of the edges represents the number of patients supporting the corresponding edges. The figures in the middle show the results of network deconvolution, where the numbers on the edge indicate deconvoluted weights. (E and G) The accuracy of TEDG framework on linear model and branching model at β = 0.2.

https://doi.org/10.7554/eLife.02869.006

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).

Figure 3 with 2 supplements see all
Evolutionary network analysis of CLL genetic lesions.

(A) Network representing the sequential order of genomic alterations in CLL. The nodes in the network represent genetic alterations and the oriented edges (arrows) represent sequential events in different patients where alterations in one gene predate alterations in other genes. The size of the nodes represents the recurrence of alterations in our cohort. The thickness and the color codes of the edges represent the number of patients showing a specific connection between nodes. (B) TEDG of CLL, which is the deconvolution of ISN by removing indirect interactions, representing the optimal tree to explain observed orders in CLL patients. (C) The order of CLL alterations. We calculate both the sequential in-degree (number of arrows to a node) and out-degree (number of arrows from the node) of each genetic lesion in ISN and use the binomial test to assess the significance by assuming that the number of in-degree and out-degree is randomly distributed. Our null assumption is that if there is no preferential time ordering, there should be the same number of alterations in gene A occurring before alterations in gene B and vice versa, up to statistical fluctuations. Deviations from the null assumption indicate a preferential order in the development of the alterations. All events are sorted by fold change between out-degree and in-degree and all significant early or late events are labeled by *(p-value < 0.05) or **(p-value < 0.01).

https://doi.org/10.7554/eLife.02869.007

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.

Figure 4 with 1 supplement see all
TEDG analysis of specific genetic lesions.

(A) Network representing the sequential order of specific genomic alterations. Two alterations are connected by a directed edge if they are observed to successively appear in at least one patient. The patient ID labels the corresponding edges, which are further colored by the number of patients. Nodes are colored by the fold change between out-degree and in-degree, indicating temporal order of corresponding alterations. The size of the nodes represents the recurrence of the alteration. (B) Tumor Evolutionary Directed Graph of 32 specific mutations. Some nodes and arrows in this figure are manually merged or adjusted for the purpose of illustration.

https://doi.org/10.7554/eLife.02869.010

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).

Fitting the evolution models.

(A) Distribution of the clonal evolution pattern in the 60 non-Richter cases. Three of the 60 cases show replacement during tumor progression. (B) Scatter plot showing p-values of observing increased and decreased subclones in the 80 samples of the 60 multi-time point patients. Samples with evidence of clonal replacement are located in the right-up corner (highlighted by yellow circle). (C) Number of patients following linear vs branched evolution pattern according to SF3B1 mutational emergence and previous treatment (p-value 0.0016 by Fisher's exact test).

https://doi.org/10.7554/eLife.02869.012
Table 1

Patients with branching evolution*

https://doi.org/10.7554/eLife.02869.013
Patient IDSampling timeTreatmentIncreased alterationsDecreased alterationsDetectable alterations
#12001NoneSF3B1 (K666E), SF3B1 (K700E), del13q
2005FCSF3B1 (K666E)NoneSF3B1 (K666E), SF3B1 (K700E), del13q
2008FCR/CAM/RBENSF3B1 (K700E), del11qNoneSF3B1 (K666E), SF3B1 (K700E), del13q, del11q
2011BENCAMSF3B1 (K700E), del11q, del13q, delBIRC3SF3B1 (K666E)SF3B1 (K700E), del13q, del11q, delBIRC3
#142004RF+12, del13q
2007FC/CAMNoneNone+12, del13q
2010BENCAMTP53 (R248Q), del17p, SF3B1 (K666E)None+12, del13q, TP53 (R248Q), del17p, SF3B1 (K666E)
2012BENDOFATP53 (R248Q), del17p, SF3B1 (K666E), T12del13q+12, del13q, TP53 (R248Q), del17p, SF3B1 (K666E)
#332002NoneNOTCH1 (P2514-)
2004FCRSF3B1 (K700E)NoneNOTCH1 (P2514-), SF3B1 (K700E)
2009FCRSF3B1 (K700E)NOTCH1 (P2514-)SF3B1 (K700E)
  1. *

    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.

Figure 6 with 3 supplements see all
Association network of CLL lesions in a larger cohort.

(A) The mutation status of 1054 CLL patients, which is informative in 11 most common genetic lesions in this leukemia. (B) The association analysis of CLL genetic lesions in TEDG. Lesion pairs are connected if they are significantly co-mutated (red) or mutually exclusive (blue). (C) Box plot of growth rate of all genetic lesions.

https://doi.org/10.7554/eLife.02869.014
Table 2

Clinical features at presentation of the CLL cohorts*

https://doi.org/10.7554/eLife.02869.018
TEDG analysisSA analysis
NTotal%NTotal%
IGHV homology >98%376953.6532138038.6
del13q14357050.0682140348.6
+12197027.1205140314.6
del11q22-q2387011.413614039.7
del17p13207028.611414038.1
TP53 mutation227031.411414018.1
NOTCH1 mutation227031.4150140310.7
SF3B1 mutation177024.310014037.1
MYD88 mutation6708.65214033.7
BIRC3 mutation107014.34114032.9
BIRC3 deletion6708.67914035.6
  1. *

    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).

Table 3

Patients showing high maximal mutation frequency slope (MMFS)

https://doi.org/10.7554/eLife.02869.019
PatientMMFS*Fast growing mutationsMCF1 (%)MCF2 (%)dT (moths)TreatmentRichter syndrome transformation
#52117.7TP53 (P152+)2.111.4<1NoneYes
#689.1TP53 (R273C)0.0100.011NoneNo
#516.7SF3B1 (I704F)12.268.95RCVPYes
#376.4TP53 (R248Q)3.890.010RDHAOXYes
#575.9NOTCH1 (P2514-)0.894.212NoneYes
#474.6del13q37.6100.013NoneNo
#143.9del17p18.990.016ANo
#423.4TP53 (N239T)0.0100.016RDHAOXYes
#43.3del17p53.6100.015CLB-ONo
#543.0NOTCH1 (P2415-)100.0100.022FCONo
#222.9BIRC3 deletion0.093.929FCRNo
#202.8del13q0.0100.036CLBNo
#632.5BIRC3 (M388V)0.086.024FCRYes
#62.3del17p0.092.134CLBNo
#382.3NOTCH1 (P2514-)41.742.94CVPYes
#132.2TP53 (G136H)0.0100.045RDHAOXNo
#12.1del11q11.8100.041FCR/A/BRNo
  1. *

    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.

  2. 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:

MCF={2×MAF,MAF<0.51,MAF0.5.

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=K(MAF)n[1+K(MAF)n],

where K and n are parameters. To find optimal parameters, such that MCF ≈ 2 × MAF, we optimize the following objective function:

minK,nF(K,n)=[00.5(K(MAF)n[1+K(MAF)n]2×MAF)2dMAF].

A grid-search method was applied to exhaust potential values of K and n. Particularly, for K ∈ {20,⋯, 29} 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 K^=16 and n^=2. So for mutations of MYD88 and SF3B1,

MCF=16(MAF)2[1+16(MAF)2].

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 xi(i = 1,⋯, 4), respectively. So

MAFmut=x3+x42x1+x2+2x3+x4=x3+x42(x2+x4)=MCFmut2MCFdel

and then MCF of mutation depends on MAF of mutation, as well as MCF of the deletion:

MCFmut=MAFmut(2MCFdel).

To smooth it, we optimized

minK,nF(K,n;MCFdel)=[00.5(K(MAF)n[1+K(MAF)n]2×MAF+MAF×MCFdel)2dMAF].

Here if we use K^ and n^ to represent the optimal solution given MCFdel, for mutations of TP53 and BIRC3,

MCF=K^(MAF)n^[1+K^(MAF)n^].

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(ψtk|ψtk1,,ψt1)=P(ψtk|ψtk1), where ψtk is the observed mutation profile at time tk, which depended on the mutation profile at time tk−1. To define the transition probability, we focused on two simple models of four mutations x1,x2,x3,x4.

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

{π0=ϕ,π1=(x1),π2=(x1,x2),π3=(x1,x2,x3),π4=(x1,x2,x3,x4)}.

According to Nordling's multi-mutation model (Nordling, 1953), the transition probability was defined as

P(ψtk=πj|ψtk1=πi)={[1ef·(tktk1)]ji[1ef·(tktk1)]ji+1,j>i;ef·(tktk1),j=i;0,j<i.,

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 x3 and x4 are independently following the mutation of x2. All possibilities of mutation status are

{π0=ϕ,π1=(x1),π2=(x1,x2),π3=(x1,x2,x3),π4=(x1,x2,x4)}.

According to Nordling's multi-mutation model (Nordling, 1953), when j ≤ 2, the transition probability is the same as linear model:

P(ψtk=πj|ψtk1=πi)={[1ef·(tktk1)]ji[1ef·(tktk1)]ji+1,j>i;ef·(tktk1),j=i;0,j<i.

When j > 2, the transition probability is

P(ψtk=πj|ψtk1=πi)={[1ef·(tktk1)]3i,i<3;1,i=j;0,otherwise.

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) dini is equal to out-degree (the number of outgoing arrows) douti for each node i, and then i(dinidouti)2douti 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 Gdir 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 Gdir = Gobs(I + βGobs)−1, where, Gobs 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 ith mutation in jth patient at the kth time point tk

Zij(tk)=fij(tk)fij(tk1)σseqfij(tk)+fij(tk1)2(1fij(tk)+fij(tk1)2),

where f is the mutation frequency of some particular gene mutation and σseq is its standard variation. Similarly, for the ith copy number change in the jth patient at the kth time point tk

Zij(tk)=fij(tk)fij(tk1)σFISHfij(tk)+fij(tk1)2(1fij(tk)+fij(tk1)2),

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 ith mutation in the jth patient at the kth time point tk, mutation frequency slope is defined as sij(tk)=zij(tk)tktk1, where z is the Z-scores defined above. The growth rate of ith mutation in the jth patient at the kth time point tk is defined as gij(tk)=max{sij(tk),0}. Furthermore, for the jth patient, MMFS is defined as sj=maxi{maxk[sij(tk)]}.

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

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
    World health organization classification of tumours, pathology and genetics of tumours of haematopoietic and lymphoid tissues
    1. HK Müller-Hermelink
    2. E Montserrat
    3. D Catovsky
    4. E Campo
    5. NL Harris
    6. H Stein
    (2008)
    180–182, World health organization classification of tumours, pathology and genetics of tumours of haematopoietic and lymphoid tissues, Lyon, France, IARC.
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
    Monoclonal tumor origin is an underlying misconception of the RESIC approach
    1. BL Parsons
    (2011)
    E15, Proceedings of the National Academy of Sciences of USA, 108, author reply E16, 10.1073/pnas.1017998108.
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54

Decision letter

  1. Emmanouil T Dermitzakis
    Reviewing Editor; University of Geneva Medical School, Switzerland

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.

https://doi.org/10.7554/eLife.02869.022

Author response

We found the comments extremely helpful, and according to the suggestions we have made several major revisions summarized as follows:

A) We introduce the concept of mutation cell frequency (MCF) to adjust mutation allele fractions using copy-number information. The calculation of MCF not only eliminates potential bias caused by “two-hit” alterations but also unify next generation sequencing (NGS) and fluorescence in situ hybridization (FISH) for data integration. After the adjustment, all our main results are consistent as shown in the following.

B) To address the major comment #1, we introduce a phylogenetic tree method to infer timing order of mutations. Our analysis shows that the phylogenetic tree method and TEDGs differ in capturing the effect of competition between different subclones. However, the main results inferred from both methods are consistent, indicating the robustness of the overall analysis (see detailed response below).

C) To address the concern that the growth rates of initial clonal events are biased, we only consider the increase rates of early subclonal events (mutation frequency <20% at diagnosis). In the new analysis this conclusion remains.

D) To clarify the confusion of the sample size in our study, we generate a new figure with a clear statement of how we select patients, and how many samples of each patient. To show how many patients are necessary to reveal the ground truth, we simulate the evolutionary and sampling process by an artificial example. The results indicate 30 patients will be enough to reconstruct a simple evolutionary graph.

All comments are addressed in detail as follows:

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? (Note: comment #2 is addressed first because other comments will depend on the adjustment of mutation data.)

Author response image 1
Definition of mutation cell frequency.

MAF: mutation allele frequency; MCF: mutation cell frequency. The black lines within the circles represent DNA copies, and the crosses represent mutations. The contingent table shows the difference between MAF and MCF.

This is a good suggestion. To consider potential bias caused by copy number abnormalities, we introduce mutation cell frequency (MCF) instead of Mutation Allele Frequency (MAF). As shown in Author response image 1, MCF represents the fraction of cells containing alterations. Different from MAF, MCF of a mutation is not affected by copy numbers. To infer MCF of different types of alterations, we apply the following strategy:

i) The percentages of cells harboring the copy number abnormality (i.e. +12, del11q, del13q_x1, del13q_x2, del17p, delBIRC3, delRB1_x1, delRB1_x2) are estimated by FISH analysis.

Author response image 2
Optimization of Hill function by grid search method.

z-axis indicates the objective function, x-axis and y-axis are parameters of Hill function.

ii) Genes that are not frequently deleted or amplified in CLL, such as NOTCH1, MYD88, and SF3B1 are considered as diploid. So the MCF should be close to two folds of MAF. However, MAF of some mutations may exceed 0.5, because of the random noise introduced in PCR process or because sometimes the diploid assumption does not hold. MCF could be defined in a simple Piecewise function:

MCF={2×MAF,MAF<0.51,MAF0.5.

But this function is not smooth and it is not able to distinguish MAF=0.5 and MAF=0.9. To smooth this function, a Hill function strategy is introduced.

MCF=K(MAF)n[1+K(MAF)n],

Assume where K and n are parameters of the Hill function. To find optimal K and n such that MCF2×MAF, we optimize the following objective function:

minK,nF(K,n)=[00.5(K(MAF)n[1+K(MAF)n]2×MAF)2dMAF].
Author response image 3
Optimal correlation between MAF and MCF.

A grid search method is applied to exhaust potential values of K and n. Particularly, for K{20,,29}, and n{1,,10}, F is calculated by the numerical integration of the above equation. Author response image 2 shows the optimal solution is K^=16 and n^=2. So for mutations of NOTCH1, MYD88, and SF3B1:

MCF=16(MAF)2[1+16(MAF)2].

This Hill function, which is used to smooth the piecewise function, shows a more reasonable correlation between MAF and MCF (Author response image 3).

iii) Genes that are frequently deleted in CLL, such as TP53 and BIRC3, are analyzed in the model without considering copy-neutral LOH and homozygous deletions (Author response image 4). We assume four different genotypes: wild type (upper left), heterozygous deletion (upper right), mutation (lower left), and both mutation and deletion (lower right). The cell fractions are represented by xi(i=1,,4), respectively. So

MAFmut=x3+x42x1+x2+2x3+x4=x3+x42(x2+x4)=MCFmut2MCFdel

and then MCF of mutation depends on MAF of mutation, as well as MCF of the deletion:

MCFmut=MAFmut(2MCFdel).
Author response image 4
Genotypes of two-hit model.

Black lines within the circles represent copies of DNA, and crosses on the lines indicates point mutations. χ1, χ2, χ3, and χ4 indicate proportion of the cells with corresponding genotype.

To smooth it, we optimize:

minK,nF(K,n;MCFdel)=[00.5(K(MAF)n[1+K(MAF)n]2×MAF+MAF×MCFdel)2dMAF].

Assume K^ and n^ are optimal solution given MCFdel, for mutations of TP53 and BIRC3,

MCF=K^(MAF)n^[1+K^(MAF)n^].

The above methods have been applied to adjust the data reported in this study before the TEDG analysis. New analysis based on adjusted data shows that: 1) all results of the TEDG analysis are totally identical to the previous version; 2) the growth rate of the late events is still significantly higher than that of early events (p-value=0.0259). This adjustment not only reduces the potential bias caused by copy number alterations, but also addresses another concern mentioned in Minor comment #4 about the integration of FISH data and NGS data.

To modify our manuscript:

A) We have updated all relevant data in this manuscript, updated the description of TEDG framework in Figure 1A, and added the following sentence to the legend of the Figure 1:

“Then we use mutation cell frequency (MCF) to adjust and unify the data (middle panel).”

B) We have added the following sentence in Results, “Tumor Evolutionary Directed Graphs”:

“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 abnormality (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) for quantification of genetic lesions (Materials and methods, Figure 1–figure supplement 1).”

C) We have put Author response images 1–4 in Figure 1–figure supplement 1, and added expanded the Materials and methods accordingly.

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.

Author response image 5
Relative timing of mutations of 70 patients.

Each column represents one patient with at least two time points. Magenta (MCF>5%, present) and blue bars (absent), which are ordered by time information, indicate the mutation status of the corresponding alteration. For one patient, if the present of alteration A is earlier than B, we claim A predates B.

In our work, the relative timing of mutations is inferred by longitudinal data with multiple time points. If the presence of mutation A (allele fraction >5%) predates that of mutation B, mutation A is considered earlier than mutation B (Author response image 5).

Author response image 6
Phylogenetic trees of CLL patients.

Twenty-one out of all CLL patients contain the change of mutation status during disease progression. The phylogenetic trees are constructed based on mutation status of the driver genes. Green balls are normal cells, while all the others are cancer cells with particular alterations.

To explain this method clearly in the manuscript, we put Author response image 5 as Figure 1–figure supplement 2, and add the following sentence in Results, “Tumor Evolutionary Directed Graphs”:

“Particularly, we include CNA represented in >5% of leukemic cells by FISH and mutations with MCF>5% (see examples of CLL patients in Figure 1– supplement 2). First, if a given genetic lesion is observed to be temporally earlier than another lesion, we connected them with a directed edge to represent their sequential order of development (Figure 1A).”

We agree that the construction of phylogenetic trees is a good option to infer temporal orders of gene mutations by assuming common ancestors are early events, but limited to the fact that only very few genes have been considered in this study, it is difficult to reconstruct an accurate phylogenetic tree of patients. To approximate the phylogenetic trees, we define the pairwise distances between samples by considering both the number of common mutations and the number of different mutations, and make use of the complete linkage distance method to construct the phylogenetic trees of wild type cell, tumor cell at diagnosis, and tumor cell at different stages of relapse (Author response image 6). Other methods, such as neighbor-joining and parsimony, produce similar results, as expected from the small number of branches and informative sites (driver mutations that vary within a tumor) and, as described below, branching assessment using bootstrap is generically non-significant.

Author response image 7
Change of mutation frequency of patient #1.

There are four samples for patient #1. Mutation frequencies of different alterations change according to the progression of the disease.

In most of cases, phylogenetic tree method indicates the same order of mutations as our longitudinal method used in TEDG. However, the order of gene mutations inferred by longitudinal method and phylogenetic tree method are not totally identical. For example, the phylogenetic tree of patient #1 (Author response images 6 and 7) indicates del13q and delRB1 are early events, and then two branches are independently developed with one harbors SF3B1(K700E) and del11q, and the other harbors SF3B1(K666E). The first branch then develops SF3B1(K666E) and delBIRC3. The phylogenetic tree method does not predict the order of delBIRC3 and the mutation of SF3B1(K666E), but according to Author response images 5 and 7, in patient #1, mutation of SF3B1(K666E) predates delBIRC3. The temporal order of delBIRC3 and mutation SF3B1(K666E) may not occur in the same subclone, but it indicates the competition between two subclones, suggesting subclone with delBIRC3 may have stronger fitness in this particular patient. Therefore, different from the phylogenetic tree method, the longitudinal analysis of TEDG presents the temporal orders of mutations, not only the results of the disease progression within a subclone but also consequences of competition among different clones.

To compare the results between the phylogenetic trees method and TEDG, we have applied both methods in our CLL data. The results show that all significantly early or late events with p-value<0.01 (**) reported by TEDG (before or after adjustment of MAF) are still significant (Author response image 8A). NOTCH1, reported significant early in TEDG, is very early (not significant due to limited number of samples). Furthermore, we rank all events based on fold change between indegree and outdegree by phylogenetic tree method and TEDG separately. Strikingly, the rank of all events based on those two methods is significantly related with Pearson correlation more than 0.9 (Author response image 8C). In the TEDG network analysis, the backbone of CLL evolution is slightly changed, but two main branches are the same, indicating two major types of CLL patients suffering T12 and del13q independently (Author response image 8). The actual topologies of the phylogeny inferred network and TEDG are very similar (six out of ten edges are in common and share directionality, p-value<0.0001 by Fisher’s exact test). It is interesting to observe that the differences are in the order of del17p and TP53 mutations and the indirect association of SF3B1.

Author response image 8
The comparison between TEDG (B) and phylogenetic tree method (A).

The two methods are compared in terms of the order of mutations and the TEDG networks. * indicates p-value<0.05, and ** indicates p-value<0.01. p-value in (C) is calculated by Fisher’s exact test.

We would like also to emphasize that phylogenetic methods usually require higher number of alterations than the ones considered in our manuscript. Standard statistical techniques to assess the branching for each tree rely on bootstrapping segregating sites. 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 described above 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.

In summary, we introduce a phylogenetic tree method to compare with our longitudinal method, and show that although the two methods are slightly different in capturing the effect of competition of different subclones, the main results inferred from both methods are almost the same, indicating the robustness of our overall analysis.

To modify our manuscript, we add a new paragraph in the Discussion section:

“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 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.”

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.

To address this concern, we only consider the growth rate of subclonal mutations. Particularly, clonal mutations with MCF greater than 20% are eliminated from the comparison. As shown in Author response image 9, late events are significantly faster than early events with p-value =0.0483 (Wilcoxon rank-sum test), indicating late events grow fast is not only because of the limitations on early events.

To clarify this point in the manuscript, we have added Author response image 9 as Figure 6–figure supplement 3, and the following sentence in Results, “Rate of Allele Frequency Change”:

“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, 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.0483 (Figure 6–figure supplement 3).

Author response image 9
Comparison of the growth rate of subclonal alterations between early and late events.

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 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?

We do not know the ground truth, and this is the main purpose of our and other methods to infer the order of alterations from longitudinal data. However, we have performed simulations where we start with the ground truth, a model of how mutations accumulate and a sampling strategy. In particular, we tested the TEDG method in an artificial example by separately simulating linear evolution and branching evolution of cancer as one-step Markov process. The simulation shows that TEDG’s accuracy in the linear model reaches 80% when N is 30 (Figure 2E). Similarly, the accuracy reaches 90% when N is 30 in branching model (Figure 2F-G, Materials and methods). Simulations allow estimating the power of our data and method to assess the order of alterations and the different evolutionary trajectories.

Author response image 10
Summary of longitudinal data in 70 patients.

(A) The 70 patients are selected from a big cohort of 1,403 CLL patients with no-bias screening. (B) 70 patients are ranked according to their minimal cell frequency at diagnosis. Patient with minimal cell frequency less than 20 are in red, the others are in green.

Regarding the number of samples considered in each analysis, we are sorry for the confusion of the presentation of sample size in Figure 3–figure supplement 1 and Figure 4. To clarify the samples used in our study, we present our data in Author response image 10. No biases are present in the case series we have utilized in the manuscript. Indeed, we have collected 14,03 newly diagnosed CLL patients between the year of 2001 and the year of 2012, during which 202 patients contain samples of multiple time pints. We systematically genotyped all the sequential samples collected at each progression (for progressive cases) and at the last follow-up (for non- progressive cases) from consecutive cases recruited in our database. Out of those 202 patients, 70 are 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 (Author response image 10A). Our study is mainly focusing on those 70 patients. The mutation information, as well as the number of sequential samples of each patient is shown in Author response image 5. In total, we have 50 patients with two time points, 16 with three time points, and 4 with four time points. In Figure 4, we did not say there are 32 samples. We used 70 CLL patients to analyze evolution network of 32 alterations (one gene can have more alterations by considering mutation sites).

To answer the question that “in what percentage of sample was a polyclonal tumor population or a subclonal driver identified at the time of diagnosis”, we have added the distribution of minimal cell frequency at diagnosis as shown in Author response image 10B. Half (35 out of 70) of the patients have at least one subclonal mutation with cell frequency less than 20% at diagnosis.

To modify the manuscript:

A) We have added the following sentence in Results, “Tumor Evolutionary Directed Graphs”:

“To test TEDG method and also to show how many patients are required to approximate the ground truth, we employ the TEDG method in an artificial example by both simulating linear evolution and branching evolution of cancer, where the longitudinal data is generated by one-step Markov process and Nordling’s multi-mutation model.”

B) We have included Author response image 10 to replace Figure 3–figure supplement 1, and added the following sentence in Results, “Tumor Evolutionary Directed Graph of CLL”:

“Half (35 out of 70) of the patients have at least one subclonal genetic lesion with cell frequency less than 20% at diagnosis (Figure 3–figure supplement 1).”

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

We have addressed the remaining concerns about the justification of defining Mutation Cell Fraction (MCF) and the application to FISH analysis. Following the suggestions: we have 1) further clarified the definition of MCF in manuscript, 2) described the variations of MAF (VAF) caused by sampling error based on dilution experiments, 3) included FACS analysis to consider tumor purity to calibrate MCF, 4) justified the Hill function form of MCF by comparing with more simplistic definitions, and 5) compared the calibrated MCF with previously defined Cancer Cell Fraction (CCF) by ABSOLUTE 1.0.6 (Landau et al., 2013).

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.

1) To clarify the definition of MCF, we have revised manuscript as follows:

Excerpt from Results, “Tumor Evolutionary Directed Graphs”:

“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).”

Excerpt from Materials and methods, “Adjustment of Mutation Frequency”:

“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.”

2) To estimate the variation of MAF (VAF) caused by sampling error, we have employed a bioinformatic approach to call mutations of low abundance out of the background error noise of deep sequencing data. Based on dilution experiments, we calibrated for systematic biases that lead to sequencing errors and derived the depth distribution of sequencing errors to be negative binomial. We determined that the sequencing depth of 1,000x is sufficient for a highly sensitive detection of mutations with allele frequency >1% out of the background error noise and variance of the estimate is small when allele frequency is larger than the detection limit. A detailed account of calibration of estimates using ultra-deep sequencing has been described in a related publication (Figure S3 of Rossi et al., 2014, referenced in the manuscript).

Author response image 11
The impact of sampling error on MCF.

Confidence interval of MCF is inferred based on the variation of MAF estimated by dilution experiments.

To show how the sampling error affects MCF, we calculate the 95% confidence interval of MCF based on variation of MAF described by Rossi et al., 2014. Particularly, the 95% confidence interval of MAF at 0.5%, 2.5%, 5%, and 25% are used to estimate upper and lower bound of MCF with the formula described in Materials and methods, “Adjustment of Mutation Allele Frequency”, and Figure 1–figure supplement 1. The length of 95% confidence interval of MCF are respectively 0.1%, 0.6%, 1.2%, and 2%, indicating that the sampling error slightly fluctuates the value of MCF; however, only 0.2% (1 out of 415) of the detected mutations, which is on the boundary of our cutoff, will be affected by this fluctuation, not changing the overall results (Author response image 11).

To modify the manuscript, we have added the following sentence at the end of the Materials and methods section, “Estimation of the Variant Allele Frequency”:

“The sequencing depth in this study is on average 1200x, which is sufficient for a highly sensitive detection of mutations with allele frequency >1% out of the background error noise (Rossi et al., 2014).”

3) Our estimation of MCF in the first revision of manuscript was based on two assumptions: 1) there are no copy neutral LOHs (cn-LOH) and 2) tumor purity is very high. The two assumptions may limit the application of this approach in cn-LOH regions or low-purity samples.

The reviewers have suggested considering neighboring SNPs to address the concern of cn-LOH. Unfortunately, we do not have whole-exome or SNP array data for most of our longitudinal cases. However, whole-exome sequencing data are available for 10 out of 70 patients, including (a) four patients with paired tumor and normal control (published recently in Messina et al., Blood, 2014 and discussed below) and (b) six patients with only tumor samples (published recently in Fabbri et al., 2013, JEM, and discussed below) in our cohort. To check whether the driver genes (i.e., TP53, SF3B1, BIRC3, MYD88, NOTCH1) are affected by cn-LOH in those ten patients, we have checked the neighboring SNPs in IGV viewer, and find none of our driver genes in the regions of cn-LOH. In a recent study, Pfeifer et al. have analyzed 70 CLL patients, and listed the regions of cn-LOH in Table 4 of their paper (Pfeifer et al, 2007, Blood). We compare all our diver genes with their results, and find none of our genes in their list. In summary, based on limited cases with heterozygosity information in the list of driver genes selected in these studies and in published results by an independent group in an extended cohort, we reason that cn-LOH is not a major effect in the targeted genes of our current study.

To consider tumor purity, we have included FACS analysis of CD19/CD5 expression in the revised manuscript. CD19+CD5+ markers represent a gold standard to define the representation of CLL cells in a diagnostic or research sample, as these tumors characteristically and invariably express these molecules on their surface, in contrast with normal peripheral B cells.

To modify the manuscript, we have added the following sentences in Materials and methods, “Adjustment of Mutation Frequency”:

“To consider tumor content we performed CD19/CD5 FACS analysis to quantify the fraction of tumor cells.”

“The MCF values were then divided by the fraction of CD19+CD5+ cells based on FACS.”

We have also added the following new section in Materials and methods:

“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 isothyocyanate-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.”

Author response image 12
Justification of MCF.

X-axis indicates fraction of CD19+CD5+ cells assessed by FACS analysis, and y-axis indicates maximal mutation fraction of all targeted driver genes of each sample calculated by different methods. One blue dot represents one sample, and contours indicate the density of dots. A suitable calculation of maximal driver mutation fraction will approximate but not exceed the fraction of cancer nuclei. The upper red line indicates CD19+CD5+ cell fraction, and the lower red line indicates a 20% lower interval of it. Apparently, tumor purities of 55 samples are properly assessed by the Hill function MCF, which is better than both MAF without adjustment (10 samples) and simple piecewise MCF (27 samples).

4) To further justify the Hill function form of MCF, we have used data from FACS analysis to compare it to more simplistic methods. Particularly, given one sample, the maximal mutation fraction of all drivers is separately calculated based on three approaches: MAF, piecewise MCF, and Hill function MCF. A suitable calculation of this score will approximate but not exceed the fraction of cancer nuclei estimated by FACS analysis of the fraction of CD19+CD5+ cells. The comparison shows that the Hill function form of MCF is more robust to assess the tumor purity (Author Response Image 12).

Author response image 13
Comparison of MCF and CCF.

The samples with WES available show that our MCF falls within the estimates of ABSOLUTE. ABSOLUTE underestimates tumor purity in sample #37 as compared to CD19/CD5 FACS and our MCF estimate.

To modify the manuscript, we have used Author response image 12 as Figure1–figure supplement 1E, and added the following sentence in Materials and methods, “Adjustment of Mutation Frequency”:

“Additionally, by comparing to 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).”

5) As mentioned by the reviewers, several approaches for correcting MAF have been described recently. ABSOLUTE (Carter et al., 2012, Nat Biotech) infers tumor purity and average cancer cell ploidy from the analysis of SNP array and somatic DNA alterations. ABSOLUTE 1.0.5 (Landau et al., 2013, Cell) and 1.0.6 were further released to include the calculation of cancer cell fraction (CCF) for each sSNV. More recently, cloneHD (Fischer et al., 2014 Cell Reports) was published to reconstruct the subclonal structure of a population from sequencing data. All these approaches have been proved to be robust in the analysis of cancer genome, but all methods require genome-wide information, such as SNP array, whole-exome or whole-genome sequencing, which are only available for a few samples in our analysis. We want to emphasize that compared to the other approaches this manuscript presents the analysis of a few targeted alterations in large number of longitudinal samples (70 longitudinal patients).

Author response image 14
ABSOLUTE report of case #37.

Several models fit the data with different purity estimates.

To compare our method to genome-wide methods, we have analyzed previously published whole-exome sequencing (WES) and SNP array data of four of our patients: #57, #4, #37, and #59 (Fabbri et al., 2013, JEM, and Messina et al., 2014, Blood. SNP array data for the other six patients are not available). We ran ABSOLUTE 1.0.6 with the parameters min.ploidy=1, and max.ploidy=3 (given that all samples in the research of Landau et al. (Landau et al., 2013, Cell) were estimated to have near-diploid DNA content). CCFs (cancer cell fraction as defined by ABSOLUTE) are calculated. By comparing MCF with the confidence interval of CCF, for mutations of 3 out of 4 cases, the results are consistent (Author response image 13). In patient #37, however, according to ABSOLUTE, MCF of a TP53 mutation (MAF=0.34) is close to 100% because the tumor content is estimated to be ∼40%, which is not consistent with our FACS analysis (CD19+CD5+ cell fraction is 95%). This indicates that in this discrepant case, ABSOLUTE underestimates the tumor purity (40%) compared to FACS (95%). We were extremely curious about this underestimation and upon further investigation we find that ABSOLUTE reports different models that fit the data, and in fact there are several solutions that coincide with our estimate of MCF and the FACS analysis (Author response image 14).

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.

To address this concern, we have included FACS analysis. We performed the analysis with corrected tumor content based on CD19/CD5 FACS analysis and all the major results are the same as expected from the high tumor content of these samples (average 83%). Specifically, we updated all MCFs dividing by the fraction of CD19+CD5+ cells. After this correction, MCFs of a few number alterations, which are not considered before, exceed the cutoff of 5%, including del13q (delRB1) at the first time point of patients #13 and #58, and NOTCH1 mutation at the second time point of patient #34. Accordingly, Figure1–figure supplement 2, Figure 3A were updated (edges in the network plot, indegree, outdegree, etc.). Importantly, this modification did not change TEDG (Figure 3B) or the order of mutations (Figure 3C), indicating our analysis is robust to this adjustment.

As described earlier in response to Major comment 1, we have found that cn-LOH is not a major effect in the targeted genes of our current study.

To modify the manuscript, we have updated this data to all analysis in our manuscript, including the MCFs in Table 3, and modifications in figures: change of some arrows in Figure 4, indegree/outdegree in Figure3–figure supplement 2 and Figure4–figure supplement 1, boxplots in Figure 6C and Figure 6–figure supplement 3, etc.

Also we have added the following sentence in Materials and methods, “Adjustment of Mutation Frequency”:

“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.”

Additionally, we have also added a supplementary table with detailed data of CD19/CD5 as Supplementary file 2.

https://doi.org/10.7554/eLife.02869.023

Article and author information

Author details

  1. Jiguang Wang

    1. Department of Biomedical Informatics, Columbia University, New York, United States
    2. Department of Systems Biology, Columbia University, New York, United States
    3. Center for Computational Biology and Bioinformatics, Columbia University, New York, United States
    Contribution
    JW, Performed statistical analysis, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Hossein Khiabanian and Davide Rossi
    For correspondence
    jw2983@columbia.edu
    Competing interests
    The authors declare that no competing interests exist.
  2. Hossein Khiabanian

    1. Department of Biomedical Informatics, Columbia University, New York, United States
    2. Department of Systems Biology, Columbia University, New York, United States
    3. Center for Computational Biology and Bioinformatics, Columbia University, New York, United States
    Contribution
    HK, Performed statistical analysis, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Jiguang Wang and Davide Rossi
    Competing interests
    The authors declare that no competing interests exist.
  3. Davide Rossi

    Division of Hematology, Department of Translational Medicine, Amedeo Avogadro University of Eastern Piedmont, Novara, Italy
    Contribution
    DR, Performed statistical analysis, Acquisition of data, Drafting or revising the article
    Contributed equally with
    Jiguang Wang and Hossein Khiabanian
    Competing interests
    The authors declare that no competing interests exist.
  4. Giulia Fabbri

    1. Institute for Cancer Genetics, Columbia University, New York, United States
    2. Department of Pathology and Cell Biology, Columbia University, New York, United States
    Contribution
    GF, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  5. Valter Gattei

    Clinical and Experimental Onco-Hematology Unit, Centro di Riferimento Oncologico, Aviano, Italy
    Contribution
    VG, Collected patients and samples, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  6. Francesco Forconi

    1. Cancer Sciences Unit, Cancer Research UK Centre, University of Southampton, Southampton, United Kingdom
    2. Haematology Department, Southampton University Hospital Trust, Southampton, United Kingdom
    Contribution
    FF, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  7. Luca Laurenti

    Institute of Hematology, Catholic University of the Sacred Heart, Rome, Italy
    Contribution
    LL, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  8. Roberto Marasca

    Division of Hematology, Department of Oncology and Hematology, University of Modena and Reggio Emilia, Modena, Italy
    Contribution
    RM, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  9. Giovanni Del Poeta

    Department of Hematology, Tor Vergata University, Rome, Italy
    Contribution
    GDP, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  10. Robin Foà

    Department of Cellular Biotechnologies and Hematology, Sapienza University, Rome, Italy
    Contribution
    RF, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  11. Laura Pasqualucci

    1. Institute for Cancer Genetics, Columbia University, New York, United States
    2. Department of Pathology and Cell Biology, Columbia University, New York, United States
    Contribution
    LP, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  12. Gianluca Gaidano

    Division of Hematology, Department of Translational Medicine, Amedeo Avogadro University of Eastern Piedmont, Novara, Italy
    Contribution
    GG, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  13. Raul Rabadan

    1. Department of Biomedical Informatics, Columbia University, New York, United States
    2. Department of Systems Biology, Columbia University, New York, United States
    3. Center for Computational Biology and Bioinformatics, Columbia University, New York, United States
    Contribution
    RR, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    rr2579@cumc.columbia.edu
    Competing interests
    The authors declare that no competing interests exist.

Funding

National Center for Advancing Translational Sciences (UL1 TR000040)

  • Jiguang Wang

National Institutes of Health (1 U54 CA121852-05)

  • Hossein Khiabanian
  • Raul Rabadan

National Institutes of Health (1R01CA185486-01)

  • Jiguang Wang
  • Raul Rabadan

National Institutes of Health (1R01CA179044-01A1)

  • Jiguang Wang
  • Raul Rabadan

Associazione Italiana per la Ricerca sul Cancro (Special Program Molecular Clinical Oncology, 5x1000, No. 10007)

  • Robin Foà
  • Gianluca Gaidano

Associazione Italiana per la Ricerca sul Cancro (My First AIRC Grant 2012)

  • Davide Rossi

Ministero dell'Istruzione, dell'Università e della Ricerca (Futuro in Ricerca 2002 and 2012)

  • Davide Rossi

Istituto Nazionale di Fisica Nucleare (PRIN 2008)

  • Gianluca Gaidano

Istituto Nazionale di Fisica Nucleare (PRIN 2009)

  • Davide Rossi

Ministero della Salute (2008 and 2010)

  • Davide Rossi

Alexander and Margaret Stewart Trust

  • Raul Rabadan

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

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.

Ethics

Human subjects: The study was approved by the institutional ethical committee of the Azienda Ospedaliero-Universiataria Maggiore della Carita 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.

Reviewing Editor

  1. Emmanouil T Dermitzakis, Reviewing Editor, University of Geneva Medical School, Switzerland

Publication history

  1. Received: March 24, 2014
  2. Accepted: December 10, 2014
  3. Accepted Manuscript published: December 11, 2014 (version 1)
  4. Version of Record published: January 28, 2015 (version 2)

Copyright

© 2014, Wang et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,635
    Page views
  • 377
    Downloads
  • 6
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Neuroscience
    Hiroaki Koseki et al.
    Research Article Updated