Quantifying chromosomal instability from intratumoral karyotype diversity using agent-based modeling and Bayesian inference

  1. Andrew R Lynch
  2. Nicholas L Arp
  3. Amber S Zhou
  4. Beth A Weaver
  5. Mark E Burkard  Is a corresponding author
  1. Carbone Cancer Center, University of Wisconsin-Madison, United States
  2. McArdle Laboratory for Cancer Research, University of Wisconsin-Madison, United States
  3. Department of Cell and Regenerative Biology, University of Wisconsin, United States
  4. Division of Hematology Medical Oncology and Palliative Care, Department of Medicine University of Wisconsin, United States

Abstract

Chromosomal instability (CIN)—persistent chromosome gain or loss through abnormal mitotic segregation—is a hallmark of cancer that drives aneuploidy. Intrinsic chromosome mis-segregation rate, a measure of CIN, can inform prognosis and is a promising biomarker for response to anti-microtubule agents. However, existing methodologies to measure this rate are labor intensive, indirect, and confounded by selection against aneuploid cells, which reduces observable diversity. We developed a framework to measure CIN, accounting for karyotype selection, using simulations with various levels of CIN and models of selection. To identify the model parameters that best fit karyotype data from single-cell sequencing, we used approximate Bayesian computation to infer mis-segregation rates and karyotype selection. Experimental validation confirmed the extensive chromosome mis-segregation rates caused by the chemotherapy paclitaxel (18.5 ± 0.5/division). Extending this approach to clinical samples revealed that inferred rates fell within direct observations of cancer cell lines. This work provides the necessary framework to quantify CIN in human tumors and develop it as a predictive biomarker.

Editor's evaluation

The authors have developed a framework to quantify rates of chromosomal instability (CIN) in human tumors by fitting karyotype distributions inferred from low-depth DNA-sequencing to in silico models of CIN with karyotype selection pressures, sweeping through parameter space. This is particularly useful for the development of biomarkers for CIN, which is associated with cancer metastasis and drug resistance.

https://doi.org/10.7554/eLife.69799.sa0

eLife digest

DNA contains all the information that cells need to function. The DNA inside cells is housed in structures called chromosomes, and most healthy human cells contain 23 pairs. When a cell divides, all chromosomes are copied so that each new cell gets a complete set. However, sometimes the process of separating chromosomes is faulty, and new cells may get incorrect numbers of chromosomes during cell division. Cancer cells frequently exhibit this behavior, which is called chromosomal instability’, or CIN.

Chromosomal instability affects many cancer cells with varying severity. In cancers with high chromosomal instability, the number of chromosomes may change almost every time the cells divide. These cancers are often the most aggressive and difficult to treat.

Scientists can estimate chromosomal instability by counting differences in the number of chromosomes across many cells. However, many cells that are missing chromosomes die, resulting in inaccurate measures of chromosomal instability. To find a solution to this problem, Lynch et al. counted chromosomes in human cells with different levels of chromosomal instability and created a computer model to work out the relationship between chromosomal instability and chromosome number.

The model could account for both living and dead cells, which gave more accurate results. Lynch et al. then confirmed the accuracy of their approach by using it on a group of cells treated with a chemotherapy drug that causes a known level of chromosomal instability. They also used existing data from breast and bowel cancer, which revealed that levels of chromosomal instability varied between one mistake per three to twenty cell divisions.

Lower levels of chromosomal instability can be linked to a better prognosis for cancer patients, but it currently cannot be measured reliably. These results may help to reveal the causes of chromosomal instability and the role it has in cancer. If this method is successfully applied to patient samples, it could also improve our ability to predict how each cancer will progress and may lead to better treatments.

Introduction

Chromosomal instability (CIN) is characterized by persistent whole-chromosome gain and loss through mis-segregation during cell division. Genome instability is a hallmark of cancer (Hanahan and Weinberg, 2011) and one type, CIN, is the principal driver of aneuploidy, a feature of ~80% of solid tumors (Hancock et al., 2004; Knouse et al., 2017; Weaver and Cleveland, 2006). CIN potentiates tumorigenesis (Foijer et al., 2017; Levine et al., 2017; Silk et al., 2013) and associates with therapeutic resistance (Ippolito et al., 2020; Lee et al., 2011; Lukow et al., 2020; Pavelka et al., 2010), metastasis (Bakhoum et al., 2018) and poor survival outcomes (Bakhoum et al., 2011; Denu et al., 2016; Jamal-Hanjani et al., 2017). Thus, CIN is an important characteristic of cancer biology. Despite its importance, CIN has not emerged as a clinical biomarker, in part because it is challenging to quantify.

Although CIN has classically been characterized as binary—tumors either have it or not—recent evidence highlights the importance of the rate of chromosome mis-segregation and the specific aneuploidies it produces. For example, clinical outcomes partially depend on aneuploidy of specific chromosomes (Davoli et al., 2013; Sheltzer et al., 2017; Vasudevan et al., 2020). Further, higher levels of CIN suppress tumor growth when they surpass a critical threshold, thought to be due to lethal loss of essential genes and irregular expression due to imbalanced gene dosage (Funk et al., 2021; Silk et al., 2013; Weaver and Cleveland, 2008; Zasadil et al., 2014). Moreover, baseline CIN may predict chemotherapeutic response to paclitaxel (Janssen et al., 2009; Swanton et al., 2009) and is proposed to both promote detection by or evasion from the immune system (Davoli et al., 2017; Santaguida et al., 2017). No single or standardized analytically valid measure of CIN has emerged and this gap has precluded its clinical validation as a prognostic or predictive biomarker.

Prior measures of CIN use various means to compare levels in tumors or populations, but do not establish a standardized quantitative rate. These prior measures include histologic analysis of mitotic defects (Bakhoum et al., 2011; Jin et al., 2020), fluorescence in situ hybridization (FISH) with probes to detect individual chromosomes (Thompson and Compton, 2008), and gene-expression methodologies like CIN scores (Carter et al., 2006). While these methods are readily accessible, they have significant drawbacks for clinical application. FISH and mitotic visualization approaches are laborious. Direct visualization of mitotic defects to measure CIN is only possible in the most proliferative tumors where enough cells are captured in short-lived mitosis. FISH typically quantifies only a subset of chromosomes, which will be misleading if there is bias toward specific chromosome gains/losses (Dumont et al., 2020). While gene expression scores are proposed as indirect measures of CIN, they are not specific to CIN and correlate highly with proliferation and structural aneuploidy (Carter et al., 2006; Sheltzer, 2013).

Single-cell sequencing promises major advances in quantitative measures of CIN by displaying cell-cell variation for each chromosome across hundreds of cells (Navin et al., 2011; Wang et al., 2014). However, selection poses another complication. To date, single-cell analyses have identified surprisingly low cell-cell karyotype variation, even when mitotic errors are directly observed by microscopy (Bolhaqueiro et al., 2019; Gao et al., 2016; Kim et al., 2018; Nelson et al., 2020; Wang et al., 2014). These observations highlight the confounding role of selection against aneuploid karyotypes in measuring CIN in human tumors. Indeed, selection reduces karyotype variance in cancer cell populations that directly exhibit mitotic errors (Gerstung et al., 2020; Ippolito et al., 2020; Lukow et al., 2020). Here, we seek to overcome gap by modeling chromosomal instability and explicitly considering the evolutionary selection of aneuploid cells, to derive a quantitative measure.

We describe a quantitative framework to measure CIN by sampling population structure and cell-cell karyotypic variance in human tumors, accounting for selection on aneuploid karyotypes. We built our framework on the use of phylogenetic topology measures to quantify underlying evolutionary processes (Mooers and Heard, 1997); in this case to quantify CIN from both the diversity and the aneuploid phylogeny within a tumor. Using an agent-based model of CIN, we determine how distinct types and degrees of selective pressure shape the karyotype distribution and population structure of tumor cells at different rates of chromosome mis-segregation. We then use this in silico model as a foundation for parameter inference to provide a quantitative estimate of CIN as the numerical rate of chromosome mis-segregation per cell division. We apply this model to quantify CIN caused by the chemotherapeutic paclitaxel in culture. Next, using existing single-cell whole-genome sequencing data (scDNAseq), we measure CIN in cancer biopsy and organoid samples. As a whole, this work provides a framework to quantify CIN in human tumors, a first step toward developing CIN as a prognostic and predictive biomarker.

Results

A framework for modeling CIN and karyotype selection

To assess intratumoral CIN via cell-cell karyotype heterogeneity, we considered how selection on aneuploid karyotypes impacts observed chromosomal heterogeneity within a tumor. By modeling fitness of aneuploid cells, we observe chromosomal variation in a population of surviving cells. The selective pressure of diverse and specific aneuploidies on human cells has not been, to our knowledge, directly measured. Therefore, we employ previously developed models of selection.

In models of CIN, fit karyotypes are selected while unfit aneuploid karyotypes are eliminated over time (Ippolito et al., 2020; Ravichandran et al., 2018; Sheltzer et al., 2017; Vasudevan et al., 2020). We use two previously proposed models of aneuploidy-associated cellular fitness, as well as hybrid and neutral selection models. The Gene Abundance model is based on the relatively low incidence of aneuploidy in normal tissues and assumes cellular fitness declines as the cell’s karyotype diverges from a balanced euploid karyotype (Sheltzer and Amon, 2011; Zhu et al., 2012). When an individual chromosome diverges from euploid balance (2 N, 3 N, 4 N, for example), its contribution to cellular fitness is weighted by its abundance of genes (Figure 1—figure supplement 1A, left). Alternatively, the Driver Density model assumes that each chromosome’s contribution to cellular fitness is weighted by its ratio of Tumor suppressor genes, Oncogenes, and Essential genes (TOEs)(Davoli et al., 2013; Laughney et al., 2015). For example, Driver Density selection will favor loss of chromosomes with many tumor suppressors and favor gain of chromosomes replete with oncogenes and essential genes (Figure 1—figure supplement 1A, right). The hybrid averaged model accounts for both karyotypic balance and TOE densities (Figure 1—figure supplement 1A, middle). Using these fitness models, we assigned chromosome scores to reflect each chromosome’s value to cellular fitness (Figure 1—figure supplement 1B, Table 1), the sum of which represent the total fitness value for the cell, relative to a value of 1 for a euploid cell. Further, we scaled the impact of cell fitness with a scaling factor, S, ranging from 0 (no selection) to 100 (high selection). While these models are approximations, they are nevertheless useful to estimate how mis-segregation and selective pressure cooperate to mold karyotypes in the cell population.

Table 1
Base chromosome-specific fitness scores for individual models.
Selection model
CHR ARMGene AbundanceDriver DensityHybrid
1p0.04780162–0.00240180.02269992
1q0.043403210.032443620.03792341
2p0.027336550.029357170.02834686
2q0.042440540.039432670.0409366
3p0.023104120.032896950.02800053
3q0.02997560.054167360.04207148
4p0.012381950.017849090.01511552
4q0.031817960.029013240.0304156
5p0.011784430.042811660.02729805
5q0.037876150.019499340.02868775
6p0.025577190.023986190.02478169
6q0.025543990.000116250.01283012
7p0.01795880.098892840.05842582
7q0.032315890.069333140.05082451
8p0.015917280.027695640.02180646
8q0.02549420.058614270.04205423
9p0.01301266–0.00129410.00585929
9q0.025726570.047026810.03637669
10 p0.0112201–0.0364218–0.0126008
10q0.027502530.011426880.01946471
11 p0.019618580.038186210.0289024
11q0.036299360.018987840.0276436
12 p0.01425750.05515510.0347063
12q0.036598120.062737860.04966799
13 p000
13q0.02333649–0.01015390.00659128
14 p1.66E-0508.30E-06
14q0.037925940.025574390.03175016
15 p000
15q0.037013060.02065660.02883483
16 p0.023834420.043347360.03359089
16q0.01900446–0.00714440.00593005
17 p0.01548573–0.00859750.00344414
17q0.035535860.043634740.0395853
18 p0.006273960.005336970.00580547
18q0.01434049–0.0263632–0.0060113
19 p0.021593720.053714160.03765394
19q0.028133250.005503380.01681831
20 p0.00896280.043510250.02623653
20q0.015269960.049935930.03260295
21 p0.0023236900.00116185
21q0.01233215–0.00330920.00451147
22 p0.0001327806.64E-05
22q0.02297134–0.00515810.0089066
Xp0.0155521300.00777606
Xp0.0249962700.01249813

We employed these selection models in an agent-based model of exponential population growth wherein each cell has its own karyotype (Figure 1 and Figure 1—figure supplement 1). Briefly, simulations started with 100 euploid cells and were run in discrete time steps with variable rates of selective pressure, S, and rates of chromosome mis-segregation (Pmisseg, see definitions in Table 2). The rate—or probability—of mis-segregation events, Pmisseg, is the measure of CIN. During each time step, cells have a Pdivision ( = 0.5 for euploid) chance of dividing. Each dividing cell has a Pmisseg chance of improper segregation of each chromosome. Segmental chromosome breaks occur with a probability Pbreak, set at 0 or 0.5. After division, fitness (F) of each daughter is assessed. Cells are removed from the population if any given chromosome has copy number 0 or >6. The Pdivision value of the remaining viable cells is adjusted by the cell’s fitness under selection (FS). Due to computational limitations, pseudo-Moran or Wright-Fisher models are employed to limit the modeled cell population (Figure 1—figure supplement 1C, D). These limits did not significantly affect the measures extracted from these populations (Figure 1—figure supplement 2). Thus, these models simulate an evolving population of aneuploid cells under given rates of CIN, Pmisseg, and models and strength of selection.

Table 2
Parameters varied during agent-based modeling.
ParameterDescription
PmissegProbability of mis-segregation per chromosome per division
PbreakProbability of chromosome breakage after mis-segregation
PdivisionProbability of cellular division per time step
SMagnitude of selective pressure on aneuploid karyotypes
Figure 1 with 2 supplements see all
A framework for modeling CIN and karyotype selection.

(A) Chromosome arm scores for each model of karyotype selection. Gene Abundance scores are derived from the number of genes per chromosome arm normalized to the number of all genes. Chromosome arms 13 p and 15 p did not have an abundance score and were set to 0. Driver Density scores come from the pan-cancer chromosome arm scores derived in Davoli et al., 2013, and normalized to the sum of chromosome arm scores for chromosomes 1-22,X. Chromosome arms 13 p, 14 p, 15 p, 21 p, 22 p, and chromosome X did not have driver scores and were set to 0. Hybrid model scores are set to the average of the Driver and Abundance models. The neutral model (not displayed) is performed with all cell’s fitness constitutively equal to 1 regardless of karyotype. (B) Framework for the simulation of and selection on cellular populations with CIN. Cells divide (Pdivision starts at 0.5 in the exponential pseudo-Moran model and is constitutively equal to 1 for the constant Wright-Fisher model) and probabilistically mis-segregate chromosomes (Pmisseg ∈ [0, 0.001… 0.05]). After, cells experience selection under one of the selection models, altering cellular fitness and the probability (Pdivision) a cell will divide again (green check). Additionally, cells wherein the copy number of any chromosome falls to zero or surpasses 6 are removed (red x). After this, the cycle repeats. See Materials and methods for further details.

Evolutionary dynamics is imparted by CIN

To understand the interplay between CIN and selection, we simulated 100 steps of cell growth with CIN under each selection model. We varied the rate of CIN (Pmisseg,c ∈[0, 0.001… 0.05] per chromosome; or 0–2.3 chromosome mis-segregations per division) and selective pressure ranging from none to heavy selection (S ∈[0, 2… 100]). As expected, the simulated cell number increases rapidly to the pseudo-Moran cap of 3000, where it remains (Figure 2A). As displayed in Figure 2B, diversity of the cell population, expressed as mean karyotypic variance increases over time, but also depends on mis-segregation rate, and selection levels (Figure 2B). As expected, high mis-segregation rates (Pmisseg, Y axis) and low selection (S = 0; top row) enhance the variance of the population. Further, without selection (S = 0; top row) all models returned comparable profiles over time, resembling neutral selection. However, when selective pressure is applied (S > 0), the distinct profiles appear. The abundance model (first column) negatively selects against all aneuploid karyotypes and yields low heterogeneity that increases modestly with mis-segregation rate. With the Driver model (second column), there is a sharp increase in heterogeneity even at low mis-segregation rates, as this model favors specific aneuploid states that maximizes oncogenes and minimizes tumor suppressors. The Hybrid model falls between the other two. Results were not specific to the pseudo-Moran process of capping at 3000 cells—dynamics were similar in the constant-population Wright-Fisher model (Figure 2—figure supplement 1A, B). These data illustrate how CIN and selection operate together to shape the karyotype diversity in the cell population.

Figure 2 with 2 supplements see all
Evolutionary dynamics imparted by CIN.

