Marc J. Williams

Identification of neutral tumor evolution across cancer types

Despite extraordinary efforts to profile cancer genomes, interpreting the vast amount of genomic data in the light of cancer evolution remains challenging. Here we demonstrate that neutral tumor evolution results in a power-law distribution of the mutant allele frequencies reported by next-generation sequencing of tumor bulk samples. We find that the neutral power law fits with high precision 323 of 904 cancers from 14 types and from different cohorts. In malignancies identified as evolving neutrally, all clonal selection seemingly occurred before the onset of cancer growth and not in later-arising subclones, resulting in numerous passenger mutations that are responsible for intratumoral heterogeneity. Reanalyzing cancer sequencing data within the neutral framework allowed the measurement, in each patient, of both the in vivo mutation rate and the order and timing of mutations. This result provides a new way to interpret existing cancer genomic data and to discriminate between functional and non-functional intratumoral heterogeneity.

Quantification of subclonal selection in cancer from bulk sequencing data

Subclonal architectures are prevalent across cancer types. However, the temporal evolutionary dynamics that produce tumor subclones remain unknown. Here we measure clone dynamics in human cancers by using computational modeling of subclonal selection and theoretical population genetics applied to high-throughput sequencing data. Our method determined the detectable subclonal architecture of tumor samples and simultaneously measured the selective advantage and time of appearance of each subclone. We demonstrate the accuracy of our approach and the extent to which evolutionary dynamics are recorded in the genome. Application of our method to high-depth sequencing data from breast, gastric, blood, colon and lung cancer samples, as well as metastatic deposits, showed that detectable subclones under selection, when present, consistently emerged early during tumor growth and had a large fitness advantage (>20%). Our quantitative framework provides new insight into the evolutionary trajectories of human cancers and facilitates predictive measurements in individual tumors from widely available sequencing data.This analysis uses computational modeling of clonal selection to measure evolutionary dynamics in primary human cancers. The method employs high-throughput sequencing data and simultaneously measures the selective advantage and time of appearance of subclones.

Somatic POLE exonuclease domain mutations are early events in sporadic endometrial and colorectal carcinogenesis, determining driver mutational landscape, clonal neoantigen burden and immune response.

Genomic instability, which is a hallmark of cancer, is generally thought to occur in the middle to late stages of tumourigenesis, following the acquisition of permissive molecular aberrations such as TP53 mutation or whole genome doubling. Tumours with somatic POLE exonuclease domain mutations are notable for their extreme genomic instability (their mutation burden is among the highest in human cancer), distinct mutational signature, lymphocytic infiltrate, and excellent prognosis. To what extent these characteristics are determined by the timing of POLE mutations in oncogenesis is unknown. Here, we have shown that pathogenic POLE mutations are detectable in non‐malignant precursors of endometrial and colorectal cancer. Using genome and exome sequencing, we found that multiple driver mutations in POLE‐mutant cancers show the characteristic POLE mutational signature, including those in genes conventionally regarded as initiators of tumourigenesis. In POLE‐mutant cancers, the proportion of monoclonal predicted neoantigens was similar to that in other cancers, but the absolute number was much greater. We also found that the prominent CD8+ T‐cell infiltrate present in POLE‐mutant cancers was evident in their precursor lesions. Collectively, these data indicate that somatic POLE mutations are early, quite possibly initiating, events in the endometrial and colorectal cancers in which they occur. The resulting early onset of genomic instability may account for the striking immune response and excellent prognosis of these tumours, as well as their early presentation.

Evolutionary history of human colitis-associated colorectal cancer

