1. Cancer Biology
  2. Computational and Systems Biology

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
https://doi.org/10.7554/eLife.69799
Share this article
Cite this article
  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
  2. Download BibTeX
  3. Download .RIS

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
Download asset Open asset
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
Download asset Open asset
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
Download asset Open asset
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.

Figure 4
Download asset Open asset
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
Download asset Open asset
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
Download asset Open asset
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-x−p|x−p
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×TOEcx−p
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...

R∼U[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...

R∼U[0,1]
Mis−segregatechromosomecifPmisseg,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...

R∼U[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>R∼U[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.
    https://www.ncbi.nlm.nih.gov/bioproject/PRJNA725515

References

  81. Website
    1. Wilensky U
    (1999) NetLogo
    Accessed March 11, 2020.
    http://ccl.northwestern.edu/netlogo/

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.

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

  • 2,193
    views
  • 274
    downloads
  • 24
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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

Share this article

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

Categories and tags

Research organism

Further reading

    1. Cancer Biology
    2. Computational and Systems Biology

    Luminal epithelial cells integrate variable responses to aging into stereotypical changes that underlie breast cancer susceptibility

    Rosalyn W Sayaman, Masaru Miyano ... Mark LaBarge
    Research Article

    Effects from aging in single cells are heterogenous, whereas at the organ- and tissue-levels aging phenotypes tend to appear as stereotypical changes. The mammary epithelium is a bilayer of two major phenotypically and functionally distinct cell lineages: luminal epithelial and myoepithelial cells. Mammary luminal epithelia exhibit substantial stereotypical changes with age that merit attention because these cells are the putative cells-of-origin for breast cancers. We hypothesize that effects from aging that impinge upon maintenance of lineage fidelity increase susceptibility to cancer initiation. We generated and analyzed transcriptomes from primary luminal epithelial and myoepithelial cells from younger <30 (y)ears old and older >55y women. In addition to age-dependent directional changes in gene expression, we observed increased transcriptional variance with age that contributed to genome-wide loss of lineage fidelity. Age-dependent variant responses were common to both lineages, whereas directional changes were almost exclusively detected in luminal epithelia and involved altered regulation of chromatin and genome organizers such as SATB1. Epithelial expression of gap junction protein GJB6 increased with age, and modulation of GJB6 expression in heterochronous co-cultures revealed that it provided a communication conduit from myoepithelial cells that drove directional change in luminal cells. Age-dependent luminal transcriptomes comprised a prominent signal that could be detected in bulk tissue during aging and transition into cancers. A machine learning classifier based on luminal-specific aging distinguished normal from cancer tissue and was highly predictive of breast cancer subtype. We speculate that luminal epithelia are the ultimate site of integration of the variant responses to aging in their surrounding tissue, and that their emergent phenotype both endows cells with the ability to become cancer-cells-of-origin and represents a biosensor that presages cancer susceptibility.

    1. Cancer Biology

    Metastasis of colon cancer requires Dickkopf-2 to generate cancer cells with Paneth cell properties

    Jae Hun Shin, Jooyoung Park ... Alfred LM Bothwell
    Research Article

    Metastasis is the leading cause of cancer-related mortality. Paneth cells provide stem cell niche factors in homeostatic conditions, but the underlying mechanisms of cancer stem cell niche development are unclear. Here, we report that Dickkopf-2 (DKK2) is essential for the generation of cancer cells with Paneth cell properties during colon cancer metastasis. Splenic injection of Dkk2 knockout (KO) cancer organoids into C57BL/6 mice resulted in a significant reduction of liver metastases. Transcriptome analysis showed reduction of Paneth cell markers such as lysozymes in KO organoids. Single-cell RNA sequencing analyses of murine metastasized colon cancer cells and patient samples identified the presence of lysozyme positive cells with Paneth cell properties including enhanced glycolysis. Further analyses of transcriptome and chromatin accessibility suggested hepatocyte nuclear factor 4 alpha (HNF4A) as a downstream target of DKK2. Chromatin immunoprecipitation followed by sequencing analysis revealed that HNF4A binds to the promoter region of Sox9, a well-known transcription factor for Paneth cell differentiation. In the liver metastatic foci, DKK2 knockout rescued HNF4A protein levels followed by reduction of lysozyme positive cancer cells. Taken together, DKK2-mediated reduction of HNF4A protein promotes the generation of lysozyme positive cancer cells with Paneth cell properties in the metastasized colon cancers.

    1. Cancer Biology
    2. Computational and Systems Biology

    Unveiling the influence of tumor and immune signatures on immune checkpoint therapy in advanced lung cancer

    Nayoung Kim, Sehhoon Park ... Myung-Ju Ahn
    Research Article

    This study investigates the variability among patients with non-small cell lung cancer (NSCLC) in their responses to immune checkpoint inhibitors (ICIs). Recognizing that patients with advanced-stage NSCLC rarely qualify for surgical interventions, it becomes crucial to identify biomarkers that influence responses to ICI therapy. We conducted an analysis of single-cell transcriptomes from 33 lung cancer biopsy samples, with a particular focus on 14 core samples taken before the initiation of palliative ICI treatment. Our objective was to link tumor and immune cell profiles with patient responses to ICI. We discovered that ICI non-responders exhibited a higher presence of CD4+ regulatory T cells, resident memory T cells, and TH17 cells. This contrasts with the diverse activated CD8+ T cells found in responders. Furthermore, tumor cells in non-responders frequently showed heightened transcriptional activity in the NF-kB and STAT3 pathways, suggesting a potential inherent resistance to ICI therapy. Through the integration of immune cell profiles and tumor molecular signatures, we achieved an discriminative power (area under the curve [AUC]) exceeding 95% in identifying patient responses to ICI treatment. These results underscore the crucial importance of the interplay between tumor and immune microenvironment, including within metastatic sites, in affecting the effectiveness of ICIs in NSCLC.