(A) Population growth curve in the absence of selective pressure (Pmisseg = 0.001, S = 0, n = 3 simulations). The steady state population in null selection conditions is 3000 cells. (B) Heatmaps depicting dynamics of karyotype diversity as a function of time (steps), mis-segregation rate (Pmisseg), and selection (S) under each model of selection. Columns represent the same model; rows represent the same selection level. Mean karyotype diversity (MKV) is measured as the variance of each chromosome averaged across all chromosomes 1–22, and chromosome X. Low and high MKV are shown in white and blue respectively (n = 3 simulations for every combination of parameters). (C) Population growth under each model, varying Pmisseg and S. Pmisseg∈ [0.001, 0.022, 0.050] translate to about 0.046, 1, and 2.3 mis-segregations per division respectively for diploid cells. (D) Dynamics of the average ploidy (total # chromosome arms / 46) of a population while varying Pmisseg and S. (E) Dynamics of ploidy under each model for diploid and tetraploid founding populations. Pmisseg∈ [0.01, 0.02] translate to about 0.46 and 0.92 mis-segregations for diploid cells and 0.92 and 1.84 mis-segregations for tetraploid cells. (F) Fitness (FS) over time for diploid and tetraploid founding populations evolved under each model. (G) Karyotype diversity dynamics for diploid and tetraploid founding populations. MKV is normalized to the mean ploidy of the population at each time step. Plotted lines in C-G are local regressions of n = 3 simulations.

High levels of selection against aneuploid cells are expected to impede cell growth. To visualize this, we quantified the population of viable cells with distinct models (Figure 2C). As expected with the Abundance model at S = 10 and S = 100, cells proliferated more slowly with higher rates of mis-segregation. By contrast, the Driver model saw no growth defect as they favored specific aneuploid states that are easily reached with missegregation. As before, the Hybrid model, is intermediate, and findings are not impacted by pseudo-Moran or Wright-Fisher restrictions on cell number (Figure 2—figure supplement 1C).

To further assess model dynamics, we examined time-course of average cellular ploidy—the number of chromosomes divided by 23. In many cases, the mean ploidy of the populations tend to increase over time (Figure 2D, Figure 2—figure supplement 1D), particularly in the absence of selection (S = 0; top). This is likely due to a higher permissiveness to chromosome gains than losses in our model (since cells ‘die’ with nullisomy or any chromosome >6, the optimum is 3.0). With selection (S = 10; S = 100 rows), the models diverge. In the abundance model, populations remain near diploid. With the Driver model, the average ploidy increases more rapidly due to favoring aneuploidy states that favor high oncogenes and low tumor suppressors, consistent with previous computational models built on chromosome-specific driver densities (Davoli et al., 2013; Laughney et al., 2015). Under the Hybrid model, ploidy increases modestly. Similar effects are seen with the constant-population Wright-Fisher model (Figure 2—figure supplement 1D). In sum, selection and mis-segregation cooperate to shape the aneuploid karyotypes diversity, cell proliferation and average ploidy in a population of cells, or a human tumor. Further, sampling karyotypes in a cell population does not allow direct determination of mis-segregation rates, as their diversity is influenced by other factors such as selective pressure, selection modality, and time.

In some tumors, genome doubling occurs early in tumor initiation relative to other copy number changes (Bielski et al., 2018; Gerstung et al., 2020). Genome doubling is accomplished, for example, by endoreduplication, by failed cytokinesis, or by cell-cell fusion. Genome doubling buffers against loss of chromosomes and thereby favors aneuploidy. To determine how genome doubling impacts evolution in our model, we compared diploid and tetraploid founders (Figure 2E–G). Both diploids and tetraploids tend to converge toward the near-triploid state (ploidy ~3), as observed in many human cancers (Carter et al., 2012), although this is restrained to a degree with the Abundance and Hybrid models. Compared with diploid cells, tetraploidy buffered against the negative effects of cellular fitness in the Abundance model, despite generating similar levels of diversity over time (Figure 2F and G)— this is more pronounced when comparing Pmisseg = 0.1 in tetraploids versus Pmisseg = 0.2 in diploids to match the number of chromosome mis-segregations per division. This is consistent with the idea that tetraploidy serves as an intermediate enabling a near-triploid karyotype that is common in many cancers (Bielski et al., 2018; López et al., 2020). By contrast, in the Driver model, tetraploidy did not provide a selective advantage to high-CIN tumors (Figure 2F). Similar fitness, karyotype diversities, and ploidy increases were obtained with a Wright-Fisher model of population growth (Figure 2—figure supplement 1E-G, Figure 2—figure supplement 2).

Taken together, the agent-based model recapitulates expected key aspects of tumor evolution, lending credence to our model. Further, they illustrate the difficulty of inferring mis-segregation rates directly from assessing variation in karyotypes in human cancer. Nevertheless, this model provides a framework to incorporate selection to measure CIN through quantitative inference from the observed karyotypes, as we will demonstrate.

Long-term karyotype diversity depends profoundly on selection modality

Some current measures of CIN are derived from karyotype diversity in the population. Yet, our model suggests that selection pressure will profoundly shape this diversity. To further understand the nature of karyotype diversity under selection, we evaluated their long-term dynamics, whether they exhibit clonality, and whether populations simulated under each model converge on a common karyotype.

We simulated diploid and tetraploid populations for 3000 time steps at a fixed mis-segregation rate, in an experimentally reported range, allowing for fragmentation of chromosome arms (Pmisseg = 0.003, Pbreak = 0.5) (Bakhoum et al., 2009; Bolhaqueiro et al., 2019; Weaver et al., 2007) and S ∈ [1,25] (Figure 3A). We visualized copy-number heatmaps indicating karyotypes of sampled cells from the population. As expected, population diversity is limited under the Abundance model (Figure 3B). Even after 3000 time steps, only a small number of unique alterations and sub-clonal alterations ( + 13 p/–15 p/–22 p) existed, likely passenger alterations as they offer no fitness advantage in this model. Moreover, the karyotype average of 1500 cells across five replicates resembled a diploid karyotype (Figure 3C, row 1), indicating that the Abundance model provides stabilizing selection around the euploid karyotype. In fact, populations simulated under this model with elevated selection (S = 25) quickly reach a low, steady-state level of karyotype diversity and fitness while those with the unmodified selection values (S = 1) take a longer time to reach this steady-state and have similar levels of karyotype diversity and fitness as the other models (Figure 3—figure supplement 1). To identify any contingencies that may affect these associations, we performed the same simulation using several variants of our model. We found this steady state to be consistent for tetraploid cells as well as when we eased the upper ploidy constraint from nc c = 6 to an extreme nc c = 10, when we imposed a severe, 90% fitness reduction for all cells with a haploidy, and when we simulated populations under the Wright-Fisher model (Figure 3C, rows 2–4).

Figure 3 with 1 supplement see all
Karyotype diversity depends profoundly on selection modality.

(A) Simulation scheme to assess long-term dynamics of karyotype evolution and karyotype convergence. (B) Heatmaps depicting the chromosome copy number profiles of a subset (n = 30 out of 300 sampled cells) of the simulated population with early CIN over time under each model of karyotype selection. (C) Average heatmaps (lower) show the average copy number across the 5 replicates for (1) the Exponential Psuedo-Moran (Base), (2) the base model with the upper copy number limit set to 10, (3) the base model that invokes a FM x 0.1 penalty for any cell with a haploid chromosome, (4) and the Constant Population-Size Wright-Fisher model. Pmisseg = 0.003; S = 25 (except Neutral model; S = 0); ploidy = 2.

The Driver Density and Hybrid models generate much more diversity (Figure 3B) but nevertheless converge by 3,000 timesteps (Figure 3—figure supplement 1). Without selection (neutral model), there is high diversity and no convergence over time. Taken together, these demonstrate a high dependence on the model of selection. However, the models are not highly dependent on ploidy constraints, haploid penalties, or on selection of Pseudo-Moran or Wright-Fisher restriction of cell numbers. Taken together, long-term populations are strongly shaped by the model of karyotype selection for a given Pmisseg, but relatively insensitive to other particular features of the model. This justifies our approach henceforth of varying only the selection model, the degree of selection (S), and Pmisseg to infer parameters from data via phylogenetic topology and Bayesian inference.

Topological features of simulated phylogenies delineate CIN rate and karyotype selection

Given a model capable of recapitulating diversity and selective pressures, next we wish to infer Pmisseg as a measure of CIN from an observed population of cells. Phylogenetic trees provide insights into evolutionary processes of genetic diversification and selection. Moreover, the topology of the phylogenetic tree has been used as a quantitative measure of the underlying evolutionary processes (Colijn and Plazzotta, 2018; Dayarian and Shraiman, 2014; Manceau et al., 2015; Neher et al., 2014; Scott et al., 2020).

Here, chromosome mis-segregation gives rise to karyotype heterogeneity, and the population of cells is then shaped by selection. To evaluate this, we use chromosome copy number-based phylogenetic reconstruction, since mutation rates are not high enough in tumors to reliably infer cellular relationships, particularly with low-copy sequencing. Once phylogenies are reconstructed from simulated and experimental populations, the topological features phylogenies can be compared. These features include ‘cherries’—two tips that share a direct ancestor—and ‘pitchforks—a clade with three tips (Figure 4A). Additionally, we considered a broader metric of topology, the Colless index, which measures the imbalance or asymmetry of the entire tree. To understand how these measures are affected by selection in simulated populations, we reconstructed phylogenies from 300 random cells from each population simulated with a range of selective pressures taken at 60 time steps (~30 divisions under Hybrid selection; Figure 4B). As seen previously, aneuploidy and mean karyotypic variance (MKV) decrease with selective pressure, a trend that is robust at high mis-segregation rates (Figure 4C). By contrast, Colless indices increase with mis-segregation rates and selective pressures, as the resulting variation and selection generate phylogenetic asymmetry. Accordingly, this imbalance is apparent in phylogenetic reconstructions of simulated populations (Figure 4D). Cherries, by contrast, decrease with selection due to selection against many aneuploidies (Figure 4C). Pitchforks seemed less informative. Therefore, we tentatively selected 4 phylogenetic parameters that can retain information about chromosome missegregation—aneuploidy, MKV, Colless, and Cherries.

Topological features of simulated phylogenies delineate CIN rate and karyotype selection.

(A) Quantifiable features of karyotypically diverse populations. Heterogeneity between and within karyotypes is described by MKV and aneuploidy (inter- and intra-karyotype variance, see Materials and methods). We also quantify discrete topological features of phylogenetic trees, such as cherries (tip pairs) and pitchforks (3-tip groups), and a whole-tree measure of imbalance (or asymmetry), the Colless index. (B) Scheme to test how CIN and selection influence the phylogenetic topology of simulated populations. (C) Computed heterogeneity (aneuploidy and MKV) and topology (Colless index, cherries, pitchforks) summary statistics under varying Pmisseg and S values. MKV is normalized to the average ploidy of the population. Topological measures are normalized to population size. Spearman rank correlation coefficients (r) and p-values are displayed (n = 8 simulations). (D) Representative phylogenies for each hi/low CIN, hi/low selection parameter combination and their computed summary statistics. Each phylogeny represents n = 50 out of 300 cells for each simulation. (E) Dimensionality reduction of all simulations for each hi/low CIN, hi/low selection parameter combination using measures of karyotype heterogeneity only (left; MKV and aneuploidy) or measures of karyotype heterogeneity and phylogenetic topology (right; MKV, aneuploidy, Colless index, cherries, and pitchforks).

To characterize how well the four measures retain information about the simulation parameters, we performed dimensionality reduction with measures of karyotype heterogeneity alone (MKV and aneuploidy) alone and adding Colless and cherries—measures of phylogenetic topology (Figure 4E). This analysis indicates that when considering heterogeneity alone simulations performed under high CIN/high selection (yellow) and low CIN/low selection (red) associate closely, meaning these measures of heterogeneity are not sufficient to distinguish these disparate conditions (Figure 4E, left). These similarities arise because high selection can mask the heterogeneity expected from high CIN. By contrast, combining measures of heterogeneity with those of phylogenetic topology can discriminate between simulations with disparate levels of CIN and selection (Figure 4E, right). This provides further evidence that measures of heterogeneity alone are not sufficient to infer CIN due to the confounding effects of selection, particularly when the nature of selection is unclear or can vary. Together these results indicate that phylogenetic topology preserves information about underlying levels of selective pressure and rates of chromosome mis-segregation. Further, phylogenetic topology of single-cell populations may be a suitable way to correct for selective pressure when estimating the rate of chromosome mis-segregation from measures of karyotype diversity.

Experimental chromosome mis-segregation measured by Bayesian inference

To experimentally validate quantitative measures of CIN, we generated a high rate of chromosome mis-segregation with a clinically relevant concentration of paclitaxel (Taxol) over 48 hr (Figure 5A). We treated CAL51 breast cancer cells with either a DMSO control or 20 nM paclitaxel, which generates widespread aneuploidy due to chromosome mis-segregation on multipolar mitotic spindles (Zasadil et al., 2014), verified in this experiment (Figure 5—figure supplement 1A). At 48 hr cells will have undergone 1–2 mitoses and, consistent with abnormal chromosome segregation, we observe broadened DNA content distributions by flow cytometry (Figure 5—figure supplement 1B). Using low-coverage scDNAseq data, we characterized the karyotypes of 36 DMSO- and 134 paclitaxel-treated cells. As expected, virtually all cells had extensive aneuploidy after paclitaxel, in contrast with low variance in the control (Figure 5B). Additionally, the mean of the resultant aneuploid karyotypes for each chromosome still resembled those of bulk-sequenced cells, highlighting that bulk-sequencing is an ensemble average, and does not detect variation in population aneuploidy, particularly with balanced mis-segregation events (Figure 5B, single-cell mean and bulk). In quantifying the absolute deviation from the modal control karyotype in each cell, and assuming a single mitosis, cells exposed to 20 nM paclitaxel mis-segregate 18.5 ± 0.5—a Pmisseg of ~0.42 considering the control’s sub-diploid modal karyotype (Figure 5C). The majority of these appeared to be whole-chromosome mis-segregations (Figure 5—figure supplement 2).

Figure 5 with 5 supplements see all
Experimental chromosome mis-segregation measured by Bayesian inference experimental scheme.

(A) Cal51 cells were treated with either DMSO or 20 nM paclitaxel for 48 hr prior to further analysis by time lapse imaging, bulk DNA sequencing, and scDNAseq. (B) Heatmaps showing copy number profiles derived from scDNAseq data, single-cell copy number averages, and bulk DNA sequencing. (C) Observed mis-segregations calculated as the absolute sum of deviations from the observed modal karyotype of the control. (D) Dimensionality reduction analysis of population summary statistics (aneuploidy, MKV, Colless index, cherries) from the first three time steps of all simulations performed under the Hybrid model. (E) 2D density plot showing joint posterior distributions from ABC analysis using population summary statistics computed from the paclitaxel-treated cells using the following priors and parameters: Growth Model = ‘exponential pseudo-Moran’, Selection Model = ‘Hybrid, initial ploidy = 2, 2 time steps, S ∈[0, 2… 100], Pmisseg∈[0, 0.005… 1.00] and a tolerance threshold of 0.05 to reject dissimilar simulation results. (see Materials and Methods). Vertical dashed line represents the experimentally observed mis-segregation rate. White + represents the mean of inferred values.

In this instance, we were able to estimate mis-segregation rate by calculating absolute deviation from the modal karyotype after a single aberrant cell division. However, such an analysis would not be possible for long-term experiments, or real tumors, where new aneuploid cells may be subject to selection. Accordingly, we sought to infer the parameters of this experiment—the mis-segregation rate of 18.5 chromosomes per division and low selection—using only measures of aneuploidy, variance, and phylogenetic topology. To display this, we used dimensionality reduction to ensure that observed measures from the paclitaxel-treated Cal51 population fell within the space of those observed from simulated populations over 2 steps under the Hybrid model. The experimental data mapped to those from simulations using high mis-segregation rates and relatively low selection (red point, Figure 5D). However, this comparison does not provide a quantitative measure of CIN. Instead, parameter inference via approximate Bayesian computation (ABC) is well suited for this purpose.

By deriving phylogeny metrics from simulated populations under a wide-range of distributions of evolutionary parameters, ABC identifies evolutionary parameters most consistent with the data—the posterior probability distribution. We used ABC with simulated data to infer the chromosome mis-segregation rate and selective pressure in the paclitaxel-treated cells (Csilléry et al., 2012). Importantly, this data has directly observed rates of mis-segregation, which provide a gold standard benchmark to optimize ABC inference.

One key aspect of ABC is the selection of optimal phylogenetic summary statistics. A small number of summary statistics is optimal and larger numbers impair the model (Csilléry et al., 2012). To address this, a common approach is to identify a small set of summary statistics that achieve the best inference. Here, we used the experimentally observed mis-segregation rate as a benchmark to optimally select a panel of measures for parameter inference (Figure 5—figure supplement 3) and selected the following four metrics to use concurrently in our ABC analysis: mean aneuploidy, MKV, the Colless index (a phylogenetic balance index) and number of cherries (normalized to population size). In doing so, this analysis inferred a chromosome mis-segregation rate of 0.396 ± 0.003 (or 17.4 ± 0.1 chromosomes; mean ± SE), which compares favorably with the experimentally observed rate of 18.5 ± 0.5 (Figure 5E; dashed line represents experimental rate, white ‘+’ the inferred rate). The distribution of accepted values for selection was skewed toward lower pressure (21 ± 0.4; mean ± SE), meaning that karyotype selection had little bearing on the result at this time point, consistent with the absence of selection in a 48-hr experiment.

Interestingly, the incidence of nullisomy in the simulated population was higher than in the paclitaxel-treated populations at the observed mis-segregation rate (Figure 5—figure supplement 4A). This could be due to spindle pole clustering, a recovery mechanism often seen in paclitaxel-treated cells that causes non-random chromosome mis-segregations. A posterior predictive check of the summary statistics demonstrates how each contributes to the inference of CIN rate (Figure 5—figure supplement 4B). In short, this experimental case validated ABC-derived mis-segregation rate as a measure of CIN, with an experimentally determined mis-segregation rate. Importantly, prior estimations of mis-segregation rate selective pressure were not required to develop this quantitative measure of CIN.

Together, these data indicate that combining simulated and observed metrics of population diversity and structure with a Bayesian framework for parameter inference is a flexible method of quantifying the evolutionary forces associated with CIN. Moreover, this method reveals the hitherto unreported potential extent of chromosome mis-segregation induced by a clinically relevant concentration of the successful chemotherapeutic paclitaxel consistent with the measured mis-segregation from non-pharmacologically induced multipolar divisions (Bollen et al., 2021).

Minimum sampling of karyotype heterogeneity

The cost of high-throughput DNA sequencing of single cells is often cited as a limitation to clinical implementation (Evrony et al., 2021). In part, the cost can be limited by low-coverage sequencing which is sufficient to estimate the density of reads across the genome. Further, it may be possible to minimize the number of cells that are sampled to get a robust estimate of CIN, though sampling too few cells may result in inaccurate measurements. Accordingly, we determined how sampling impacts measurement of mis-segregation rates using approximate Bayesian computation. We first took five random samples from the population of paclitaxel-treated cells each at various sample sizes (Figure 5—figure supplement 5A). We then inferred the mis-segregation rate in each sample and identified the sample size that surpasses an average of 90% accuracy and a low standard error of measurement. We found that even small sample sizes can accurately infer the mis-segregation rate, in this context, with a low standard error (Figure 5—figure supplement 5B-D). A sample size of 60 cells produced the most accurate measurement at 99.5% and a standard error of 0.008 ( ± 0.35 chromosomes). We repeated this analysis using simulated data from the Hybrid selection model and a range of mis-segregation rates spanning what is observed in cancer and non-cancer cultures (Pmisseg ≤ 0.02; see below). We again found a range of sample sizes whose inferred mis-segregation rates underestimate the known value from those simulations (n∈ [20, 40… 180]; Figure 5—figure supplement 5E,F). Across all mis-segregation rates and selective pressures, random samples of 200 cells had a median percent accuracy of 90% and median standard error of 0.0003 ( ± 0.0138 chromosomes per division). The difference in optimal sample sizes between the paclitaxel-treated population and the simulated population is notable and likely due to the presence of ‘clonal’ structures in the simulated population. While the paclitaxel treatment resulted in a uniformly high degree of aneuploidy and little evidence of karyotype selection, the simulated populations after 60 steps (~30 generations) have discrete copy number clusters that may not be captured in each random sample. To verify this, we repeated the analysis using only data from the first time step, prior to the onset of karyotype selection (Figure 5—figure supplement 5H). In this case, we found that the sample size needed to achieve a median 90% accuracy over all simulations in this context is 100 cells, at which point the standard error for Pmisseg is 0.0068 (placing measures within ±0.31 chromosomes per division; Figure 5—figure supplement 5I, J). Thus, a larger number of cells is required in the context of long-term karyotype selection than a more acute time scale, such as we see with paclitaxel.

In conclusion, we recommend using 200 cells from a single sampled site which, at biologically relevant time scales and rates of mis-segregation, provides ≥90% accuracy. These data represent, to our knowledge, the first analysis of how sample size for single-cell sequencing affects the accuracy and measurement of chromosome mis-segregation rates.

Inferring chromosome mis-segregation rates in tumors and organoids

To determine if this framework is clinically applicable, we employed previously published scDNAseq datasets derived from tumor samples and patient-derived organoids (PDO) (Bolhaqueiro et al., 2019; Navin et al., 2011). Importantly, the data from Bolhaqueiro et al. include sample-matched live cell imaging data in colorectal cancer PDOs, with direct observation of chromosome mis-segregation events to compare with inferred measures. We established our panel of measurements on these populations (Figure 6A) and used these to tune the prior distribution of time steps and the rejection threshold for ABC. In sensitivity analysis, 20 steps or greater was sufficient to establish stable estimates of Pmisseg and selection, S (Figure 6—figure supplement 1A-B)—we chose a window of 40–80 steps for further analysis. For rejection thresholds 0.05 and smaller, the inferred mis-segregation rates remained steady (Figure 6—figure supplement 1C). With these model parameters chosen, we evaluated the different selection models, and found that the Abundance model resulted in simulated data that best resembled experimental data, for both exponential and constant-population dynamics (Table 3). Given that the Abundance model is the most biologically relevant, we will use data simulated under this model in our prior dataset for inference.

Figure 6 with 5 supplements see all
Inferring chromosome mis-segregation rates in tumors and organoids Bolhaqueiro et al., 2019Navin et al., 2011.

(A) Computed population summary statistics for colorectal cancer (CRC) patient-derived organoids (PDOs) and breast biopsy scDNAseq datasets from Bolhaqueiro et al., 2019 (gold) and Navin et al., 2011 (pink). (B) Dimensionality reduction analysis of population summary statistics showing biological observations overlaid on, and found within, the space of simulated observations. Point colors show the simulation parameters and summary statistics for all simulations using the following priors and parameters: Growth Model = ‘exponential pseudo-Moran’, Selection Model = ‘Abundance’, initial ploidy = 2, time steps ∈[40, 41… 80], S ∈[0,2… 100], Pmisseg∈[0,0.001… 0.050] and a tolerance threshold of 0.05 to reject dissimilar simulation results. (see Materials and Methods). (C) 2D density plots showing joint posterior distributions of Pmisseg and S values from the approximate Bayesian computation analysis of samples 26 N (left) and 24Tb (right) from Bolhaqueiro et al., 2019. White + represents the mean of inferred values. (D) Inferred selective pressures and mis-segregation rates from each scDNAseq dataset (mean and SEM of accepted values). (E) Predicted mis-segregation rates in CRC PDOs and a breast biopsy plotted with approximated mis-segregation rates observed in cancer (blue triangle) and non-cancer (red circle) models (primarily cell lines) from previous studies (Table 5; see Materials and methods). The predicted mis-segregation rates in these cancer-derived samples fall within those observed in cancer cell lines and above those of non-cancer cell lines. (F) Pearson correlation of predicted mis-segregation rates and predicted selective pressures in CRC PDOs from Bolhaqueiro et al., 2019. (G) Pearson correlation of predicted mis-segregation rates and the incidence of observed segregation errors in CRC PDOs from Bolhaqueiro et al., 2019. Error bars represent SEM values. (H) Pearson correlation of observed incidence of segregation errors in CRC PDOs from Bolhaqueiro et al., 2019 to the ploidy-corrected prediction of the observed incidence of segregation errors. These values assume the involvement of 1 chromosome per observed error and are calculated as the (predicted mis-segregation rate) x (mean number of chromosomes observed per cell) x 100. Dotted line = 1:1 reference.

Table 3
Model selection.
SampleGrowt ModelSelectio ModelPPBF (Ho Neutral)PmissegSSteps
7Texponential pseudo-MoranAbundance0.621Inf0.0033 ± 1e-0560.5416 ± 0.205359.8475 ± 0.0937
7Texponential pseudo-MoranDriver0.14Inf0.001 ± 1e-0549.6557 ± 0.238958.7002 ± 0.0943
7Texponential pseudo-MoranHybrid0.239Inf8e-04 ± 1e-0549.3428 ± 0.237758.5789 ± 0.0935
7Texponential pseudo-MoranNeutral0NA9e-04 ± 5e-050 ± 057.7994 ± 0.6728
7Tconstant Wright-FisherAbundance0.985Inf0.0062 ± 2e-0569.7026 ± 0.172459.9318 ± 0.0937
7Tconstant Wright-FisherDriver0NA0.0012 ± 1e-0548.2881 ± 0.238457.5239 ± 0.0933
7Tconstant Wright-FisherHybrid0.015Inf9e-04 ± 1e-0550.7803 ± 0.235958.2514 ± 0.0941
7Tconstant Wright-FisherNeutral0NA9e-04 ± 5e-050 ± 058.7803 ± 0.6701
U1Texponential pseudo-MoranAbundance0.5821999e-04 ± 1e-0556.8672 ± 0.216859.9906 ± 0.0937
U1Texponential pseudo-MoranDriver0.113390.001 ± 1e-0549.6611 ± 0.238958.6886 ± 0.0944
U1Texponential pseudo-MoranHybrid0.156548e-04 ± 1e-0549.3658 ± 0.237558.569 ± 0.0935
U1Texponential pseudo-MoranNeutral0.14919e-04 ± 5e-050 ± 057.7102 ± 0.67
U1Tconstant Wright-FisherAbundance0.6542900.001 ± 1e-0561.4358 ± 0.202960.0021 ± 0.0937
U1Tconstant Wright-FisherDriver0.115510.0012 ± 1e-0548.2767 ± 0.238357.5267 ± 0.0934
U1Tconstant Wright-FisherHybrid0.115519e-04 ± 1e-0550.8033 ± 0.235858.2507 ± 0.0941
U1Tconstant Wright-FisherNeutral0.11519e-04 ± 5e-050 ± 058.7803 ± 0.6701
U2Texponential pseudo-MoranAbundance0.6282510.0054 ± 1e-0559.4269 ± 0.210859.8349 ± 0.0935
U2Texponential pseudo-MoranDriver0.079320.0027 ± 2e-0550.1513 ± 0.239657.4538 ± 0.0934
U2Texponential pseudo-MoranHybrid0.166660.0022 ± 2e-0548.7779 ± 0.241357.7078 ± 0.0934
U2Texponential pseudo-MoranNeutral0.12710.0021 ± 7e-050 ± 056.8535 ± 0.6619
U2Tconstant Wright-FisherAbundance0.91828170.0112 ± 3e-0569.7222 ± 0.170360.0655 ± 0.0934
U2Tconstant Wright-FisherDriver0.00140.0027 ± 2e-0548.7794 ± 0.238956.4812 ± 0.0919
U2Tconstant Wright-FisherHybrid0.0641960.0022 ± 1e-0550.9564 ± 0.237957.1161 ± 0.0925
U2Tconstant Wright-FisherNeutral0.01710.0022 ± 1e-040 ± 057.7898 ± 0.6841
U3Texponential pseudo-MoranAbundance0.5821990.0029 ± 1e-0560.9557 ± 0.209159.8273 ± 0.0938
U3Texponential pseudo-MoranDriver0.113390.001 ± 1e-0549.6707 ± 0.238958.6986 ± 0.0944
U3Texponential pseudo-MoranHybrid0.156548e-04 ± 1e-0549.3754 ± 0.237658.5711 ± 0.0935
U3Texponential pseudo-MoranNeutral0.14919e-04 ± 5e-050 ± 057.7102 ± 0.67
U3Tconstant Wright-FisherAbundance0.736Inf0.0052 ± 2e-0569.8357 ± 0.171359.932 ± 0.0934
U3Tconstant Wright-FisherDriver0.13Inf0.0012 ± 1e-0548.2864 ± 0.238357.5385 ± 0.0934
U3Tconstant Wright-FisherHybrid0.134Inf9e-04 ± 1e-0550.8219 ± 0.235758.2482 ± 0.0941
U3Tconstant Wright-FisherNeutral0NA9e-04 ± 5e-050 ± 058.8567 ± 0.6676
14Texponential pseudo-MoranAbundance0.5821999e-04 ± 1e-0556.8672 ± 0.216859.9906 ± 0.0937
14Texponential pseudo-MoranDriver0.113390.001 ± 1e-0549.6614 ± 0.23958.695 ± 0.0944
14Texponential pseudo-MoranHybrid0.156548e-04 ± 1e-0549.3716 ± 0.237558.5632 ± 0.0935
14Texponential pseudo-MoranNeutral0.14919e-04 ± 5e-050 ± 057.7102 ± 0.67
14Tconstant Wright-FisherAbundance0.6542900.0011 ± 1e-0562.8579 ± 0.207560.0029 ± 0.0936
14Tconstant Wright-FisherDriver0.115510.0012 ± 1e-0548.2967 ± 0.238357.5295 ± 0.0934
14Tconstant Wright-FisherHybrid0.115519e-04 ± 1e-0550.8274 ± 0.235758.2478 ± 0.0941
14Tconstant Wright-FisherNeutral0.11519e-04 ± 5e-050 ± 058.8567 ± 0.6676
16Texponential pseudo-MoranAbundance0.5821990.002 ± 1e-0561.2401 ± 0.202859.9109 ± 0.0935
16Texponential pseudo-MoranDriver0.113390.001 ± 1e-0549.6539 ± 0.238958.7006 ± 0.0943
16Texponential pseudo-MoranHybrid0.156548e-04 ± 1e-0549.3611 ± 0.237658.574 ± 0.0935
16Texponential pseudo-MoranNeutral0.14919e-04 ± 5e-050 ± 057.7994 ± 0.6728
16Tconstant Wright-FisherAbundance0.6542900.0038 ± 1e-0569.8456 ± 0.170159.9523 ± 0.0936
16Tconstant Wright-FisherDriver0.115510.0012 ± 1e-0548.261 ± 0.238457.5233 ± 0.0933
16Tconstant Wright-FisherHybrid0.115519e-04 ± 1e-0550.7713 ± 0.235958.2554 ± 0.0941
16Tconstant Wright-FisherNeutral0.11519e-04 ± 5e-050 ± 058.7803 ± 0.6701
19Taexponential pseudo-MoranAbundance0.7113130.004 ± 1e-0560.6391 ± 0.207459.7801 ± 0.0934
19Taexponential pseudo-MoranDriver0.038170.0028 ± 2e-0550.2185 ± 0.239957.3764 ± 0.0934
19Taexponential pseudo-MoranHybrid0.135590.0022 ± 3e-0548.3823 ± 0.24257.5368 ± 0.0935
19Taexponential pseudo-MoranNeutral0.11610.0022 ± 9e-050 ± 056.5955 ± 0.6549
19Taconstant Wright-FisherAbundance0.97117600.0075 ± 2e-0569.3863 ± 0.173559.956 ± 0.0938
19Taconstant Wright-FisherDriver000.0028 ± 2e-0548.8413 ± 0.239256.4529 ± 0.0917
19Taconstant Wright-FisherHybrid0.0263150.0023 ± 1e-0550.8588 ± 0.238357.1031 ± 0.0925
19Taconstant Wright-FisherNeutral0.00410.0023 ± 1e-040 ± 057.9522 ± 0.6869
19Tbexponential pseudo-MoranAbundance0.7273200.0036 ± 1e-0560.5885 ± 0.208559.829 ± 0.0938
19Tbexponential pseudo-MoranDriver0.03130.001 ± 1e-0549.6622 ± 0.238958.6929 ± 0.0944
19Tbexponential pseudo-MoranHybrid0.127568e-04 ± 1e-0548.5237 ± 0.232258.9663 ± 0.0931
19Tbexponential pseudo-MoranNeutral0.11619e-04 ± 5e-050 ± 057.7102 ± 0.67
19Tbconstant Wright-FisherAbundance0.979473200.0068 ± 2e-0569.5697 ± 0.17359.9232 ± 0.0935
19Tbconstant Wright-FisherDriver000.0012 ± 1e-0548.2786 ± 0.238357.5433 ± 0.0934
19Tbconstant Wright-FisherHybrid0.029829e-04 ± 1e-0550.8162 ± 0.235758.2495 ± 0.0941
19Tbconstant Wright-FisherNeutral0.00119e-04 ± 5e-050 ± 058.8376 ± 0.669
24Taexponential pseudo-MoranAbundance0.7313210.0036 ± 1e-0560.5303 ± 0.208259.8208 ± 0.0938
24Taexponential pseudo-MoranDriver0.029130.001 ± 1e-0549.6703 ± 0.238958.6938 ± 0.0944
24Taexponential pseudo-MoranHybrid0.125558e-04 ± 1e-0549.3669 ± 0.237658.5778 ± 0.0935
24Taexponential pseudo-MoranNeutral0.11619e-04 ± 5e-050 ± 057.7102 ± 0.67
24Taconstant Wright-FisherAbundance0.979473460.0068 ± 2e-0569.6173 ± 0.17359.933 ± 0.0934
24Taconstant Wright-FisherDriver000.0012 ± 1e-0548.2789 ± 0.238357.5377 ± 0.0934
24Taconstant Wright-FisherHybrid0.029569e-04 ± 1e-0550.8229 ± 0.235758.2524 ± 0.0941
24Taconstant Wright-FisherNeutral0.00119e-04 ± 5e-050 ± 058.8567 ± 0.6676
24Tbexponential pseudo-MoranAbundance0.682940.0046 ± 1e-0560.2602 ± 0.208459.8073 ± 0.0936
24Tbexponential pseudo-MoranDriver0.054230.0031 ± 3e-0550.2981 ± 0.239957.2927 ± 0.0934
24Tbexponential pseudo-MoranHybrid0.149650.0025 ± 4e-0548.3833 ± 0.24457.4236 ± 0.0936
24Tbexponential pseudo-MoranNeutral0.11810.0025 ± 0.000130 ± 056.7229 ± 0.6579
24Tbconstant Wright-FisherAbundance0.95477300.0215 ± 0.0001133.6703 ± 0.296259.9064 ± 0.0937
24Tbconstant Wright-FisherDriver020.003 ± 2e-0548.7528 ± 0.239356.4175 ± 0.0918
24Tbconstant Wright-FisherHybrid0.0393180.0024 ± 2e-0550.7006 ± 0.238957.107 ± 0.0925
24Tbconstant Wright-FisherNeutral0.00610.0024 ± 0.000110 ± 058.0318 ± 0.6822
26Nexponential pseudo-MoranAbundance0.5821990.0021 ± 1e-0560.9877 ± 0.203159.9205 ± 0.0934
26Nexponential pseudo-MoranDriver0.113390.001 ± 1e-0549.6389 ± 0.238958.7018 ± 0.0944
26Nexponential pseudo-MoranHybrid0.156548e-04 ± 1e-0549.3389 ± 0.237758.5755 ± 0.0935
26Nexponential pseudo-MoranNeutral0.14919e-04 ± 5e-050 ± 057.7994 ± 0.6728
26Nconstant Wright-FisherAbundance0.6542900.0039 ± 1e-0569.794 ± 0.170459.9547 ± 0.0935
26Nconstant Wright-FisherDriver0.115510.0012 ± 1e-0548.2849 ± 0.238457.5175 ± 0.0933
26Nconstant Wright-FisherHybrid0.115519e-04 ± 1e-0550.737 ± 0.235958.2609 ± 0.0941
26Nconstant Wright-FisherNeutral0.11519e-04 ± 5e-050 ± 058.7803 ± 0.6701
9Texponential pseudo-MoranAbundance0.6852990.0044 ± 1e-0560.2829 ± 0.208659.7955 ± 0.0936
9Texponential pseudo-MoranDriver0.052230.0029 ± 2e-0550.2323 ± 0.239857.3657 ± 0.0934
9Texponential pseudo-MoranHybrid0.147640.0022 ± 3e-0548.3829 ± 0.242257.5193 ± 0.0936
9Texponential pseudo-MoranNeutral0.11710.0023 ± 9e-050 ± 056.6083 ± 0.6581
9Tconstant Wright-FisherAbundance0.95892990.0087 ± 2e-0569.6836 ± 0.172459.926 ± 0.0937
9Tconstant Wright-FisherDriver010.0028 ± 2e-0548.8394 ± 0.239256.4465 ± 0.0917
9Tconstant Wright-FisherHybrid0.0373600.0023 ± 1e-0550.8477 ± 0.238457.0952 ± 0.0925
9Tconstant Wright-FisherNeutral0.00510.0023 ± 1e-040 ± 057.9427 ± 0.687
PolyB1exponential pseudo-MoranAbundance0.6352610.0053 ± 1e-0559.5088 ± 0.210459.8379 ± 0.0935
PolyB1exponential pseudo-MoranDriver0.076310.0028 ± 2e-0550.2364 ± 0.239857.4025 ± 0.0934
PolyB1exponential pseudo-MoranHybrid0.164670.0022 ± 3e-0548.6949 ± 0.241957.6322 ± 0.0934
PolyB1exponential pseudo-MoranNeutral0.12410.0022 ± 9e-050 ± 056.5955 ± 0.6549
PolyB1constant Wright-FisherAbundance0.92534820.0111 ± 3e-0570.2557 ± 0.16960.042 ± 0.0936
PolyB1constant Wright-FisherDriver0.00140.0028 ± 2e-0548.8194 ± 0.239156.4451 ± 0.0917
PolyB1constant Wright-FisherHybrid0.0612280.0023 ± 1e-0550.895 ± 0.238157.1073 ± 0.0925
PolyB1constant Wright-FisherNeutral0.01410.0023 ± 1e-040 ± 057.9809 ± 0.6861
PolyB2exponential pseudo-MoranAbundance0.6032180.0059 ± 1e-0558.6612 ± 0.21259.7835 ± 0.0937
PolyB2exponential pseudo-MoranDriver0.086310.0038 ± 4e-0550.2948 ± 0.239457.0217 ± 0.093
PolyB2exponential pseudo-MoranHybrid0.17610.004 ± 7e-0548.9466 ± 0.247257.28 ± 0.0942
PolyB2exponential pseudo-MoranNeutral0.14110.0033 ± 0.000220 ± 056.5732 ± 0.6597
PolyB2constant Wright-FisherAbundance0.89312770.0301 ± 1e-043.0543 ± 0.016559.9142 ± 0.0936
PolyB2constant Wright-FisherDriver0.00340.0034 ± 3e-0548.7328 ± 0.239656.3664 ± 0.0917
PolyB2constant Wright-FisherHybrid0.069980.0027 ± 2e-0550.3534 ± 0.240557.1445 ± 0.0928
PolyB2constant Wright-FisherNeutral0.03610.0026 ± 0.000140 ± 058.1592 ± 0.6741

Having confirmed the summary statistics from these samples were within the space of the simulation data with our chosen priors (Figure 6B), we performed ABC analysis on these datasets to infer rates of chromosome mis-segregation and levels of selection pressure and display the joint posterior distributions as 2D density plots (Figure 6C and D; Figure 6—figure supplements 2 and 3). Figure 6C illustrates the results for two individual colon organoid lines, showing the distribution of parameters used for simulations that gave the most similar results. With ABC, inferred parameters fall within rates of mis-segregation of about 0.001–0.006. Applied to a near-diploid cell, this translates to a range of about 5–38% of cell divisions having one chromosome mis-segregation. Importantly, these inferred rates of chromosome mis-segregation fall within the range of approximated per chromosome rates experimentally observed in cancer cell lines and human tumors (Figure 6E;Table 4, Table 5; Bakhoum et al., 2014; Bakhoum et al., 2011; Bakhoum et al., 2009; Dewhurst et al., 2014; Nicholson et al., 2015; Orr et al., 2016; Thompson and Compton, 2008; Worrall et al., 2018; Zasadil et al., 2014). Higher inferred mis-segregation rates tended to coincide with lower inferred selection experienced in these samples (Figure 6F). Posterior distributions in these samples were skewed toward high selection (S) indicating the presence stabilizing selection in all cases, where the average of the distributions of some samples were slightly lower or higher (Figure 6—figure supplement 3).

Table 4
Model selection with selective pressure constrained to S = 1.
SampleGrowth ModelSelection ModelPPBF (Ho Neutral)PmissegSSteps
7Texponential pseudo-MoranAbundance0.27419e-04 ± 5e-051 ± 058.2452 ± 0.6646
7Texponential pseudo-MoranDriver0.23819e-04 ± 5e-051 ± 058.4745 ± 0.6725
7Texponential pseudo-MoranHybrid0.2619e-04 ± 5e-051 ± 058.586 ± 0.6668
7Texponential pseudo-MoranNeutral0.22819e-04 ± 6e-051 ± 058.5446 ± 0.6791
7Tconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.8089 ± 0.6627
7Tconstant Wright-FisherDriver0.2419e-04 ± 6e-051 ± 058.1783 ± 0.6771
7Tconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.0924 ± 0.6742
7Tconstant Wright-FisherNeutral0.24519e-04 ± 7e-051 ± 058.7516 ± 0.6787
U1Texponential pseudo-MoranAbundance0.27519e-04 ± 5e-051 ± 058.2452 ± 0.6646
U1Texponential pseudo-MoranDriver0.23919e-04 ± 5e-051 ± 058.4745 ± 0.6725
U1Texponential pseudo-MoranHybrid0.25819e-04 ± 5e-051 ± 058.586 ± 0.6668
U1Texponential pseudo-MoranNeutral0.22819e-04 ± 6e-051 ± 058.5446 ± 0.6791
U1Tconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.8089 ± 0.6627
U1Tconstant Wright-FisherDriver0.2419e-04 ± 6e-051 ± 058.1783 ± 0.6771
U1Tconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.1592 ± 0.6715
U1Tconstant Wright-FisherNeutral0.24519e-04 ± 7e-051 ± 058.7516 ± 0.6787
U2Texponential pseudo-MoranAbundance0.27610.0021 ± 8e-051 ± 057.3057 ± 0.653
U2Texponential pseudo-MoranDriver0.23510.0024 ± 0.000111 ± 057.7452 ± 0.6634
U2Texponential pseudo-MoranHybrid0.26410.0021 ± 7e-051 ± 058.1274 ± 0.654
U2Texponential pseudo-MoranNeutral0.22510.0024 ± 0.000111 ± 057.8758 ± 0.6772
U2Tconstant Wright-FisherAbundance0.26910.0023 ± 1e-041 ± 058.3439 ± 0.6532
U2Tconstant Wright-FisherDriver0.23310.0023 ± 9e-051 ± 057.4777 ± 0.693
U2Tconstant Wright-FisherHybrid0.26310.0023 ± 1e-041 ± 057.8662 ± 0.6683
U2Tconstant Wright-FisherNeutral0.23610.0025 ± 0.000121 ± 057.1433 ± 0.6655
U3Texponential pseudo-MoranAbundance0.27519e-04 ± 5e-051 ± 058.1624 ± 0.6643
U3Texponential pseudo-MoranDriver0.23919e-04 ± 5e-051 ± 058.4554 ± 0.6736
U3Texponential pseudo-MoranHybrid0.25819e-04 ± 5e-051 ± 058.586 ± 0.6668
U3Texponential pseudo-MoranNeutral0.22819e-04 ± 6e-051 ± 058.6178 ± 0.6777
U3Tconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.7611 ± 0.6614
U3Tconstant Wright-FisherDriver0.2419e-04 ± 6e-051 ± 058.1783 ± 0.6771
U3Tconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.0955 ± 0.674
U3Tconstant Wright-FisherNeutral0.24519e-04 ± 7e-051 ± 058.7516 ± 0.6787
14Texponential pseudo-MoranAbundance0.27519e-04 ± 5e-051 ± 058.1624 ± 0.6643
14Texponential pseudo-MoranDriver0.23919e-04 ± 5e-051 ± 058.4554 ± 0.6736
14Texponential pseudo-MoranHybrid0.25819e-04 ± 5e-051 ± 058.586 ± 0.6668
14Texponential pseudo-MoranNeutral0.22819e-04 ± 6e-051 ± 058.5446 ± 0.6791
14Tconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.8089 ± 0.6627
14Tconstant Wright-FisherDriver0.2419e-04 ± 6e-051 ± 058.1783 ± 0.6771
14Tconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.0924 ± 0.6739
14Tconstant Wright-FisherNeutral0.24519e-04 ± 7e-051 ± 058.7516 ± 0.6787
16Texponential pseudo-MoranAbundance0.27419e-04 ± 5e-051 ± 058.2452 ± 0.6646
16Texponential pseudo-MoranDriver0.23819e-04 ± 5e-051 ± 058.4745 ± 0.6725
16Texponential pseudo-MoranHybrid0.2619e-04 ± 5e-051 ± 058.586 ± 0.6668
16Texponential pseudo-MoranNeutral0.22810.001 ± 6e-051 ± 058.6274 ± 0.6789
16Tconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.8089 ± 0.6627
16Tconstant Wright-FisherDriver0.2419e-04 ± 6e-051 ± 058.1783 ± 0.6771
16Tconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.1051 ± 0.6742
16Tconstant Wright-FisherNeutral0.24519e-04 ± 7e-051 ± 058.7516 ± 0.6787
19Taexponential pseudo-MoranAbundance0.27310.0021 ± 8e-051 ± 057.4045 ± 0.6565
19Taexponential pseudo-MoranDriver0.24310.0024 ± 0.000111 ± 057.8025 ± 0.663
19Taexponential pseudo-MoranHybrid0.26110.0022 ± 8e-051 ± 057.9108 ± 0.65
19Taexponential pseudo-MoranNeutral0.22210.0025 ± 0.000121 ± 057.9331 ± 0.6777
19Taconstant Wright-FisherAbundance0.2710.0024 ± 0.000111 ± 058.2866 ± 0.6566
19Taconstant Wright-FisherDriver0.23310.0023 ± 1e-041 ± 057.8185 ± 0.6927
19Taconstant Wright-FisherHybrid0.26110.0023 ± 1e-041 ± 058.0478 ± 0.6705
19Taconstant Wright-FisherNeutral0.23710.0025 ± 0.000121 ± 057.2261 ± 0.6669
19Tbexponential pseudo-MoranAbundance0.27519e-04 ± 5e-051 ± 058.1624 ± 0.6643
19Tbexponential pseudo-MoranDriver0.23919e-04 ± 5e-051 ± 058.4554 ± 0.6736
19Tbexponential pseudo-MoranHybrid0.25819e-04 ± 5e-051 ± 058.586 ± 0.6668
19Tbexponential pseudo-MoranNeutral0.22819e-04 ± 6e-051 ± 058.5796 ± 0.6796
19Tbconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.7611 ± 0.6614
19Tbconstant Wright-FisherDriver0.2419e-04 ± 6e-051 ± 058.1178 ± 0.679
19Tbconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.1592 ± 0.6715
19Tbconstant Wright-FisherNeutral0.24519e-04 ± 7e-051 ± 058.7516 ± 0.6787
24Taexponential pseudo-MoranAbundance0.27519e-04 ± 5e-051 ± 058.1624 ± 0.6643
24Taexponential pseudo-MoranDriver0.23919e-04 ± 5e-051 ± 058.4554 ± 0.6736
24Taexponential pseudo-MoranHybrid0.25819e-04 ± 5e-051 ± 058.586 ± 0.6668
24Taexponential pseudo-MoranNeutral0.22819e-04 ± 6e-051 ± 058.6656 ± 0.6783
24Taconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.7611 ± 0.6614
24Taconstant Wright-FisherDriver0.2419e-04 ± 6e-051 ± 058.1783 ± 0.6771
24Taconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.1592 ± 0.6715
24Taconstant Wright-FisherNeutral0.24519e-04 ± 7e-051 ± 058.7516 ± 0.6787
24Tbexponential pseudo-MoranAbundance0.27310.0023 ± 0.000111 ± 057.0446 ± 0.6526
24Tbexponential pseudo-MoranDriver0.24210.0025 ± 0.000121 ± 057.551 ± 0.6661
24Tbexponential pseudo-MoranHybrid0.26410.0022 ± 9e-051 ± 057.9108 ± 0.6512
24Tbexponential pseudo-MoranNeutral0.22210.0026 ± 0.000131 ± 057.7516 ± 0.6758
24Tbconstant Wright-FisherAbundance0.26710.0024 ± 0.000131 ± 058.379 ± 0.6601
24Tbconstant Wright-FisherDriver0.23710.0024 ± 1e-041 ± 057.7357 ± 0.6922
24Tbconstant Wright-FisherHybrid0.25710.0023 ± 1e-041 ± 057.9045 ± 0.6718
24Tbconstant Wright-FisherNeutral0.23910.0025 ± 0.000121 ± 057.2643 ± 0.6726
26Nexponential pseudo-MoranAbundance0.27419e-04 ± 5e-051 ± 058.2452 ± 0.6646
26Nexponential pseudo-MoranDriver0.23919e-04 ± 5e-051 ± 058.4045 ± 0.6706
26Nexponential pseudo-MoranHybrid0.2619e-04 ± 5e-051 ± 058.586 ± 0.6668
26Nexponential pseudo-MoranNeutral0.22710.001 ± 7e-051 ± 058.6815 ± 0.6776
26Nconstant Wright-FisherAbundance0.25919e-04 ± 6e-051 ± 058.8089 ± 0.6627
26Nconstant Wright-FisherDriver0.23919e-04 ± 6e-051 ± 058.1783 ± 0.6771
26Nconstant Wright-FisherHybrid0.25719e-04 ± 5e-051 ± 059.1178 ± 0.6745
26Nconstant Wright-FisherNeutral0.24510.001 ± 7e-051 ± 058.6879 ± 0.6762
9Texponential pseudo-MoranAbundance0.27410.0021 ± 8e-051 ± 057.3854 ± 0.6574
9Texponential pseudo-MoranDriver0.24210.0024 ± 0.000111 ± 057.8025 ± 0.663
9Texponential pseudo-MoranHybrid0.26110.0022 ± 8e-051 ± 057.9108 ± 0.65
9Texponential pseudo-MoranNeutral0.22210.0025 ± 0.000121 ± 057.9522 ± 0.6787
9Tconstant Wright-FisherAbundance0.26910.0024 ± 0.000111 ± 058.2866 ± 0.6566
9Tconstant Wright-FisherDriver0.23310.0023 ± 1e-041 ± 057.9076 ± 0.6927
9Tconstant Wright-FisherHybrid0.26110.0023 ± 1e-041 ± 058.1115 ± 0.6708
9Tconstant Wright-FisherNeutral0.23610.0025 ± 0.000121 ± 057.2261 ± 0.6669
PolyB1exponential pseudo-MoranAbundance0.27410.0021 ± 8e-051 ± 057.4045 ± 0.6565
PolyB1exponential pseudo-MoranDriver0.24310.0024 ± 0.000111 ± 057.7102 ± 0.6622
PolyB1exponential pseudo-MoranHybrid0.26110.0022 ± 8e-051 ± 057.9459 ± 0.6512
PolyB1exponential pseudo-MoranNeutral0.22210.0025 ± 0.000111 ± 057.9522 ± 0.6776
PolyB1constant Wright-FisherAbundance0.27110.0023 ± 0.000111 ± 058.2834 ± 0.6575
PolyB1constant Wright-FisherDriver0.23110.0023 ± 9e-051 ± 057.6656 ± 0.6949
PolyB1constant Wright-FisherHybrid0.26110.0023 ± 1e-041 ± 057.9713 ± 0.6668
PolyB1constant Wright-FisherNeutral0.23710.0025 ± 0.000121 ± 057.207 ± 0.6674
PolyB2exponential pseudo-MoranAbundance0.27210.0027 ± 2e-041 ± 056.8471 ± 0.6544
PolyB2exponential pseudo-MoranDriver0.24510.0029 ± 0.000211 ± 057.3312 ± 0.6609
PolyB2exponential pseudo-MoranHybrid0.26310.0024 ± 0.000111 ± 057.9204 ± 0.6466
PolyB2exponential pseudo-MoranNeutral0.22110.0029 ± 0.000171 ± 057.4236 ± 0.6784
PolyB2constant Wright-FisherAbundance0.26810.0025 ± 0.000131 ± 058.2484 ± 0.6616
PolyB2constant Wright-FisherDriver0.23510.0026 ± 0.000141 ± 057.5796 ± 0.6897
PolyB2constant Wright-FisherHybrid0.25710.0026 ± 0.000151 ± 058.1115 ± 0.6741
PolyB2constant Wright-FisherNeutral0.2410.0027 ± 0.000141 ± 057.379 ± 0.6701
Table 5
Approximate reported per chromosome mis-segregation rates.
1st AuthorDOIModelTumor?StatisticAssessmentApproximate observed frequency %Aprrox modal chromosome # (ATCC)Approximate mis-segregation rate (per chromosome)
Bakhoumhttps://doi.org/10.1158/1078-0432.CCR-11-2049Tumor-TMATumorReportedLagging/Bridging31.3460.00680
Orrhttps://doi.org/10.1016/j.celrep.2016.10.030U2OSTumorApprox. MeanLagging32.5460.00707
Orrhttps://doi.org/10.1016/j.celrep.2016.10.030HeLaTumorApprox. MeanLagging22820.00268
Orrhttps://doi.org/10.1016/j.celrep.2016.10.030SW-620TumorApprox. MeanLagging22.5500.00450
Orrhttps://doi.org/10.1016/j.celrep.2016.10.030RPE1Non-tumorApprox. MeanLagging2.5460.00054
Orrhttps://doi.org/10.1016/j.celrep.2016.10.030BJNon-tumorApprox. MeanLagging8460.00174
Nicholsonhttps://doi.org/10.7554/eLife.05068AmniocyteNon-tumorApprox. MeanLagging0460.00000
Nicholsonhttps://doi.org/10.7554/eLife.05068DLD1TumorApprox. MeanLagging1460.00022
Dewhursthttps://doi.org/10.1158/2159-8290.CD-13-0285HCT116-DiploidTumorApprox. MeanLagging/Bridging23450.00511
Dewhursthttps://doi.org/10.1158/2159-8290.CD-13-0285HCT116-TetraploidTumorApprox. MeanLagging/Bridging50900.00556
Bakhoumhttps://doi.org/10.1038/ncb1809U2OSTumorReportedLagging460.01000
Zasadilhttps://doi.org/10.1126/scitranslmed.3007965CAL51TumorApprox. MeanLagging0.5440.00011
Thompsonhttps://doi.org/10.1083/jcb.200712029RPE1Non-tumorApprox. MeanAcute aneuploidy via FISH460.00025
Thompsonhttps://doi.org/10.1083/jcb.200712029HCT116-DiploidTumorApprox. MeanAcute aneuploidy via FISH450.00025
Thompsonhttps://doi.org/10.1083/jcb.200712029HT29TumorApprox. MeanAcute aneuploidy via FISH710.00250
Thompsonhttps://doi.org/10.1083/jcb.200712029Caco2TumorApprox. MeanAcute aneuploidy via FISH960.00900
Thompsonhttps://doi.org/10.1083/jcb.200712029MCF-7TumorApprox. MeanAcute aneuploidy via FISH820.00700
Bakhoumhttps://doi.org/10.1016/j.cub.2014.01.019HCT116-DiploidTumorApprox. MeanLagging6450.00133
Bakhoumhttps://doi.org/10.1016/j.cub.2014.01.019DLD1TumorApprox. MeanLagging2460.00043
Bakhoumhttps://doi.org/10.1016/j.cub.2014.01.019HT29TumorApprox. MeanLagging14710.00197
Bakhoumhttps://doi.org/10.1016/j.cub.2014.01.019SW-620TumorApprox. MeanLagging12500.00240
Bakhoumhttps://doi.org/10.1016/j.cub.2014.01.019MCF-7TumorApprox. MeanLagging17820.00207
Bakhoumhttps://doi.org/10.1016/j.cub.2014.01.019HeLaTumorApprox. MeanLagging13820.00159
Worrallhttps://doi.org/10.1016/j.celrep.2018.05.047BJNon-tumorApprox. MeanUnspecified Error5460.00109
Worrallhttps://doi.org/10.1016/j.celrep.2018.05.047RPE1Non-tumorApprox. MeanUnspecified Error5460.00109

To confirm the relevance of the inferred scalar exponent we performed our model selection scheme using only the simulation data with unmodified fitness values (S = 1; Table 4). In this case, we found that the inferred mis-segregation rates for most samples fell well below the expected range found in cancer cell lines (Figure 6E). Additionally, when we inferred mis-segregation rates and selection in the early timepoint of longitudinally sequenced organoid clones from Bolhaqueiro et al., 2019, the composition of the resultant populations simulated using these inferred characteristics better resembled the late-timepoint organoid data than those with unmodified selection values (S = 1; Figure 6—figure supplements 4 and 5).

As further validation for mis-segregation rates, we compared these inferred rates from CRC PDOs with those directly measured in live imaging from Bolhaqueiro et al., 2019. Although mis-segregation cannot be directly inferred from microscopy, diversity should correlate with the observed rate of mitotic errors. There was a strong correlation but for two outliers—14T and U1T (Figure 6G). In fact, when adjusting to the same scale and correcting for cell ploidy, these data follow a strong positive linear trend with a slightly lower slope than a 1:1 correlation, which could reflect an overestimation of mis-segregation rates in the microscopy data (Figure 6H). Particularly with lagging chromosomes, despite a chromosome’s involvement in an observed segregation defect, it may end up in the correct daughter cell. Overall, these results indicate that the inferred measures using approximate Bayesian computation and scDNAseq account for selection and provide a quantitative measure of CIN.

Discussion

The clinical assessment of mutations, short indels, and microsatellite instability in human cancer determined by short-read sequencing currently guide clinical care. By contrast, CIN is highly prevalent, yet has remained largely intractable to clinical measures. Single-cell DNA sequencing now promises detailed karyotypic analysis across hundreds of cells, yet selective pressure suppresses the observed karyotype heterogeneity within a tumor. Optimal clinical measurement of CIN may be achieved with scDNAseq, but must additionally account for selective pressure, which reduces karyotype heterogeneity.

Despite the major limitations with current measures of CIN, emerging evidence hints at its utility as a biomarker to predict benefit to cancer therapy. For example, CIN measures appear to predict therapeutic response to paclitaxel (Janssen et al., 2009; Scribano et al., 2021; Swanton et al., 2009). Nevertheless, existing measures of CIN have had significant limitations. FISH and histological analysis of mitotic abnormalities are limited in quantifying specific chromosomes or requiring highly proliferative tumor types, such as lymphomas and leukemia. Gene expression profiles are proposed to correlate with CIN among populations of tumor samples (Carter et al., 2006), although they happen to correlate better with tumor proliferation (Sheltzer, 2013); in any case, they are correlations across populations of tumors, not suitable as an individualized diagnostic. We conclude that scDNAseq is the most complete and tractable measure of cellular karyotypes, and sampling at least 200 cells, coupled with computational models and ABC, promises to offer the best measure of tumor CIN.

Computational modeling of aneuploidy and CIN has been used to explore evolution in the context of numerical CIN and karyotype selection (Elizalde et al., 2018; Gao et al., 2016; Gusev et al., 2001; Gusev et al., 2000; Laughney et al., 2015; Nowak et al., 2002). Gusev and Nowak lay the foundation for mathematical modeling of CIN. While Gusev focused on the karyotypic outcomes of CIN, Nowak considered the effects of CIN-inducing mutations and the subsequent rate of LOH. Neither considered the individual fitness differences between specific karyotypes (Gusev et al., 2001; Gusev et al., 2000; Nowak et al., 2002). This was improved in Laughney et al., 2015 and Elizalde et al., 2018 where the authors leveraged the chromosome scores derived in Davoli et al., 2013, which enable the inclusion of oncogenes and tumor suppressors in models of CIN as we have done. These studies have provided important insights such as the role of whole-genome doubling as an evolutionary bridge to optimized chromosome stoichiometry. Yet the populations derived in these studies tend to vary to a greater degree than observed with scDNAseq, as they do not model strong selection against aneuploidy. Further, they do not attempt to use their models to measure CIN in biological samples. Here, we build on these models by considering, in addition to the selection on driver genes, the stabilizing selection wrought by chromosomal gene abundance. Further, we consider that the magnitude of selection pressure may not be a constant and implement a modifier to tune selection in our models. Lastly, we use our models as a quantitative measure of CIN that accounts for this selection.

Previous studies using single-cell sequencing identified surprisingly low karyotypic variance in human tumors including breast cancer (Gao et al., 2016; Kim et al., 2018; Wang et al., 2014) and colorectal and ovarian cancer organoids (Bolhaqueiro et al., 2019; Nelson et al., 2020). It has been difficult to understand these findings in the light of widespread CIN in human cancer (Sheltzer and Amon, 2011; Silk et al., 2013; Vasudevan et al., 2020; Weaver et al., 2007; Weaver and Cleveland, 2009). The best explanation of this apparent paradox is selection, which moderates karyotypic variance. Accounting for this, we can infer rates of chromosome mis-segregation in tumors or PDOs well within the range of rates observed microscopically in cancer cell lines. Additionally, no previous work, to our knowledge, has estimated the required sample size to infer CIN from scDNAseq data.

As described by others (Dewhurst et al., 2014; López et al., 2020), and consistent with our findings, early emergence of polyploid cells can markedly reduce apparent selection, leading to an elevated karyotype diversity over time. While we do not explicitly induce chance of whole genome doubling (WGD) events in simulations, populations that begin either diploid or tetraploid converge on near-triploid karyotypes over time, consistent with the notion that WGD can act as an evolutionary bridge to highly aneuploid karyotypes. Notably, our analysis indicates the samples with apparent polyploidy experienced among the lowest levels of karyotype selection.

In some early studies, CIN is considered a binary process—present or absent. We assumed that CIN measures are scalar, not binary, and measure this by rate of chromosome mis-segregation per division. A scalar is appropriate if, for example, there was a consistent probability of chromosome mis-segregation per division. However, we recognize that some mechanisms may not well adhere to this simplified model of CIN. For example, tumors with centrosome amplification may at times undergo bipolar division without mis-segregation, or, at other times, a multipolar division with extensive mis-segregation. Further, it is possible that some mechanisms may have correlated mis-segregations such that a daughter cell that gains one chromosome is more likely to gain other chromosomes, rather than lose them. Another possibility is that CIN could result in the mis-regulation of genes that further modify the rate of CIN. Our model does not yet account for punctuated behavior or changing rates of CIN. Furthermore, while recent studies have reported non-random mis-segregation of chromosomes (Dumont et al., 2020; Worrall et al., 2018), we did not incorporate these biases into our models as these studies do not reach consensus on which chromosomes are more frequently mis-segregated, which may be model-dependent.

Our approach reconstructs phylogenetic trees via copy number variation (CNV) analysis. This approach may be suboptimal given the selection on aneuploid states, and could be particularly problematic in the setting of convergent evolution. It is possible that this method results in low accuracy of the reconstructed phylogenies. Alternative approaches are possible, but would likely require re-design of the scDNAseq assay to include spiked-in primers that span highly polymorphic regions on each chromosome. If this were done, these sequences could be read in all cells and single-nucleotide polymorphisms could track individual maternal and paternal chromosomes, allowing a means of reconstructing cell phylogeny independent of CNVs. Despite this limitation, our phylogenetic reconstructions did seem to allow inference of CIN measures consistent with directly observed rates of chromosome mis-segregation in our taxol-induced CIN model as well as several independent cancer PDO models and cell lines.

A final limitation of our approach is we used previous estimates of cellular selection in our agent-based model and used these selection models to infer quantitative measures of CIN. While this approach seems to perform well in estimates of mis-segregation rates, we recognize that the selection models do not necessarily represent the real selective pressures on distinct aneuploidies. Future investigations are necessary to measure the selective pressure of distinct aneuploidies—a project that is now within technological reach. Selective pressures could also be influenced by cell type (Auslander et al., 2019; Dürrbaum et al., 2014; Sack et al., 2018; Starostik et al., 2020), tumor cell genetics (Foijer et al., 2014; Grim et al., 2012; López-García et al., 2017; Simões-Sousa et al., 2018; Soto et al., 2017), and the microenvironment (Hoevenaar et al., 2020).

In summation, we developed a theoretical and experimental framework for quantitative measure of chromosomal instability in human cancer. This framework accounts for selective pressure within tumors and employs Approximate Bayesian Computation, a commonly used analysis in evolutionary biology. Additionally, we determined that low-coverage single-cell DNA sequencing of at least 200 cells from a human tumor sample is sufficient to get an accurate ( > 90% accuracy) and reproducible measure of CIN. This work sets the stage for standardized quantitative measures of CIN that promise to clarify the underlying causes, consequences, and clinical utility of this nearly universal form of genomic instability.

Materials and methods

Agent-based modeling

Agent-based models were implemented using the agent-based platform, NetLogo 6.0.4 (Wilensky, 1999).

Underlying assumptions for models of CIN and karyotype selection

Request a detailed protocol

Chromosome mis-segregation rate is defined as the number of chromosome missegregation events that occur per cellular division.

Cell division always results in 2 daughter cells.

Pmisseg,c is assigned uniformly for each cell in a population and for each chromosome.

Cells die when the copy number of any chromosome is equal to 0 or exceeding 6 unless otherwise noted.

Steps are based on the rate of division of euploid cells. We assume a probability of division (Pdivision) of 0.5, or half of the population divides every step, for euploid populations. This probabilistic division is to mimic the asynchrony of cellular proliferation and to allow for positive selection, where some cells may divide more rapidly than their euploid ancestors.

No chromosome is more likely to mis-segregate than any other.

Chromosome-arm scores

Gene abundance scores

Request a detailed protocol

The R package biomaRt v.2.46.3 was used to pull the chromosome arm location for each gene in Ensembl’s ‘Human genes’ dataset (GRCh38.p13). The number of genes on each chromosome arm were enumerated and Abundance scores were generated by normalizing the number of genes on each chromosome arm by the sum of all enumerated genes across chromosomes. Chromosome arms with no recorded genes were given a score of 0.

Driver density scores

Request a detailed protocol

Arm-level ‘TSG-OG-Ess’ scores derived in Davoli et al., 2013 were adapted for our purposes. These values were derived from a pan-cancer analysis (TCGA) of the frequency of mutation of these genes and their location in the genome. These scores correlate with the frequency with which chromosomes are found to be amplified in the genome. We adapted these scores by normalizing the published ‘TSG-OG-Ess’ score for each chromosome arm by the sum of all Charm scores. Chromosome arms with no published Charm score were given a score of 0. We refer to these as TOE scores for our purposes.

Hybrid scores

Request a detailed protocol

Chromosome arm scores for the Hybrid selection model are the average of the chromosome arm’s Gene Abundance and Driver Density scores.

Implementing karyotype selection

In each model, numerical scores are assigned to each chromosome, the sum of which represents the fitness of the karyotype (Figure 1B). At each simulation time step, fitness is re-calculated for each cell based on its updated karyotype. These fitness values determine if they undergo mitosis in the next round. However, the modality of selection changes how those karyotypes are assessed. Here, we implement four separate karyotype selection models (1) gene abundance, (2) driver density, (3) a hybrid gene abundance and driver density, and (4) neutral selection. The scores that are generated in each produce a fitness value (F) that can then be subjected to pressure (S) as described above.

Selection on gene abundance

Request a detailed protocol

The Gene Abundance selection model relies on the concept of gene dosage stoichiometry where the aneuploid karyotypes are selected against and that the extent of negative selection scales with the severity of aneuploidy and the identity and gene abundance on the aneuploid chromosomes (Sheltzer and Amon, 2011). Chromosome arm fitness contribution scores (fc) are taken as the chromosome arm scores derived above (section 2.1) and the sum of these scores is 1. These base values are then modified under the gene abundance model to generate a contextual fitness score (CFSGA,c) at each time step such that…

CFSGA,c=fc-fc×|nc-xp|xp
F=c=146CFSGA,c

… where X¯p is the average ploidy of the population and nc is the chromosome copy number. In this model, the fitness contribution of a chromosome declines as its distance from the average ploidy increases and that the magnitude of this effect is dependent on the size of the chromosome.

Selection on driver density

Request a detailed protocol

The Driver Density modality relies on assigned fitness values to chromosomes based on their relative density of tumor suppressor genes, essential genes, and oncogenes. Chromosome arm fitness contribution scores (fc) are taken as the chromosome arm scores derived above (section Driver density scores) and are employed such that…

CFSTOE,c=nc×TOEcxp
F=c=146CFSTOE,c

This selection model benefits cells that have maximized the density of oncogenes and essential genes to tumor suppressors through chromosome mis-segregation.

Hybrid selection

Request a detailed protocol

The hybrid model relies on selection on both gene abundance and driver densities. CFSTOE,c and CFSGA,c are both calculated and averaged such that…

F = c=146CFSGA,c+CFSTOE,c2

Neutral selection

Request a detailed protocol

When populations are grown under neutral selection, the fitness of each cell is constitutively set to 1 regardless of the cells’ individual karyotypes.

F = 1

Scaling selection pressure

Request a detailed protocol

Within each model of karyotype selection, the magnitude of selective pressure upon any karyotype, with fitness F, can be scaled by applying the scalar exponent S to produce a modified fitness score FM. Thus…

FM = FS

For example, in the Gene Abundance model of karyotype selection, an otherwise diploid cell with three copies of chromosome 1 in a diploid population will have a F value of 0.954. Under selection-null conditions (S = 0)…

FM=FS=0.9540=1

… the fitness of the aneuploid cell is equivalent to that of a euploid cell. Under conditions of high selection (S = 50)…

FM=FS=0.95450=0.097

…fitness of the aneuploid cell is ~10% that of the euploid cell and thus divides ~10% as frequently.

Modeling growing and constant population dynamics

To accommodate different population size dynamics, we implemented our model using either growing, pseudo-Moran limited population dynamics and constant-size populations with approximated Wright-Fisher population dynamics.

Simulating CIN in exponentially growing populations with pseudo-Moran limits

Request a detailed protocol

Populations begin with 100 founder cells with a euploid karyotype of integer value X¯p and the simulation is initiated.

CFS values are calculated for each chromosome in a cell according to the chosen karyotype selection model.

Cellular fitness is calculated based on CFS values.

Selective pressure (S) is applied to fitness (F) values to modify cellular fitness (FM).

Cells are checked to see if any death conditions are met and if the population limit is met. Cells die if any chromosome arm copy (nc) is less than 1 or greater than 6 (unless otherwise indicated). We implemented population size limits in a pseudo-Moran fashion to reduce computational constraints. If the population size is 3000 cells or greater, a random half of the population is deleted.

Cells probabilistically divide if their fitness is greater than a random float (R) between 0 and 2. Thus...

RU[0,1]

If a cell does not divide, it restarts the cycle from CFS values are calculated for each chromosome in a cell according to the chosen karyotype selection model. If a cell divides, mis-segregations may occur.

Each copy (nc) of each chromosome (c) has an opportunity to mis-segregate probabilistically. For each chromosome copy, a mis-segregation occurs if a random float (R) between 0–1 falls below Pmisseg. Thus...

RU[0,1]
MissegregatechromosomecifPmisseg,c>R

If a chromosome copy is not mis-segregated, the next chromosome copy is tested. If a chromosome copy is mis-segregated, chromosome arms may be segregated separately (i.e. a reciprocal, arm-level CNA) if a random float (R) between 0 and 1 falls below Pbreak. Thus...

RU[0,1]
BreakchromosomecifPmisseg,c>R

The karyotype of the cell is modified according to the results of the mis-segregation sequence above. When the mis-segregation sequence is complete, a clone of the initial cell with any reciprocal copy number alterations to its karyotype is created.

The simulation ends if it reaches 100 steps and data are exported. Otherwise, the simulation continues from CFS values are calculated for each chromosome in a cell according to the chosen karyotype selection model.

Simulating CIN in constant-size populations with approximated Wright-Fisher dynamics

Request a detailed protocol

We approximated constant-size Wright-Fisher dynamics in our model by re-initiating the population at each time step and randomly drawing from the previous generation’s distribution of chromosome copy numbers for each chromosome in each cell of the new population. Because the exponential pseudo-Moran model relies on proliferation rates across over-lapping generations to enact karyotype selection, such a method would not be useful here. To accommodate karyotype selection in this model, we employed an additional baseline death rate of about 20% (Sottoriva et al., 2015) that increases for cells with lower fitness and decreases for cells with higher fitness (see section 4.2.9). In this way, the karyotypes of the cells that die are removed from the pool of karyotypes that are drawn upon in the subsequent generation. CIN is simulated in this model as follows:

Populations begin with 4,500 founder cells and the simulation is (re-)initiated. The population begins with a euploid karyotype of integer value X¯p if the population is being created for the first time.

Cells divide every step, regardless of fitness.

Chromosomes are mis-segregated in the same fashion as the exponential pseudo-Moran model above (sections 4.1.8–4.1.10).

The simulation ends if it reaches 100 steps and data are exported. Otherwise, the simulation continues from 4.2.1.

CFS values are calculated for each chromosome in a cell according to the chosen karyotype selection model.

Cellular fitness is calculated based on CFS values.

Selective pressure (S) is applied to fitness (F) values to modify cellular fitness (FM).

Cells are checked to see if any death conditions are met and if the population limit is met. Cells die if any chromosome arm copy (nc) is less than 1 or greater than 6 (unless otherwise indicated).

Additionally, the cells’ fitness values and a random float (R) between 0 and 5 are used to determine if they die. In this way, a cell with a fitness of 1 has a 20% baseline death rate. Thus, cells die if…

1FS+0.001>RU[0,5]

After determining cell death, the copy number distributions of each cells’ chromosome arm (c) are individually stored.

The cycle repeats from 4.2.1. However, the re-initated population will have its chromosome arm copy numbers drawn from the previous generation’s stored chromosome arm copy number distributions.

Analysis of population diversity and topology in biological and simulated data

Request a detailed protocol

Phylogenetic trees were reconstructed from chromosome copy number profiles from live and simulated cells by calculating pairwise Euclidean distance matrices and performing complete-linkage clustering in R (R Development Core Team, 2021). Phylogenetic tree topology measurements were performed in R using the package phyloTop v2.1.1 (Kendall et al., 2018). Sackin and Colless indices of tree imbalance were calculated, normalizing to the number of tree tips. Cherry and pitchfork number were also normalized to the size of the tree. MKV is taken as the variance of individual chromosomes taken across the population, averaged across all chromosomes, then normalized to the average ploidy of the population. Average aneuploidy is calculated as the variance within a single cell’s karyotype averaged across the population.

Approximate bayesian computation

Request a detailed protocol

Approximate Bayesian computation was used for parameter inference of experimental data from simulated data. For this we employed the the “abc” function in the R package abc v2.1 (Csilléry et al., 2010). In short, a set of simulation parameters, θi, is sampled from the prior distribution. This set of parameters corresponds to a set of simulated summary statistics, S(yi), in this case phylogenetic tree shapes, which can be compared to the set of experimental summary statistics, S(yo). The Euclidean distance between the experimental and simulated summary statistics can then be calculated (dS(yi),S(yo)). A threshold, T, is then selected—0.05 in our case—which rejects the lower 1 T sets of simulation parameters that correspond. The remaining parameters represent those that gave summary statistics with the highest similarity to the experimental summary statistics. These represent the posterior distribution of accepted parameters.

Bayesian model selection was performed using the “postpr” function in the same R package using tolerance threshold of 0.05 and rejection sampling method. This was used to calculate the posterior probability of each selection model within each growth model and the Bayes factor for each selection model with neutral selection as the null hypothesis. Bayes factors > 5 were considered substantial evidence of the alternative hypothesis.

Sliding window analysis to tune time-steps for approximate Bayesian computation

Request a detailed protocol

We chose which simulation time steps to use for approximate Bayesian computation on organoid and biopsy data by repeating the inference using a sliding window of prior datasets with a width of 11 time steps (i.e. parameters from steps ∈ [0–10], [10-20], …, [91-100]) to see if the posterior distributions would stabilize over time. We then chose simulations from 40 to 80 time steps as our prior dataset as this range provided both a stable inference and is centered around 60 time steps (analogous to 30 generations, estimated to generate a 1 cm palpable mass of ~1 billion cells).

Cell cultivation procedures

Request a detailed protocol

Cal51 cells expressing stably integrated RFP-tagged histone H2B and GFP-tagged a-tubulin were generated as previously described (Zasadil et al., 2014). Cells were maintained at 37 ºC and 5% CO2 in a humidified, water-jacketed incubator and propagated in Dulbecco’s Modified Eagle’s Medium (DMEM) – High Glucose formulation (Cat #: 11965118) supplemented with 10% fetal bovine serum and 100 units/mL penicillin-streptomycin. Paclitaxel (Tocris Bioscience, Cat #: 1097/10) used for cell culture experiments was dissolved in DMSO. The Cal51 cells were obtained from the DSMZ-German Collection of Microorganisms and Cell Cultures and were free from mycoplasma contamination prior to study. Karyotype analysis confirms the near-diploid characteristic of the cell line and the presence of both fluorescent markers suggests they are free of other contaminating cell lines.

Time-lapse fluorescence microscopy

Request a detailed protocol

Cal51 cells were transduced with lentivirus expressing mNeonGreen-tubulin-P2A-H2B-FusionRed. A monoclonal line was treated with 20 nM paclitaxel for 24, 48, or 72 hr before timelapse analysis at 37 oC and 10% CO2. Five 2 µm z-plane images were acquired using a Nikon Ti-E inverted microscope with a cMos camera at 3-min intervals using a 40 X/0.75 NA objective lens and Nikon Elements software.

Flow cytometric analysis and cell sorting

Cells were harvested with trypsin, passed through a 35 μm mesh filter, and rinsed with PBS prior to fixation in ice cold 80% methanol. Fixed cells were stored at –80 ºC until analysis and sorting at which point fixed cells were resuspended in PBS containing 10 μg/ml DAPI for cell cycle analysis.

Flow cytometric analysis

Request a detailed protocol

Initial DNA content and cell cycle analyses were performed on a 5 laser BD LSR II. Doublets were excluded from analysis via standard FSC/SSC gating procedures. DNA content was analyzed via DAPI excitation at 355 nm and 450/50 emission using a 410 nm long pass dichroic filter.

Fluorescence activated cell sorting

Request a detailed protocol

Cell sorting was performed using the same analysis procedures described above on a BD FACS AriaII cell sorter. In general, single cells were sorted through a 130 μm low-pressure deposition nozzle into each well of a 96-well PCR plate containing 10 μl Lysis and Fragmentation Buffer cooled to 4 ºC on a Eppendorf PCR plate cooler. Immediately after sorting PCR plates were centrifuged at 300 x g for 60 s. For comparison of single-cell sequencing to bulk sequencing, 1000 cells were sorted into each ‘bulk’ well. The index of sorted cells was retained allowing for the post hoc estimation of DNA content for each cell.

Low-coverage single-cell whole genome sequencing

Request a detailed protocol

Initial library preparation for low-coverage scDNAseq was performed as previously described (Leung et al., 2016) and adapted for low coverage whole genome sequencing instead of high coverage targeted sequencing. Initial genome amplification was performed using the GenomePlex Single Cell Whole Genome Amplification Kit and protocol (Sigma Aldrich, Cat #: WGA4). Cells were sorted into 10 μl pre-prepared Lysis and Fragmentation buffer containing Proteinase K. DNA was fragmented to an average of 1 kb in length prior to amplification. Single cell libraries were purified on a 96-well column plate (Promega, Cat #: A2271). Library fragment distribution was assessed via agarose gel electrophoresis and concentrations were measured on a Nanodrop 2000. Sequencing libraries were prepared using the QuantaBio sparQ DNA Frag and Library Prep Kit. Amplified single-cell DNA was enzymatically fragmented to ~250 bp, 5’-phosphorylated, and 3’-dA-tailed. Custom Illumina adapters with 96 unique 8 bp P7 index barcodes were ligated to individual libraries to enable multiplexed sequencing (Leung et al., 2016). Barcoded libraries were amplified following size selection via AxygenAxyPrep Mag beads (Cat #: 14-223-152). Amplified library DNA concentration was quantified using the Quant-iT Broad-Range dsDNA Assay Kit (Thermo, Cat #: Q33130). Single-cell libraries were pooled to 15 nM and final concentration was measured via qPCR. Single-end 100 bp sequencing was performed on an Illumina HiSeq2500.

Single-cell copy number sequencing data processing

Request a detailed protocol

Single-cell DNA sequence reads were demultiplexed using unique barcode index sequences and trimmed to remove adapter sequences. Reads were aligned to GRCh38 using Bowtie2. Aligned BAM files were then processed using Ginkgo to make binned copy number calls. Reads are aligned within 500 kb bins and estimated DNA content for each cell, obtained by flow cytometric analysis, was used to calculate bin copy numbers based on the relative ratio of reads per bin (Garvin et al., 2015). We modified and ran Ginkgo locally to allow for the analysis of highly variable karyotypes with low ploidy values (see Code and Data Availability). Whole-chromosome copy number calls were calculated as the modal binned copy number across an individual chromosome. Cells with fewer than 100,000 reads were filtered out to ensure accurate copy number calls (Baslan et al., 2015). Cells whose predicted ploidy deviated more than 32% from the observed ploidy by FACS were also filtered out. The final coverage for the filtered dataset was 0.03 (5). Single cell data extracted from Navin et al., 2011 were separated into their individual clones and depleted of euploid cells. Single cell data from Bolhaqueiro et al., 2019 were filtered to include only the aneuploid data that fell within the ploidies observed in the study (see Code and Data Availability).

Review and approximation of mis-segregation rates from published Studies

Request a detailed protocol

We reviewed the literature to extract per chromosome rates of mis-segregation for cell lines and clinical samples. Some studies publish these rates. For those that did not, we estimated these rates by approximating the plotted incidence of segregation errors thusly:

Approximatemissegregrationrateperchromosome=Observed%frequencyoferrorsperdivision/100Total#modalchromosomesinsample

Modal chromosome numbers were either taken from ATCC where available or were assumed to equal 46. Observed % frequencies were approximated from published plots. Approximated rates assume that 1 chromosome is mis-segregated at a time.

Data availability

Single-cell DNA sequencing data from this study has been deposited in NCBI SRA (PRJNA725515). All data and scripts used for modeling and analysis have been deposited in OSF at https://osf.io/snrg3/.

The following data sets were generated
    1. Lynch AR
    2. Arp NL
    3. Zhou AS
    4. Weaver BA
    5. Burkard ME
    (2021) Open Science Framework
    Quantifying chromosomal instability from intratumoral karyotype diversity using agent- based modeling and Bayesian inference.
    https://doi.org/10.17605/OSF.IO/SNRG3
    1. Lynch AR
    2. Arp NL
    3. Zhou AS
    4. Weaver BA
    5. Burkard ME
    (2021) NCBI BioProject
    ID PRJNA725515. Quantifying chromosomal instability from intratumoral karyotype diversity Quantifying chromosomal instability from intratumoral karyotype diversity.

References

  1. Website
    1. Wilensky U
    (1999) NetLogo
    Accessed March 11, 2020.

Decision letter

  1. Adèle L Marston
    Reviewing Editor; University of Edinburgh, United Kingdom
  2. Anna Akhmanova
    Senior Editor; Utrecht University, Netherlands
  3. Trevor A Graham
    Reviewer; Queen Mary, University of London, United Kingdom

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Quantifying chromosomal instability from intratumoral karyotype diversity using agent- based modeling and Bayesian inference" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Anna Akhmanova as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Trevor A Graham (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. Modelling.

Three different models to compute cell fitness based on karyotypic alteration are explored. The construction of all these models feels a little arbitrary, and the assumptions and evolutionary dynamics in each scenario should be more comprehensively explored. Specifically:

a) In the TOE model, fitness is inversely proportional to average ploidy, so it seems higher ploidies are always selected against. Is this a reasonable assumption? Why is it necessary to divide by the average ploidy?

b) In all models, is the simulated population always "out of equilibrium". If the simulations ran for longer would an "optimal karyotype" be established. Relatedly, the dynamics appear to be strongly influenced by the copy number >6 being lethal – chromosomes (in the TOE model) which are beneficial to be gained might tend to increase copy number to 5 whereas deleterious gains reduce copy number to 1 and the population rests on that "precipice". How reasonable are these "boundary conditions" and do the dynamics change significantly if they are relaxed?

c) Gains and losses appear to be treated equivalently – again is this reasonable? Especially in the TOE model where TSG gains and OG losses (and vice versa) have differing consequences. (see also point 8 below).