Objective IBD confers an increased lifetime risk of developing colorectal cancer (CRC), and colitis-associated CRC (CA-CRC) is molecularly distinct from sporadic CRC (S-CRC). Here we have dissected the evolutionary history of CA-CRC using multiregion sequencing. Design Exome sequencing was performed on fresh-frozen multiple regions of carcinoma, adjacent non-cancerous mucosa and blood from 12 patients with CA-CRC (n=55 exomes), and key variants were validated with orthogonal methods. Genome-wide copy number profiling was performed using single nucleotide polymorphism arrays and low-pass whole genome sequencing on archival non-dysplastic mucosa (n=9), low-grade dysplasia (LGD; n=30), high-grade dysplasia (HGD; n=13), mixed LGD/HGD (n=7) and CA-CRC (n=19). Phylogenetic trees were reconstructed, and evolutionary analysis used to reveal the temporal sequence of events leading to CA-CRC. Results 10/12 tumours were microsatellite stable with a median mutation burden of 3.0 single nucleotide alterations (SNA) per Mb, ~20% higher than S-CRC (2.5 SNAs/Mb), and consistent with elevated ageing-associated mutational processes. Non-dysplastic mucosa had considerable mutation burden (median 47 SNAs), including mutations shared with the neighbouring CA-CRC, indicating a precancer mutational field. CA-CRCs were often near triploid (40%) or near tetraploid (20%) and phylogenetic analysis revealed that copy number alterations (CNAs) began to accrue in non-dysplastic bowel, but the LGD/HGD transition often involved a punctuated ‘catastrophic’ CNA increase. Conclusions Evolutionary genomic analysis revealed precancer clones bearing extensive SNAs and CNAs, with progression to cancer involving a dramatic accrual of CNAs at HGD. Detection of the cancerised field is an encouraging prospect for surveillance, but punctuated evolution may limit the window for early detection.

Reply: Uncertainties in tumor allele frequencies limit power to infer evolutionary pressures.

Reply: Uncertainties in tumor allele frequencies limit power to infer evolutionary pressures

Reply to ‘Revisiting signatures of neutral tumor evolution in the light of complexity of cancer genomic data’

Williams et al. reply — Balaparya and De1 question the applicability of the power-law neutral-evolution model to adequately describe the pattern of subclonal somatic mutations in bulk cancer sequencing data. The authors’ letter focuses on the issues of the inherent noise in next-generation sequencing data, whereby random sampling of alleles, PCR amplification during library preparation, limited depth sequencing, and subclonal copy number changes may cause considerable uncertainty in variantallele frequency (VAF) measurement. The authors suggest that these errors lead to VAF measurements that, owing to overdispersion, follow a beta-binomial and not a binomial distribution. We thank Balaparya and De for the insightful comments and address their points in the response below. The issue of VAF measurement accuracy is a very important point and something that concerned us in our original study2. For this reason, we provided extensive error-propagation analysis in our original manuscript to identify the inherent biases that affect VAF estimation (Methods and equations (12)–(14) in ref. 2). We aimed at starting from the analytical form of neutral evolution (equation (7) in ref. 2) as the expected signal (S) and adding the different sources of noise (N), such as purity and allele sampling during library preparation, to generate the expected pattern S + N reported by the data. Our results demonstrate that the signature of neutral evolution is detectable with moderately high sequencing depth (≥ 100× ; Methods and Supplementary Fig. 10 in ref. 2), and we fully acknowledged that the signature of neutral evolution versus selection cannot be reasonably extracted (or rejected) from lower-depth datasets. In an effort to address Balaparya and De’s concerns, we tested the ability of our model to recover neutral evolutionary dynamics in the presence of beta-binomially distributed noise, and we found no significant differences with respect to the binomial noise used in our original manuscript, although with very high dispersion (ρ = 0.1), a degree of difference was appreciable (Fig. 1a). Moreover, we estimated the degree of overdispersion in the data that we analyzed in ref. 2 by fitting a beta-binomial model to the clonal cluster by using Markovchain Monte Carlo inference. In both the 100× whole-genome gastric cancer3 and whole-exome colon cancer4 data, we estimated the dispersion parameter ρ to be < 0.005 (Fig. 1b,c, respectively), a value notably 10× lower than postulated by Balaparya and De (Fig. 1c,d in ref. 1). Given that as ρ→ 0, the beta-binomial distribution converges to a binomial distribution, we argue that using a binomial distribution to model noise in sequencing data was appropriate in our original analysis. Balaparya and De also suggest that, because copy number alterations affect VAF distributions, very strict thresholds are necessary to ensure that regions analyzed with our method are truly diploid. This is an important point, and we concur that the original threshold of absolute log R ratio ≤ 0.5 may have been too lenient. To test the effect of this confounding factor, we reanalyzed the TCGA pan-cancer dataset by using the new publicly accessible Genomic Data Commons Data Portal (see URLs) variant calls, which were not available at the time of our original manuscript.