The modelling assumes exponential growth to a relatively small number of cells (4500) and then randomly kills half the cells to reinitiate exponential growth from 2500 cells. This regime will influence the evolutionary dynamics of the system: the random killing will cause emerging clones to often go extinct and could exacerbate the influence of drift in the system. It also effects the influence of selection, (see for example: https://www.nature.com/articles/ng.3214). Alternative growth dynamics could be implemented – such as a Wright Fisher type model of either constant or growing populations (for the construction of a growing WF see: https://pubmed.ncbi.nlm.nih.gov/25665006/) and the influence of the growth dynamics on karyotypic heterogeneity robustly assessed.

The modelling assumes only whole chromosome missegregations (and the bioinformatics data analysis averages over sub-chromsomal sized events). The authors could consider extending their analysis to handle part-chromosome events to better represent the biological data.

2. Bayesian inference.

It is a concern that unusual prior distributions are used in the ABC inference and this effects the reliability of the inference. Figure 5F and 6C show smoothed density plots for the prior distributions – which confusingly show density for S<0 – and the true priors might instead be a series of point masses at a handful of S values. This should be clarified.

It is likely that the posterior the alteration rate and S are interrelated (high S inferred when the alteration rate is high and vice versa) – so joint posteriors should be shown. Because of the interrelationship between the parameters, it is a concern that the current parameter estimates are inaccurate – currently the prior for S has zero mass for many S values, and the inferred value of the alteration rate will depend on which values of S are explored form the prior distribution. The inference should be repeated using continuous distribution over S, a uniform distribution is suggested.

When real data is analysed, only the hybrid model is compared to data, but their Figure 3 shows the diversity depends on the underlying model of selection. The authors should implement a model selection routine (one is available with the ABC-sysbio package: https://pubmed.ncbi.nlm.nih.gov/20591907/) to test which selection paradigm (if any) best represents the data. They could also consider comparing what is arguably the "null hypothesis" of neutral evolution (S=0) to a case with selection as part of this model selection, to quantitatively determine the evidence of selection in the data.

A single threshold for acceptance in the ABC algorithm (epsilon=0.05). The authors should show this is sufficiently small: a straightforward way is to plot the mean and variance of the posteriors as a function of epsilon.

Above, the possibility of temporal effects in the model are mentioned – how do temporal ("out of equilibrium" evolutionary dynamics) affect the inference?

When public data are analysed, the growth dynamics of the breast tissue, or colorectal organoid is very unlikely to match the exponential growth assumption that is assumed by the authors. These growth dynamics likely strongly determine the pattern of observed heterogeneity (exponential growth leads to large founder effects in the extant population) and so the influence of alterative growth models needs to be explored. This can also be done within a model selection framework.

The initial state of the model used to fit to public dataset is not specified and needs to be. The authors might also consider the influence of spatial sampling confounders (the breast dataset in particular in not a well-mixed population).

3. Single cell data analysis. The data are presented as if only whole chromosomal alterations occur (e.g. figure 5B). Is this actually the case? (certainly the assumption is false in breast and colon organoid datasets) Could relative read counts in say 1MB bins for each cell be provided in the appendix to reassure the reader that the types of genetic alterations occurring in their single cell data mirror the assumptions of the modelling.

As noted above, the authors should consider modelling part-chromosome alterations. They should be able to determine in their simulations how their assumption of only whole chromosome missegregation events, despite their being part-chromosomal copy number alterations in the data, affects the accuracy of the inferred chromosomal missegregation rate.

4. The presented framework is lacking expanded characterization and validation of selection models that are biologically relevant. The current framework simply applies a scalar exponent to already published fitness models for selection. It is unclear what this exponent mirrors biologically, beyond amplifying the selection pressures already explored in existing gene abundance and driver density models.

5. Related to point (4)., how is the CIN ON-OFF model in which CIN is turned off after so many cell divisions relevant biologically? Typically, CIN is a considered a trait that evolves later in cancer progression, that once tolerated, is ongoing and facilitates development of metastasis and drug resistance. A more relevant model to explore would be that of the effect of a whole genome duplication (WGD) event on population evolution, which is thought to facilitate tolerance of ensuing missegregation events (because reduce risk of nullisomy).

5. The authors utilize two models of karyotype fitness – a gene abundance model and driver density model – to evaluate impact of specific karyotypes on cellular fitness. They also include a hybrid model whose fitness effects are simply the average of these two models, which adds little value as only a weighted average. In silico results shows inferred missegregation rates are extremely disparate across the two primary models. And while a description of these differences is provided, the presented analyses do not make clear the most important question – which of these models is more clinically relevant? Toward this, in Figure 2F, the authors claim the three models approach a triploid state – which is unsupported by the in silico results. Clearly the driver model approaches a triploid state, as previously reported. But the abundance model does not and hybrid only slightly so, given that it is simply a weighted average of these two approaches. Because the authors have developed a Bayesian strategy for inferring which model parameters best fit observed data, it would be very useful to see which model best recapitulates karyotypes observed in cancer cell lines or patient materials.

6. Topological features of phylogenetic trees, while discriminatory, are largely dependent on accurate phylogenetic tree reconstruction. The latter requires more careful consideration of cell linkages beyond computing pairwise Euclidean distances and performing complete-linkage clustering. For example, a WGD event, would appear very far from its nearest cell ancestor in Euclidean space.

7. Experimental validation of the added selection exponential factor is imperative. Works have already shown models of karyotypic evolution without additional selection exponential coefficient can accurately recover rates of missegregation observed in human cell lines and cancers by fluorescent microscopy. Incorporation of this additional weight on selection pressure has not been demonstrated or validated experimentally. This would require experimental sampling of karyotypes longitudinally and is a critical piece of this manuscript's novelty.

8. It seems like this model treats chromosome gains and losses equivalently. Is this appropriate? Chromosome loss events are much more toxic than chromosome gain events – as evidenced by the fact that haploinsufficiency is widespread, and all autosomal monosomies are embryonically-lethal while many trisomies are compatible with birth and development. Can the authors consider a model in which losses exert a more significant fitness penalty that chromosome gains? (see also point 1c above).

9. Chromosomes do not missegregate at the same rate (PMID: 29898405). This point should be discussed, and, if feasible, incorporated into the authors' models.

10. Can the authors clarify their use of live cell imaging (e.g., in Figure 6G)? Certain apparent errors that are visible by live-cell imaging (like a lagging chromosome) can be resolved correctly and result in proper segregation. Is it appropriate to directly interfer missegregation rates as is done in this paper?

11. The authors should discuss in greater detail earlier mathematical models of CIN, including PMID: 26212324, 30204765, and 12446840. How does their approach improve on this prior work?

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

Thank you for resubmitting your work entitled "Quantifying chromosomal instability from intratumoral karyotype diversity using agent- based modeling and Bayesian inference" for further consideration by eLife. Your revised article has been evaluated by Anna Akhmanova (Senior Editor), Reviewing Editors and two of the original reviewers.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

1. The reviewers are concerned about the choice of statistics used for the ABC analysis. While these summary statistics indeed must contain some statistical signal on the rate of chromosomal instability, the current method of "phylogenetic reconstruction" which based on euclidian distances of chromosomal copy number variation is very indirect. The reviewers are not convinced that this is an optimal way to extract relevant information from the distribution of chromosomal copy number across cells. At least they are not clearly rooted in models of clonal evolution.

Why this particular set of summary statistics was chosen, how it relates to the quantity of interest (rate of chromosomal instability), how it compares to other possible ways to capture relevant information must be properly justified. The choice of these statistics should be explored in much more details in the manuscript with clear rationale given and caveats with this choice clearly stated.

2. Related to point 1 – The phylogenetic reconstruction is an unusual way to summarise the statistical data and has not been rigorously assessed. The authors should elaborate on the consequences of the possibly low accuracy of the tree reconstruction itself.

3. The manuscript is difficult to follow. For example, there is a lengthy part on validating simulations that doesn't have a very specific message and could be reduced quite a lot. Please also ensure that the manuscript narrative is accessible to the general readership of eLife, including defining specialist concepts, referring to figures (and figure supplements) in the order they appear in the text. Please also consider adding additional labels to the figures themselves to aid the reader and avoid over-simplification e.g. "Tet" "Dip" could be written in full, Figure 2D is missing a key etc.

4. More details are required on the ABC analysis, please address the additional comments of reviewer #3 below.

Reviewer #3 (Recommendations for the authors):

I focused my attention at looking at the revised ABC analysis.

I found the new analysis a little hard to follow.

It is not clear to me precisely what the "step windows" and sliding/expanding timesteps are in Fig6-S1&2. I couldn't see this explicitly explained in the methods. (I can see that the authors are trying to optimise the simulated number of cell divisions, I presume some kind of adaptive search was done)

I couldn't understand the rationale for only considering the mutation rate in these optimisations. I think that because stronger selection suppresses diversity, the strength of selection should inversely correlate with the time needed to evolve the diversity observed in the organoids.

I'm concerned about the posterior distributions shown in Figure 6-S3. In most cases the mass in the posteriors is at the extremes of the prior, which suggests that the true parameter values lie outside of the range of the prior (i.e. that the priors were too narrow). This problem appears always the case for the selection coefficient and appears to be a problem for the mutation rate in some cases too (at least it seems that the mutation rate cannot be distinguished from the boundary of 0). Possibly the same issue of a to narrow a prior applies to Figure 5F (taxol treated cell lines). The concern is that the data could be much better explained by a set of parameters located outside of the parameter space considered, which would mean that the conclusions could be incorrect.

Calculating posterior predictive distributions (analagous to Figure 6-S4&5) would help to convince the goodness of fit.

Line 479 says that a large-ish value of the acceptance threshold (epsilon=0.05) was used to retain prior data. This seems an incorrect statement: epsilon determines how many simulations are rejected but all prior data is always used. Smaller epsilons should give greater accuracy but at the cost of more computation.

https://doi.org/10.7554/eLife.69799.sa1

Author response

Essential revisions:

1. Modelling.

Three different models to compute cell fitness based on karyotypic alteration are explored. The construction of all these models feels a little arbitrary, and the assumptions and evolutionary dynamics in each scenario should be more comprehensively explored. Specifically:

a) In the TOE model, fitness is inversely proportional to average ploidy, so it seems higher ploidies are always selected against. Is this a reasonable assumption? Why is it necessary to divide by the average ploidy?

This normalization was done to account for gene balance. This fundamental premise of this model is neutral selection if all chromosomes are balanced without excess oncogenes or reduced tumor suppressors. This balance would occur equally with diploid, triploid, and tetraploid cells. By including this normalization, only imbalances which create excess oncogenes or reduce suppressors relative to the entire genome impart selection. We believe this model is most consistent with the known and specific roles of oncogenes and tumor suppressors in promoting or restraining cell proliferation.

b) In all models, is the simulated population always "out of equilibrium". If the simulations ran for longer would an "optimal karyotype" be established. Relatedly, the dynamics appear to be strongly influenced by the copy number >6 being lethal – chromosomes (in the TOE model) which are beneficial to be gained might tend to increase copy number to 5 whereas deleterious gains reduce copy number to 1 and the population rests on that "precipice". How reasonable are these "boundary conditions" and do the dynamics change significantly if they are relaxed?

We believe these boundaries are reasonable as large scale copy number alterations higher than this are rare (see PMID: 31036964 and 32054838). However, to address this, we implemented a variant of our model that considers alternative thresholds. Additionally, we agree that the CIN ON-OFF model had limited biologic relevance and removed this. To improve on this, we have changed our approach to use constant CIN for a much longer period of time (3000 time steps). We agree that WGD is a relevant phenomenon. However, others have already explicitly modeled this (see PMID: 26212324 and 32139907), so we avoid doing the same. Instead, we show that tetraploid founding cells tolerate high mis-segregation rates better than diploid founding cells.

c) Gains and losses appear to be treated equivalently – again is this reasonable? Especially in the TOE model where TSG gains and OG losses (and vice versa) have differing consequences. (see also point 8 below).

While we agree that monosomies are more detrimental than trisomies in non-cancerous tissue, this is not necessarily the case in tumors in which monosomy is often observed (see PMID: 32054838). Nevertheless, to address this critique we have now added a model variant with an additional condition in which cells experience extreme fitness penalties (90% reduction) if any chromosome is haploid. We apply this condition to all selection models and find this attenuates a ploidy increase over time in diploid cells in most selection models (see Figure 3 ‘haploid penalty’).

The modelling assumes exponential growth to a relatively small number of cells (4500) and then randomly kills half the cells to reinitiate exponential growth from 2500 cells. This regime will influence the evolutionary dynamics of the system: the random killing will cause emerging clones to often go extinct and could exacerbate the influence of drift in the system. It also effects the influence of selection, (see for example: https://www.nature.com/articles/ng.3214). Alternative growth dynamics could be implemented – such as a Wright Fisher type model of either constant or growing populations (for the construction of a growing WF see: https://pubmed.ncbi.nlm.nih.gov/25665006/) and the influence of the growth dynamics on karyotypic heterogeneity robustly assessed.

To address these concerns, we improved and more clearly detailed the prior distributions for each inference within the figure legends, we tested for karyotype convergence in each model (see Figure 3), and we demonstrate that inference under the Abundance model is robust to changes in the number of time steps included in the prior data (see Figure 6 — figure supplement 1).

The modelling assumes only whole chromosome missegregations (and the bioinformatics data analysis averages over sub-chromsomal sized events). The authors could consider extending their analysis to handle part-chromosome events to better represent the biological data.

Thank you for this critique. We have now extended the model to handle arm-level segmental aneuploidy.

2. Bayesian inference.

It is a concern that unusual prior distributions are used in the ABC inference and this effects the reliability of the inference. Figure 5F and 6C show smoothed density plots for the prior distributions – which confusingly show density for S<0 – and the true priors might instead be a series of point masses at a handful of S values. This should be clarified.

We now provide data simulated using uniform distributions of mis-segregation rate and selection.

It is likely that the posterior the alteration rate and S are interrelated (high S inferred when the alteration rate is high and vice versa) – so joint posteriors should be shown. Because of the interrelationship between the parameters, it is a concern that the current parameter estimates are inaccurate – currently the prior for S has zero mass for many S values, and the inferred value of the alteration rate will depend on which values of S are explored form the prior distribution. The inference should be repeated using continuous distribution over S, a uniform distribution is suggested.

We now provide the joint posterior distributions inferred from uniform prior distributions in Figure 6C.

When real data is analysed, only the hybrid model is compared to data, but their Figure 3 shows the diversity depends on the underlying model of selection. The authors should implement a model selection routine (one is available with the ABC-sysbio package: https://pubmed.ncbi.nlm.nih.gov/20591907/) to test which selection paradigm (if any) best represents the data. They could also consider comparing what is arguably the "null hypothesis" of neutral evolution (S=0) to a case with selection as part of this model selection, to quantitatively determine the evidence of selection in the data.

We now include neutral evolution. To address these concerns, we improved and more clearly detailed the prior distributions for each inference within the figure legends, we tested for karyotype convergence in each model (see Figure 3), and we demonstrate that inference under the Abundance model is robust to changes in the number of time steps included in the prior data (see Figure 6 — figure supplement 1).

A single threshold for acceptance in the ABC algorithm (epsilon=0.05). The authors should show this is sufficiently small: a straightforward way is to plot the mean and variance of the posteriors as a function of epsilon.

This is now addressed in the new Figure 6 — figure supplement 1.

Above, the possibility of temporal effects in the model are mentioned – how do temporal ("out of equilibrium" evolutionary dynamics) affect the inference?

To address these concerns, we improved and more clearly detailed the prior distributions for each inference within the figure legends, we tested for karyotype convergence in each model (see Figure 3), and we demonstrate that inference under the Abundance model is robust to changes in the number of time steps included in the prior data (see Figure 6 — figure supplement 1).

When public data are analysed, the growth dynamics of the breast tissue, or colorectal organoid is very unlikely to match the exponential growth assumption that is assumed by the authors. These growth dynamics likely strongly determine the pattern of observed heterogeneity (exponential growth leads to large founder effects in the extant population) and so the influence of alterative growth models needs to be explored. This can also be done within a model selection framework.

We have now concurrently modeled chromosomal instability with a constant population size by approximating constant-population Wright Fisher dynamics (see Materials and methods). We find these models produce similar results at the karyotype level, addressing concerns about the effects of growth patterns on karyotype evolution in this model.

The initial state of the model used to fit to public dataset is not specified and needs to be. The authors might also consider the influence of spatial sampling confounders (the breast dataset in particular in not a well-mixed population).

We have updated the figure legends with more detailed prior and parameter settings.

3. Single cell data analysis. The data are presented as if only whole chromosomal alterations occur (e.g. figure 5B). Is this actually the case? (certainly the assumption is false in breast and colon organoid datasets) Could relative read counts in say 1MB bins for each cell be provided in the appendix to reassure the reader that the types of genetic alterations occurring in their single cell data mirror the assumptions of the modelling.

Paclitaxel causes whole-chromosome aneuploidy through chromosome mis-segregation (Scribano et al. Sci Trans Med 2021). As requested, the binned read counts are now illustrated in the new Figure 5 — figure supplement 2.

As noted above, the authors should consider modelling part-chromosome alterations. They should be able to determine in their simulations how their assumption of only whole chromosome missegregation events, despite their being part-chromosomal copy number alterations in the data, affects the accuracy of the inferred chromosomal missegregation rate.

As requested, we have now performed simulations that allow for segmental aneuploidy.

4. The presented framework is lacking expanded characterization and validation of selection models that are biologically relevant. The current framework simply applies a scalar exponent to already published fitness models for selection. It is unclear what this exponent mirrors biologically, beyond amplifying the selection pressures already explored in existing gene abundance and driver density models.

Biologically, this mirrors the extent to which aneuploid karyotypes are selected for or against, given a particular model. It is an exponent rather than a multiplier because the selection value is already transformed to a probability of division. We now provide model selection and further validation of this method.

To address this, we have greatly expanded the models and their characterization. We now explicitly include a neutral model throughout, tested various modifications of the model (Figure 3C-E), and use ABC to enable model selection (see Table 3).

We implemented cellular fitness as the sum of normalized chromosome scores such that the fitness of euploid cells is 1 and the probability of division = 0.5. In this framework, within the ‘abundance’ model, a cell with triploidy of chromosome arm 1p would have a fitness of 0.98. With no additional selection, the probability that this cell divides is 0.98 x 0.5 = 0.49.

The published fitness models for karyotype selection do not experimentally determine how fitness relates to the probability of division within a given time. For example, there is no clear reason why (or evidence indicating) an extra copy of chromosome arm 1p would reduce the probability of division from 0.5 to precisely 0.49 for a given period. The proposed model of karyotype selection that our ‘abundance’ model is based on only stipulates that aneuploidy of larger chromosomes is more detrimental than small chromosomes. Thus, these fitness values behave as arbitrary units and, therefore, we believe that adjusting and fitting an arbitrary scaling factor to the biological data is appropriate. For example, with an additional selection of S=10, the same cell with trisomy of chromosome arm 1p would divide with a probability of FS x 0.5 = 0.9810 x 0.5 = 0.41.

We could have implemented a multiplicative framework where fitness (Fmult) is defined as the total deviation from euploid fitness (1) multiplied by a scaling factor S (Fmult = S(1 - F)). For the trisomy 1p example, the same fitness value (FS=0.9810) can be achieved multiplicatively as exponentially via 1 – (9.14 x (1 - 0.98)) ~ 0.9810. Thus, the same fitness values can be achieved through arbitrary scaling. We regret that this may have been misinterpreted because it was implemented exponentially vs multiplicatively.

To further address this critique, we have now better fitted the S values with a flat prior probability across all values, shown how it relates to Pmisseg in posterior probabilities (e.gs, Figure 6C, Table 3) and performed the separate analysis requested.

The selection values of F are in arbitrary units and so we believe a selection scaling factor is important to include in the model. For example, without additional selection, a hypothetical aneuploid cell with a trisomy resulting in F = 0.95 would be 5% less likely to divide than a euploid cell with F = 1. The exponent scales the selection such that when S = 2, the fitness of the trisomic cell is F ~ 0.9, or 10% less likely to divide. This scaling is necessary to enable both positive and negative selection in a system fitness is decided as the sum of chromosome scores. To further validate the additional weight on selection pressure we did the following:

1) We constrained the prior distribution of simulated data for our model selection to S=1 giving only the base fitness values without additional scaling. We, again, performed model selection on the data from Bolhaqueiro et al., 2019 and Navin et al., 2011 and found that, with this constrained prior dataset, we inferred mis-segregation rates (see Table 4) that were far below rates seen in cancer cell lines (see Figure 6E).

2) Given the initial clarification that reviewers were looking for longitudinal analysis, we leveraged data provided by the authors of Bolhaqueiro et al., 2019 where they sequenced single cells from 3 clones from organoid line 16T at 3 weeks and 21 weeks after seeding. We inferred mis-segregation rates and selective pressures in these clones at the 3-week timepoint. We did so under the Abundance model using the same prior distribution of steps given that the diversity of populations under the Abundance model rapidly reach a steady state. When we simulated additional populations using these inferred characteristics we found that the karyotype composition of the simulated populations most closely resembled the biological population than did populations simulated with the unmodified selection values (see Figure 6 — figure supplement 4). This lends credence to the biological relevance of scaled selective pressure vs. unmodified selective pressure.