Reply to ‘Currently available bulk sequencing data do not necessarily support a model of neutral tumor evolution’

Werner et al. reply — In their correspondence, McDonald et al.1 question our assertion that the distribution of mutations in tumor bulk sequencing data suggests an underlying neutral evolutionary process in a proportion of cancers2 and instead propose alternative explanations that incorporate subclonal selection. We agree with the authors’ demonstration that it is possible, in principle, to construct models of selection that produce patterns similar to the neutral model. However, the key issue is whether the proposed models of selection are realistic, meaningful and, most importantly, more appropriate than the null neutral model. Before examining this issue, we first note that we extensively stressed in the original manuscript2 that the majority of cases we examined were not consistent with neutral evolution (~70% appeared non-neutral), and we did specifically cite Gerlinger et al.3 as an example of data dominated by selection2. Our finding that the majority of cancers do show evidence of subclonal selection is consistent with previous literature, including the cases highlighted by McDonald et al.3,4. Arguably, clonal evolution results from the interplay of three fundamental processes: random alterations (genetic, epigenetic, etc.), random drift and nonrandom selection, the third of which is the most complex to define and model. In the established field of population genetics, extensive effort has been dedicated to modeling the first two processes without selection, the so-called neutral dynamics5–7. This includes the development of entire statistical frameworks based on neutrality, such as coalescent theory8. On the contrary, models that include selection, especially in growing populations, have been much harder to derive analytically owing to the large number of assumptions in the definition of selection, including whether selection is clone intrinsic or clone extrinsic (microenvironmentally defined) and whether the magnitude of selection is constant or fluctuates in response to population dynamics. Importantly, most models of selection describe cancer dynamics in terms of time9,10 (for example, time to fixation of a selected mutant) and therefore, although insightful, are hard to apply to cancer genomic data where temporal dynamics are often unobservable. In light of this complexity, in our study, we asked the simple question of what happens to the mutations in a growing tumor in the case where only the first two processes above, namely random mutations and drift, are operating. This leads to a relatively simple model that is analytically tractable, wherein subclonal mutations accumulate following a 1/f cumulative distribution2. We note that this is the underlying solution of the fully stochastic Luria–Delbrück model, as previously demonstrated11,12. Importantly, this model is based on the ‘null hypothesis’ of molecular evolution in cancer13–15 and predicts what the absence of subclonal selection should look like in a growing tumor. We tested this hypothesis against subclonal mutations from a large body of sequencing data and found that in about 30% of cases we could not reject this null hypothesis, at least within the resolution of the currently available data. In their correspondence, McDonald et al.1 propose a more complex scenario that includes ongoing selection and report that in some cases their model also fits the 1/f cumulative distribution. First, we examine the fit of their proposed model to the data and highlight that considering the stochastic nature of selected mutants would change the interpretation of their analysis. Second, we discuss the distinction between evaluating the power of a test and the limitations of the information content in the data to which the test is applied, in this case singlesample bulk sequencing. Third, we analyze the plausibility of the authors’ biological assumptions underlying their model. In the correspondence by McDonald et al.1, neutrality was correctly rejected in a considerable proportion of simulations with subclonal selection (R2 < 0.98; their Fig. 1b). The exact proportion of cases incorrectly classified as neutral is not reported, but a few specific examples are shown in their Fig. 1c–f. Importantly, in those cases, the mutant proportion at the time of sampling is not reported nor is the time when the mutant was introduced. Both are key factors in judging the strength of the selection signal for two reasons: (i) in the case of strong and early selection, wherein a selected mutant sweeps to fixation, the evolutionary dynamics revert to neutral, and hence accepting the null for the final tumor is correct (as all cells in the tumor bear the selected mutation, so there is no subclonal selection) and (ii) because of the inherent stochasticity of the evolutionary process, selected mutants can either occur too late to grow to a detectable size or be weakly selected such that the clonal population of the tumor remains virtually unchanged with respect to the neutral expectation. Judging from Fig. 2a, this seems to be what happens often: most mutants have fitness slightly higher than 1 (where 1 is neutral) and many have fitness even lower than 1 (should be negatively selected), but all persist in the population. In such a model, it is clear that selection is not sculpting the population by removing unfit clones and benefitting fitter ones, as any mutant—fit or unfit—seems to survive. Thus, the dynamics described in the models of McDonald et al.1 are ‘effectively neutral’, and relatedly, it is not surprising that deviations from neutrality are undetectable. We highlight that it is fundamentally important to consider the size of differentially selected subclones when considering whether or not a tumor can be classified as neutrally evolving. In the authors’ second simulation model (their Fig. 2), many clones arise very late and are therefore undetectable in the data (high frequency of red dots representing a clone size of one cell in their Fig. 2a). We argue that no test will ever be able to detect a subclone made of a single cell in a whole malignancy—and indeed, it is debatable whether a clone of size 1 can even be considered to have been selected. We discuss the detection limits imposed by current data in our original manuscript (Fig. 5)2, as well as in subsequent work16,17. To demonstrate the impact of subclone size in determining whether a tumor is classified as (effectively) neutral or not, we performed a more thorough analysis of our previous model of a stochastic branching process under selection (Fig. 1 in this letter). These simulations show that, in the presence of a subclone of detectable size in the data (for example, one that is not too small to be out of the detectable range of the variant

Reply to ‘Neutral tumor evolution?’

Impact of clonal copy number alterations In our previous study2, we assessed the cumulative variant allele frequency (VAF) distribution M(f) over the frequency range f = [0.12,0.24] to restrict our analysis to subclonal variants within a range applicable to the diverse datasets that we considered. Tarabichi and colleagues note that tumors with a tetraploid genome will have an additional ‘peak’ of clonal mutations at f ~0.25 (mutations in a single allele, Supplementary Fig. 1a), thus causing incorrect rejection of neutrality (Supplementary Fig. 1b). The integration range that we chose was based on a triploid tumor with read depth of 100× , thereby resulting in an upper threshold of 0.26 (Supplementary Methods). Although this procedure is suitable in most cases, it is not suitable for a tetraploid tumor, thus suggesting that the number of tumors consistent with neutral evolution could be larger than we reported. We show how this problem can be avoided by adjusting the range for tetraploid tumors (Supplementary Fig. 1c). We do acknowledge that the 1/f integration method is more accurate when applied to the entire VAF spectrum of subclonal mutations only. Moreover, we have recently developed a Bayesian model selection framework that compares the neutral model against models with selection, using the entire VAF distribution3. We do stress, however, that most cancers analyzed in our original manuscript were not neutral and showed signs of subclonal selection.

Functional versus non-functional intratumor heterogeneity in cancer.

ABSTRACT Next-generation sequencing data from human cancers are often difficult to interpret within the context of tumor evolution. We developed a mathematical model describing the accumulation of mutations under neutral evolutionary dynamics and showed that 323/904 cancers (∼30%) from multiple types were consistent with the neutral model of tumor evolution.