We have now concurrently modeled chromosomal instability with a constant population size by approximating constant-population Wright Fisher dynamics (see Materials and Methods). We find these models produce similar results at the karyotype level, addressing concerns about the effects of growth patterns on karyotype evolution in this model.

5. Related to point (4)., how is the CIN ON-OFF model in which CIN is turned off after so many cell divisions relevant biologically? Typically, CIN is a considered a trait that evolves later in cancer progression, that once tolerated, is ongoing and facilitates development of metastasis and drug resistance. A more relevant model to explore would be that of the effect of a whole genome duplication (WGD) event on population evolution, which is thought to facilitate tolerance of ensuing missegregation events (because reduce risk of nullisomy).

We agree that this may not be relevant biologically and have removed the CIN ON-OFF scheme and updated Figure 3 to, instead, explore the convergence of karyotypes. Additionally, we agree that the CIN ON-OFF model had limited biologic relevance and removed this. To improve on this, we have changed our approach to use constant CIN for a much longer period of time (3000 time steps). We agree that WGD is a relevant phenomenon. However, others have already explicitly modeled this (see PMID: 26212324 and 32139907), so we avoid doing the same. Instead, we show that tetraploid founding cells tolerate high mis-segregation rates better than diploid founding cells.

5. The authors utilize two models of karyotype fitness – a gene abundance model and driver density model – to evaluate impact of specific karyotypes on cellular fitness. They also include a hybrid model whose fitness effects are simply the average of these two models, which adds little value as only a weighted average. In silico results shows inferred missegregation rates are extremely disparate across the two primary models. And while a description of these differences is provided, the presented analyses do not make clear the most important question – which of these models is more clinically relevant? Toward this, in Figure 2F, the authors claim the three models approach a triploid state – which is unsupported by the in silico results. Clearly the driver model approaches a triploid state, as previously reported. But the abundance model does not and hybrid only slightly so, given that it is simply a weighted average of these two approaches. Because the authors have developed a Bayesian strategy for inferring which model parameters best fit observed data, it would be very useful to see which model best recapitulates karyotypes observed in cancer cell lines or patient materials.

We agree that the CIN ON-OFF model had limited biologic relevance and removed this. To improve on this, we have changed our approach to use constant CIN for a much longer period of time (3000 time steps). We agree that WGD is a relevant phenomenon. However, others have already explicitly modeled this (see PMID: 26212324 and 32139907), so we avoid doing the same. Instead, we show that tetraploid founding cells tolerate high mis-segregation rates better than diploid founding cells.

6. Topological features of phylogenetic trees, while discriminatory, are largely dependent on accurate phylogenetic tree reconstruction. The latter requires more careful consideration of cell linkages beyond computing pairwise Euclidean distances and performing complete-linkage clustering. For example, a WGD event, would appear very far from its nearest cell ancestor in Euclidean space.

We agree that the abundance and hybrid models are unable to approach a triploid state, in earnest, as does the driver and have made that clearer in the text and improved the figure panel in question for clarity. To address your latter point on which model best fits observed data, we have implemented a model selection scheme to do this (see Table 3). This indicates the gene abundance model as the most biologically relevant and provides evidence for stabilizing selection as the primary mode of selection occurring in the organoid and biopsy data we analyzed.

7. Experimental validation of the added selection exponential factor is imperative. Works have already shown models of karyotypic evolution without additional selection exponential coefficient can accurately recover rates of missegregation observed in human cell lines and cancers by fluorescent microscopy. Incorporation of this additional weight on selection pressure has not been demonstrated or validated experimentally. This would require experimental sampling of karyotypes longitudinally and is a critical piece of this manuscript's novelty.

The selection values of F are in arbitrary units and so we believe a selection scaling factor is important to include in the model. For example, without additional selection, a hypothetical aneuploid cell with a trisomy resulting in F = 0.95 would be 5% less likely to divide than a euploid cell with F = 1. The exponent scales the selection such that when S = 2, the fitness of the trisomic cell is F ~ 0.9, or 10% less likely to divide. This scaling is necessary to enable both positive and negative selection in a system fitness is decided as the sum of chromosome scores. To further validate the additional weight on selection pressure we did the following:

1) We constrained the prior distribution of simulated data for our model selection to S=1 giving only the base fitness values without additional scaling. We, again, performed model selection on the data from Bolhaqueiro et al., 2019 and Navin et al., 2011 and found that, with this constrained prior dataset, we inferred mis-segregation rates (see Table 4) that were far below rates seen in cancer cell lines (see Figure 6E).

2) Given the initial clarification that reviewers were looking for longitudinal analysis, we leveraged data provided by the authors of Bolhaqueiro et al., 2019 where they sequenced single cells from 3 clones from organoid line 16T at 3 weeks and 21 weeks after seeding. We inferred mis-segregation rates and selective pressures in these clones at the 3-week timepoint. We did so under the Abundance model using the same prior distribution of steps given that the diversity of populations under the Abundance model rapidly reach a steady state. When we simulated additional populations using these inferred characteristics we found that the karyotype composition of the simulated populations most closely resembled the biological population than did populations simulated with the unmodified selection values (see Figure 6 — figure supplement 4). This lends credence to the biological relevance of scaled selective pressure vs. unmodified selective pressure.

8. It seems like this model treats chromosome gains and losses equivalently. Is this appropriate? Chromosome loss events are much more toxic than chromosome gain events – as evidenced by the fact that haploinsufficiency is widespread, and all autosomal monosomies are embryonically-lethal while many trisomies are compatible with birth and development. Can the authors consider a model in which losses exert a more significant fitness penalty that chromosome gains? (see also point 1c above).

While we agree that monosomies are more detrimental than trisomies in non-cancerous tissue, this is not necessarily the case in tumors in which monosomy is often observed (see PMID: 32054838). Nevertheless, to address this critique we have now added a model variant with an additional condition in which cells experience extreme fitness penalties (90% reduction) if any chromosome is haploid. We apply this condition to all selection models and find this attenuates a ploidy increase over time in diploid cells in most selection models (see Figure 3 ‘haploid penalty’).

9. Chromosomes do not missegregate at the same rate (PMID: 29898405). This point should be discussed, and, if feasible, incorporated into the authors' models.

While this may be true in some contexts, the limited data on this topic (namely Worral et al. Cell Rep. 2018 and Dumont et al. EMBO J. 2020) do not agree on which chromosomes are mis-segregated more often. Worral suggested chromosomes 1-2 are particularly mis-segregated, whereas Dumont finds chromosome 3, 6, X are the highest. These differences may be explained by a context-dependent effects that depend on the model and mechanism of mis-segregation. Worral uses nocodazole washout to generate merotelics whereas Dumont gets mis-segregation through depleting CENP-A. It is unknown which if these mechanisms, if either, is representative of the mechanisms at play in human tumors so we decided to take a general approach assuming equivalent mis-segregation rates. However, we appreciate that this will be a question for other readers and we have now added this to the discussion.

10. Can the authors clarify their use of live cell imaging (e.g., in Figure 6G)? Certain apparent errors that are visible by live-cell imaging (like a lagging chromosome) can be resolved correctly and result in proper segregation. Is it appropriate to directly interfer missegregation rates as is done in this paper?

We did not perform this live cell imaging experiment. We cite these data as being kindly offered by the Kops laboratory and they correspond to the scDNAseq data for normal colon and CRC organoids from Bolhaqueiro et al. Nat Gen. 2019. We agree that chromosome mis-segregation rates cannot be directly inferred by imaging. As you say, lagging chromosomes may resolve and segregate to the correct daughter cell. The fundamental assumption is that, although not all lagging chromosomes mis-segregate, that specimens with higher rate of lagging chromosomes have higher rates of mis-segreation. Because there is no gold-standard measure of CIN in the literature to date, we feel it is necessary to show the correlation between the two and how the data from that study relates to the inferred rates in this study. We have made this clearer in the text.

11. The authors should discuss in greater detail earlier mathematical models of CIN, including PMID: 26212324, 30204765, and 12446840. How does their approach improve on this prior work?

We now provide a more detailed discussion on prior mathematical models, incorporating these and others.

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

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

1. The reviewers are concerned about the choice of statistics used for the ABC analysis. While these summary statistics indeed must contain some statistical signal on the rate of chromosomal instability, the current method of "phylogenetic reconstruction" which based on euclidian distances of chromosomal copy number variation is very indirect. The reviewers are not convinced that this is an optimal way to extract relevant information from the distribution of chromosomal copy number across cells. At least they are not clearly rooted in models of clonal evolution.

We agree that Euclidean distance reconstruction of CNVs is not necessarily the optimal method of phylogenetic reconstruction. The ideal approach would be to use maternal/paternal SNPs to reconstruct phylogeny. However, this is not possible with the data available from state-of-the art single-cell whole-genome sequencing (scWGS). The current state of the art for scWGS is 0.1X coverage so median sequencing coverage is only 10% of the genome. While this is amply sufficient to estimate CNVs across large genomic regions, it is not sufficient to reconstruct phylogenies by SNPs with currently available algorithms. The current work nevertheless makes a material advance over existing studies which use estimates of euclidean distances alone to infer CIN. We have expanded our discussion of this important limitation in the manuscript to acknowledge this limitation.

Why this particular set of summary statistics was chosen, how it relates to the quantity of interest (rate of chromosomal instability), how it compares to other possible ways to capture relevant information must be properly justified. The choice of these statistics should be explored in much more details in the manuscript with clear rationale given and caveats with this choice clearly stated.

Thank you for raising the question about selection of summary statistics. In Approximate Bayesian Computation (ABC), only a small number of summary statistics are employed—larger numbers impair the model (Beaumont et al. Genetics 162: 2025, 2002). To identify the optimal ones, Figure 5 Supplementary 3 evaluates the 9 possible combinations of summary statistics to infer mis-segregation rate, compared with experimental, consistent with what is thought to be the best approach (Csilléry, Katalin, et al. "Approximate Bayesian computation (ABC) in practice." Trends in ecology & evolution 25.7 (2010): 410-418.). Over 90% accuracy is achieved with 4 summary statistics—aneuploidy, MKV, Colless, and cherries—justifying the use of these statistics. This was not well described in the prior version of the manuscript—we have revised the manuscript to clearly describe the approach to select the summary statistics.

2. Related to point 1 – The phylogenetic reconstruction is an unusual way to summarise the statistical data and has not been rigorously assessed. The authors should elaborate on the consequences of the possibly low accuracy of the tree reconstruction itself.

As noted above, we have expanded the discussion about using CNVs to reconstruct phylogenetic trees, noting the possible low accuracy of this tree reconstruction.

3. The manuscript is difficult to follow. For example, there is a lengthy part on validating simulations that doesn't have a very specific message and could be reduced quite a lot. Please also ensure that the manuscript narrative is accessible to the general readership of eLife, including defining specialist concepts, referring to figures (and figure supplements) in the order they appear in the text. Please also consider adding additional labels to the figures themselves to aid the reader and avoid over-simplification e.g. "Tet" "Dip" could be written in full, Figure 2D is missing a key etc.

Thank you for this important critique. As you can see, we have substantively revised the manuscript to improve flow and remove or subsume needless detail in the text. We now ensure that figures and their supplements are in order and have improved labels.

Reviewer #3 (Recommendations for the authors):

I focused my attention at looking at the revised ABC analysis.

I found the new analysis a little hard to follow.

It is not clear to me precisely what the "step windows" and sliding/expanding timesteps are in Fig6-S1&2. I couldn't see this explicitly explained in the methods. (I can see that the authors are trying to optimise the simulated number of cell divisions, I presume some kind of adaptive search was done)

In response to reviews, we had attempted to demonstrate that the model was robust to the number of time steps. Sliding step analysis showed that whether steps 0-10, 10-20, 20-30 etc. were included in the prior, results were similar. Expanding step analysis showed that larger numbers of steps could be included without major impact on the results. Unfortunately, this was a complex analysis and we failed to present it clearly. We have now removed the expanding step windows for simplicity and attempted to simplify the text and display/explain this analysis more clearly. We also have added a relevant section describing this analysis in the methods.

I couldn't understand the rationale for only considering the mutation rate in these optimisations. I think that because stronger selection suppresses diversity, the strength of selection should inversely correlate with the time needed to evolve the diversity observed in the organoids.

Thank you—we have now added selective pressure to these optimizations and display this in Figure 6—figure supplement 1B (to replace the expanding step windows).

I'm concerned about the posterior distributions shown in Figure 6-S3. In most cases the mass in the posteriors is at the extremes of the prior, which suggests that the true parameter values lie outside of the range of the prior (i.e. that the priors were too narrow). This problem appears always the case for the selection coefficient and appears to be a problem for the mutation rate in some cases too (at least it seems that the mutation rate cannot be distinguished from the boundary of 0). Possibly the same issue of a to narrow a prior applies to Figure 5F (taxol treated cell lines). The concern is that the data could be much better explained by a set of parameters located outside of the parameter space considered, which would mean that the conclusions could be incorrect.

We agree that the mass of posteriors should not be close to the edge of the prior. We have revised Figure 5F (now Figure 5E) by expanding the priors to higher mis-segregation rates to address this concern. The updated prior distribution is now used in the revised article. In Figure 6C, the joint inference distributions, S increases asymptotically with low mis-segregation rate. The reason this occurs is because even with infinite selective pressure, at a given rate of chromosome mis-segregation, there will remain a low number of aneuploids in the population. Although this model provides a precise estimate of mis-segregation rate (i.e. CIN), the estimates of selective pressure therefore fall in a very wide band. We infer this to mean that most aneuploid clones are adversely selected against, consistent with prior work (Lukow et al. Dev Cell 2021). We have expanded the description of this in the text.

Calculating posterior predictive distributions (analagous to Figure 6-S4&5) would help to convince the goodness of fit.

Done—see Figure 5—figure supplement 4B.

Line 479 says that a large-ish value of the acceptance threshold (epsilon=0.05) was used to retain prior data. This seems an incorrect statement: epsilon determines how many simulations are rejected but all prior data is always used. Smaller epsilons should give greater accuracy but at the cost of more computation.

Thank you—we have revised this text.

https://doi.org/10.7554/eLife.69799.sa2

Article and author information

Author details

  1. Andrew R Lynch

    1. Carbone Cancer Center, University of Wisconsin-Madison, Madison, United States
    2. McArdle Laboratory for Cancer Research, University of Wisconsin-Madison, Madison, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0238-682X
  2. Nicholas L Arp

    Carbone Cancer Center, University of Wisconsin-Madison, Madison, United States
    Contribution
    Formal analysis, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8709-0667
  3. Amber S Zhou

    1. Carbone Cancer Center, University of Wisconsin-Madison, Madison, United States
    2. McArdle Laboratory for Cancer Research, University of Wisconsin-Madison, Madison, United States
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  4. Beth A Weaver

    1. Carbone Cancer Center, University of Wisconsin-Madison, Madison, United States
    2. McArdle Laboratory for Cancer Research, University of Wisconsin-Madison, Madison, United States
    3. Department of Cell and Regenerative Biology, University of Wisconsin, Madison, United States
    Contribution
    Funding acquisition, Project administration, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7830-3816
  5. Mark E Burkard

    1. Carbone Cancer Center, University of Wisconsin-Madison, Madison, United States
    2. McArdle Laboratory for Cancer Research, University of Wisconsin-Madison, Madison, United States
    3. Division of Hematology Medical Oncology and Palliative Care, Department of Medicine University of Wisconsin, Madison, United States
    Contribution
    Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review and editing
    For correspondence
    mburkard@wisc.edu
    Competing interests
    declares the following: Medical advisory board of Strata Oncology; Research funding from Abbvie, Genentech, Puma, Arcus, Apollomics, Loxo Oncology/Lilly, and Elevation Oncology. I hold patents on microfluidic device for drug testing, and for homologous recombination and super-resolution microscopy technologies. I declare all interests without adjudicating relationship to the published work
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4215-7722

Funding

National Cancer Institute (R01CA234904)

  • Mark E Burkard

National Institutes of Health (R01GM141068)

  • Mark E Burkard

National Cancer Institute (P30CA014520)

  • Mark E Burkard

National Cancer Institute (F31CA254247)

  • Andrew R Lynch

National Institutes of Health (T32HG002760)

  • Andrew R Lynch

National Institutes of Health (T32GM81061)

  • Andrew R Lynch

National Institutes of Health (T32GM008692)

  • Nicholas L Arp

National Institutes of Health (T32GM008688)

  • Amber S Zhou

National Institutes of Health (T32GM140935)

  • Nicholas L Arp

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

Acknowledgements

This study was supported by grants to MEB and BAW from the NCI (5R01CA234904). ARL was supported by the UW Cellular and Molecular Pathology (5 T32GM081061) and the UW Genomic Sciences Training Program (5T32HG002760) NIH training grants. NLA was supported by T32GM008692. ASZ. was supported in part by T32GM008688.Technical support comes from University of Wisconsin Carbone Cancer Center (UWCCC) Shared Resources funded by the UWCCC Support Grant P30 CA014520 – Flow Cytometry Core Facility (1S10RR025483-01), Cancer Informatics Shared Resource, Small Molecule Screening Facility. The authors thank the UW Biotechnology Center DNA Sequencing Facility for providing Illumina sequencing services. Special thanks go to Drs. Ana Bolhaqueiro, Bas Ponsioen, and Geert Kops for the provision of scDNAseq data for our analyses and to Dr. Caitlin Pepperell for valuable comments related to approximate Bayesian computation.

Senior Editor

  1. Anna Akhmanova, Utrecht University, Netherlands

Reviewing Editor

  1. Adèle L Marston, University of Edinburgh, United Kingdom

Reviewer

  1. Trevor A Graham, Queen Mary, University of London, United Kingdom

Publication history

  1. Preprint posted: April 27, 2021 (view preprint)
  2. Received: April 27, 2021
  3. Accepted: April 1, 2022
  4. Accepted Manuscript published: April 5, 2022 (version 1)
  5. Version of Record published: April 29, 2022 (version 2)

Copyright

© 2022, Lynch 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

  • 1,068
    Page views
  • 179
    Downloads
  • 0
    Citations

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

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)

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

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

  1. Andrew R Lynch
  2. Nicholas L Arp
  3. Amber S Zhou
  4. Beth A Weaver
  5. Mark E Burkard
(2022)
Quantifying chromosomal instability from intratumoral karyotype diversity using agent-based modeling and Bayesian inference
eLife 11:e69799.
https://doi.org/10.7554/eLife.69799

Further reading

    1. Cancer Biology
    Samantha C Schwager, Katherine Young ... Cynthia A Reinhart-King
    Research Article

    Cancer cell migration is highly heterogeneous, and the migratory capability of cancer cells is thought to be an indicator of metastatic potential. It is becoming clear that a cancer cell does not have to be inherently migratory to metastasize, with weakly migratory cancer cells often found to be highly metastatic. However, the mechanism through which weakly migratory cells escape from the primary tumor remains unclear. Here, utilizing phenotypically sorted highly and weakly migratory human breast cancer cells, we demonstrate that weakly migratory metastatic cells disseminate from the primary tumor via communication with stromal cells. While highly migratory cells are capable of single cell migration, weakly migratory cells rely on cell-cell signaling with fibroblasts to escape the primary tumor. Weakly migratory cells release microvesicles rich in tissue transglutaminase 2 (Tg2) which activate murine fibroblasts and lead weakly migratory cancer cell migration in vitro. These microvesicles also induce tumor stiffening and fibroblast activation in vivo and enhance the metastasis of weakly migratory cells. Our results identify microvesicles and Tg2 as potential therapeutic targets for metastasis and reveal a novel aspect of the metastatic cascade in which weakly migratory cells release microvesicles which activate fibroblasts to enhance cancer cell dissemination.

    1. Cancer Biology
    2. Immunology and Inflammation
    Lei Yang, Xichen Dong ... Zhenjun Wang
    Research Article

    Efficacy of immunotherapy is limited in patients with colorectal cancer (CRC) because high expression of tumor-derived transforming growth factor (TGF)-β pathway molecules and interferon (IFN)-stimulated genes (ISGs) promotes tumor immune evasion. Here, we identified a long noncoding RNA (lncRNA), VPS9D1-AS1, which was located in ribosomes and amplified TGF-β signaling and ISG expression. We show that high expression of VPS9D1-AS1 was negatively associated with T lymphocyte infiltration in two independent cohorts of CRC. VPS9D1-AS1 served as a scaffolding lncRNA by binding with ribosome protein S3 (RPS3) to increase the translation of TGF-β, TGFBR1, and SMAD1/5/9. VPS9D1-AS1 knockout downregulated OAS1, an ISG gene, which further reduced IFNAR1 levels in tumor cells. Conversely, tumor cells overexpressing VPS9D1-AS1 were resistant to CD8+ T cell killing and lowered IFNAR1 expression in CD8+ T cells. In a conditional overexpression mouse model, VPS9D1-AS1 enhanced tumorigenesis and suppressed the infiltration of CD8+ T cells. Treating tumor-bearing mice with antisense oligonucleotide drugs targeting VPS9D1-AS1 significantly suppressed tumor growth. Our findings indicate that the tumor-derived VPS9D1-AS1/TGF-β/ISG signaling cascade promotes tumor growth and enhances immune evasion and may thus serve as a potential therapeutic target for CRC.