Both normal tissue development and cancer growth are driven by a branching process of cell division and mutation accumulation that leads to intra-tissue genetic heterogeneity. However, quantifying somatic evolution in humans remains challenging. Here, we show that multi-sample genomic data from a single time point of normal and cancer tissues contains information on single-cell divisions. We present a new theoretical framework that, applied to whole-genome sequencing data of healthy tissue and cancer, allows inferring the mutation rate and the cell survival/death rate per division. On average, we found that cells accumulate 1.14 mutations per cell division in healthy haematopoiesis and 1.37 mutations per division in brain development. In both tissues, cell survival was maximal during early development. Analysis of 131 biopsies from 16 tumours showed 4 to 100 times increased mutation rates compared to healthy development and substantial inter-patient variation of cell survival/death rates.

Most cells in human tissues have a limited life span and need to be replenished for tissues to remain functional^1,2,3. This cell turnover leads to somatic evolution, with cells accumulating mutations upon which selection may act^4,5. Inter- and intra-tumour genetic heterogeneity^6,7 as well as treatment resistance^8,9 are now understood to be consequences of somatic evolutionary processes. Recent studies demonstrate somatic evolution in healthy non-cancerous tissues throughout live^{10,11,12,13,14}. Normal brain cells carry hundreds of mutations weeks after conception¹² and normal skin or esophagus cells accumulate hundreds of cancer driver mutations during adulthood^10,11.

These observations call for a better quantitative understanding of the somatic evolutionary forces in both cancerous and healthy tissues¹⁵. However, unlike species evolution, for which a timed fossil record exists^16,17, the lack of sequential human data over time due to ethical and technical limitations is a major obstacle. Furthermore, some evolutionary forces are difficult to measure. For example, the mutational burden in a tissue is the combined effect of per-cell mutation and per-cell survival rates, which remain hidden in sequencing data^18,19 (Fig. 1). Currently, we cannot independently infer these two for somatic evolution fundamental quantities from single time point sequencing data.

a Each of the seven spatially separated tissue samples (in grey) consists of thousands to millions of cells that descended from a single most recent common ancestor (MRCA) cell. The genomic make-up of the single ancestral cell is described by the mutations clonal to the bulk sample. Those appear at high variant allele frequency in the sample (bottom-left panel, in purple). The intersection of mutations in any two bulk MRCA cells corresponds to the genomic profile of another more ancestral cell. This process continues back in time until the MRCA cell of all the sampled cells is reached. b The level of genomic variation within a growing tissue (e.g. development or cancer) is the direct consequence of mutation accumulation during cell divisions, leading to a branching structure. Importantly, the most fundamental parameters, the mutation rate μ and survival rate β of cells per division that drive this process are not directly observable. c Mutation rate per division μ and cell survival rate β leave identifiable fingerprints in the observable patterns of genetic heterogeneity within a tissue. Cell divisions occur in increments of natural numbers and thus the mutational distance between any two ancestral cells is a multiple of the mutation rate μ.

Here, we show that multiple bulk or single-cell sequencing from the same patient contain recoverable information on these important quantities that can be recovered with evolutionary theory. This allows inferring in vivo cell mutation and cell survival rates in tissues of individual humans from single time point sequencing data.

We draw our inferences by defining and quantifying the distribution of mutational distances amongst multiple samples. We first discuss the required theoretical considerations and derive an analytical expression for the expected distribution of mutational distances from multi-sample sequencing data. We introduce a Bayesian sampling framework based on the mutational distance distribution, allowing us to disentangle mutation rates per cell division and cell survival/death rates. We apply this framework to whole-genome single-cell sequencing data of haematopoiesis and brain tissue and measure both evolutionary parameters during early development. Finally, we utilise multi-sample sequencing data on 16 tumours to infer patient specific evolutionary parameters in human cancers.

Results

The distribution of mutational distances

All cells in a human tissue must have descended from a most recent common ancestor cell (MRCA) that existed briefly during early development. Similarly, all cells in a sample of a tissue must have descended from a (different) MRCA that was present in that tissue at an earlier time (Fig. 1a). Mutations found in all cells of the sample (clonal mutations) were present in this MRCA. If we take multiple samples of the same tissue, we can reconstruct the mutational profile (all mutations carried by a single cell) of multiple ancestral cells (Fig. 1a). Typically, these ancestral cells differ in their exact mutational profile between one another, because mutations inevitably accumulate differently in distinct lineages (Fig. 1b). We use the differences of the mutational profiles between ancestral cells to construct the distribution of mutational distances. We define a mutational distance as the number of mutations different between any two ancestral cells (Fig. 1c). In the language of set theory, if ancestral cell 1 carries a set of mutations A and ancestral cell 2 carries a set of mutations B, then by definition, both cells must have coalesced from an earlier ancestral cell (Fig. 1a). The mutational profile of this cell is given by the intersection A ∩ B. This allows us to construct two mutational distances given by

y1=|A∖(A∩B)|andy2=|B∖(A∩B)|$$y_1=\left|A\setminus \left(A\cap B\right)\right|\mathrm{a}\mathrm{n}\mathrm{d}{y}_2=\left|B\setminus \left(A\cap B\right)\right|$$

                    (1)               

This process can be iterated for increasing combinations of samples per tumour.

We now turn to quantitative expressions for the expected distribution of mutational distances P(y). In a single division, the probability of a cell to acquire X novel mutations follows a Poisson distribution

P(X)=(μL)XX!e−μL⋅$$P\left(X\right)=\frac{(\mu L{)}^X}{X!}{e}^{-\mu L}\cdot$$

                    (2)               

Here, μ is the mutation rate (in units of base pairs per cell division) and L the size of the sequenced genome. Throughout the paper, we assume a constant mutation rate and do not consider more punctuated catastrophic events or mutational bursts. Distances between cells of a lineage may arise from more than a single cell division. Instead, double, triple and higher modes of cell division contribute to the distribution of mutational distances of multiple samples. In general, a cell accumulates X₁ + X₂ +…+ X_n number of novel mutations after n divisions, which is again Poisson distributed.

In addition, we must account for cell death or differentiation, leading to lineage loss. We therefore introduce a probability β of having two surviving lineages after a cell division and a probability 1 – β of a single surviving lineage (cell death). We can split the total of n cell divisions into r divisions that result in two surviving lineages (branching divisions) and m divisions with only a single surviving lineage (non-branching divisions). The number of non-branching events m is again a random variable, which follows a Negative Binomial distribution

P(m|r)=(r+m−1r−1)βr(1−β)m.$$P\left(m|r\right)=\left(\begin{array}{c}r+m-1\\ r-1\end{array}\right){\beta }^r(1-\beta {)}^m.$$

                    (3)               

The number of mutations acquired between two branching divisions depends jointly on the Poisson distributed number of mutations and the Negative binomial distributed number of non-branching divisions m. Formally, we can write for the total number of mutations between two branching divisions

Y=∑i=1mXi.$$Y=\sum _{i=1}^m{X}_i.$$

                    (4)               

Equation (4) is a random sum of random variables and different combinations of X and m imply the same mutational burden Y within a single cell lineage. Intuitively, a measured mutational burden in a single lineage can result from either many non-branching divisions with a low mutation rate or, alternatively a few non-branching divisions with high mutation rate. The mutational burden of a single sample is insufficient to disentangle per-cell mutation and per-cell survival/death rates.

We therefore turn to the number of mutations different between ancestral cells. Suppose two ancestral cells are separated by r branching divisions. Following from Eq. (4), we can calculate the probability distribution of the number of acquired mutations P(y|r) after r branching divisions

P(y|r)=∑i=r∞(i−1r−1)βr(1−β)i−re−iμL(iμL)yy!.$$P(y|r)=\sum _{i=r}^{\mathrm{\infty }}\left(\begin{array}{c}i-1\\ r-1\end{array}\right){\beta }^r(1-\beta {)}^{i-r}{e}^{-i\mu L}\frac{(i\mu L{)}^y}{y!}.$$

                    (5)               

Here the sum starts at r, as we need to have at least r branching divisions and runs to infinity as in principal infinitely many non-branching divisions can occur (with vanishingly low probability). Finally, we need the expected distribution of branching divisions P(r) in a growing population of cells, which follow from coalescence theory^20,21,22. For a growing population, e.g. human tissues during early development or cancer growth, we find

P(r)=exp(−e−β(r+1)β)−exp(−e−βrβ)1−exp(−e−ββ).$$P(r)=\frac{\exp \left(-\frac{e^{-\beta (r+1)}}{\beta }\right)-\exp \left(-\frac{e^{-\beta r}}{\beta }\right)}{1-\exp \left(-\frac{e^{-\beta }}{\beta }\right)}.$$

                    (6)               

We provide a more detailed derivation in the Methods. Combining Eqs. (5) and (6) we arrive at the final expression for the expected distribution of mutational distances in an exponentially growing population

P(y)=∑r=1∞∑i=r∞P(r)(i−1r−1)βr(1−β)i−re−iμL(iμL)yy!.$$P(y)=\sum _{r=1}^{\mathrm{\infty }}\sum _{i=r}^{\mathrm{\infty }}P(r)\left(\begin{array}{c}i-1\\ r-1\end{array}\right){\beta }^r(1-\beta {)}^{i-r}{e}^{-i\mu L}\frac{(i\mu L{)}^y}{y!}.$$

                    (7)               

The two evolutionary parameters of interest here, the mutation rate per cell division μ and the cell survival rate β, disentangle in Eq. (7). There are approximately four possible regimes for the distribution of mutational distances, discriminated by uni- or multimodality determined by combinations of small or large μ and β. In Fig. 2a we show four representative realisations of Eq. (7). The distribution of mutational distance is unimodal for sufficiently small mutation rate μ (bottom panels in Fig. 2a) with a single peak at the mean mutational distance μL. The per-cell survival probability β determines the weight of the distribution towards larger distances. For β = 1 the distribution is sharply located around the mean mutation rate. However, for smaller β more weight is given to larger distances and the distribution gets a fat tail. The same is true for the case of high mutation rate μ, except the distribution is multi-modal with peaks separated by multiples of the mean mutational distance μL (Fig. 2a). Again, β determines the weight to higher mutational distances with lower β causing a distribution with a long oscillating tail (top right panel in Fig. 2a). Note, the y-axes in Fig. 2a correspond to the probabilities of observing certain mutational distances. Lower probabilities require a higher resolution and therefore more sampling to resolve the exact shape of the distribution. In practice, the distribution of mutational distances is easiest to recover from data with low μ and high β (fewest number of tissue samples required), whereas most samples are required for high μ and low β (top right panel in Fig. 2a).

a The quantised nature of cell divisions leads to a characteristic predicted distribution of mutational distances across cell lineages. The shape of the distribution depends on the exact values of μ and β. Roughly four different scenarios of combinations of small and large μ and β are possible. They influence the shape of the distribution differently and thus constructing the distribution of mutational distances allows disentangling the mutation rate μ and cell survival rate β. b Spatial stochastic simulations confirm the ability of mutational distance distributions to disentangle mutation and lineage expansion rates (red area shows the spatial spread of a subclonal mutation). Dots show mutational distances inferred from 200 samples of a single stochastic computer simulation (μ = 20, β = 0.95), the dashed line is the predicted distribution based on our Eq.(7). c A Monte Carlo Markov Chain inference framework based on mutational distance distributions reliably identifies mutation and lineage expansion rates in simulations of spatial and stochastically growing tissues (two-dimensional spatial stochastic simulations, μ: Spearman Rho = 0.98, p = 4 × 10^–23; β: Spearman Rho = 0.93, p = 8 × 10^–16, Relative error: η_μ = 0.056, η_β = 0.045).

Computational validation and MCMC inference framework

We implemented stochastic spatial simulations of mutation accumulation in growing tissues using previously published code²³. Briefly, cell birth and death on a two- or three-dimensional grid was simulated using a Gillespie algorithm²⁴. During division, cells accumulate a number of new mutations drawn from a Poisson distribution. Simulations were stopped when the tissue reached ~1 million cells. This allowed us to take samples (either single cells or bulks) and construct all pairwise mutational distances of all ancestral cell lineages detectable in the samples. In Fig. 2b we show an example of the mutational distance distribution derived from 200 samples of a stochastic simulation (dots) compared to the theoretical prediction (dashed line).

We want to infer the microscopic evolutionary parameters μ and β given a measured distribution of mutational distances. This can be done by Markov chain Monte Carlo methods (MCMC). We implemented a standard Metropolis-Hastings algorithm. In brief, a random pair of parameters μ and β is drawn from uninformed uniform distributions and the likelihood of the model parameters given the data is calculated. The new set of parameters is accepted with a probability proportional to the likelihood ratio of the new and old parameter set (see Methods for more details). This framework recovers the true underlying parameters from stochastic simulations (Fig. 2c and Supplementary Figs. 17–21).

In vivo mutation and cell survival rate inference in healthy haematopoiesis during early development

We discuss the in vivo mutation accumulation in healthy haematopoiesis during early development as a first application. The cell population is growing and we expect a low mutation rate and a high per-cell survival rate during the development of early haematopoiesis^13,25. In a recent study, Lee-Six and colleagues¹³ sequenced the genome of 89 healthy haematopoietic stem cells of a single 59-year-old man and subsequently constructed the phylogeny of healthy haematopoiesis. They estimated the per-cell mutation rate to be 1.2 mutations per genome per division during early development assuming perfect cell doublings. Using the same data we construct the pairwise mutational distances of all ancestral cells limited to the 20 earliest branching events. The resulting distribution of mutational distances is shown in Fig. 3a. We then use the same MCMC framework discussed above to jointly infer the mutation and cell survival rate. The MCMC algorithm rapidly converges to a fixed set of parameters (Supplementary Fig. 17). In Fig. 3a, b we show the posterior parameter distributions after an initial burn in phase of 200 MCMC steps. In agreement with Lee-Six and colleagues, we find a mutation rate of μ=1.14+0.12−0.24$\mu ={1.14}_{-0.24}^{+0.12}$ mutations per genome per division (shown is the medium mutation rate per bp/cell-division and 95% credibility intervals inferred from the MCMC posterior parameter distribution), which corresponds to a mutation rate of μ = 3.9 × 10^–10 base pairs/division (assuming 3 × 10⁹ bp in the human genome). Furthermore, we infer a per-cell survival rate of β=0.96+0.038−0.102$\beta ={0.96}_{-0.102}^{+0.038}$, independently confirming the original assumption of almost perfect cell doubling during early development¹³.

a Mutational distance distribution inferred from 89 whole-genome sequenced healthy haematopoietic stem cells from ref. ¹³ (black dots), and best theoretical fit (grey line). Posterior parameter distribution of the MCMC inference for b the mutation rate per cell division (μL=1.14+0.12−0.24$\mu L={1.14}_{-0.24}^{+0.12}$ mutations per genome per cell division) and c the cell survival rate (β=0.96+0.038−0.102$\beta ={0.96}_{-0.102}^{+0.038}$). Median point estimates and 95% credibility intervals were taken from the posterior parameter distributions. The inferred mutation rate per cell division agrees with the original estimation of 1.2 mutations per cell division. Furthermore, our joined inference of mutation and cell survival rate confirms the original assumption of no cell death during early development of haematopoiesis.

In vivo mutation and cell survival rate inference in single neurons during development

In a recent publication, Bae et al.¹² collected single neurons from three fetuses 15 to 21 weeks post conception. Cells were expanded in culture and the whole-genome was sequenced. Here we focus on the case where 14 whole-genome sequenced single neurons were available (one fetus 17 weeks after conception). Again, we inferred all pairwise mutational differences, constructed the corresponding distribution of mutational distances (Fig. 4a) and used our MCMC framework for joint parameter estimates. The MCMC converges rapidly and we find sharply localised posterior distributions for the mutation and cell survival rate. We infer a median mutation rate of μ=1.37+0.1−0.1$\mu ={1.37}_{-0.1}^{+0.1}$ mutations per genome per division (corresponding to a mutation rate of μ = 4.6 × 10^–10 base pair/division) and a per-cell survival rate of β=0.998+0.002−0.01$\beta ={0.998}_{-0.01}^{+0.002}$. This inference agrees with Bae and colleagues original estimate of 1.3 mutations per genome per division based on a weighted average of all three fetuses, again assuming no cell death during early development. It also agrees with estimates of 1.2 mutations per division from de novo SNVs in familial trios²⁶. The almost identical mutation rates in haematopoietic and brain tissue during early development may not be surprising. We would expect the DNA duplication and repair machinery to be stable across tissues during early development. It may even remain stable throughout live, as suggested by the linear rate of mutation accumulation with age across individuals^27,28,29.

a Mutational distance distribution inferred from 14 whole-genome sequenced single neurons from ref. ¹² derived from one fetus (17 weeks past conception) (black dots), and best theoretical fit (grey line). MCMC inference for b the mutation rate per cell division (μL=1.37+0.1−0.1$\mu L={1.37}_{-0.1}^{+0.1}$ mutations per genome per cell division) and c the per-cell survival rate (β=0.998+0.002−0.01$\beta ={0.998}_{-0.01}^{+0.002}$). Median point estimates and 95% credibility intervals were taken from the posterior parameter distributions. The inferred mutation rate per cell division agrees with the original estimation of 1.3 mutations per cell division. Furthermore, our joined inference of mutation and cell survival rate confirms the original assumption of no cell death during early brain development.

In vivo mutation and cell survival rates in human tumours

We then investigated the per-cell mutation and survival rates in individual tumours. We analysed whole genome or exome sequencing of 131 biopsies from 16 tumours comprised of 1 colon adenoma, 7 colon carcinomas, 5 clear cell renal carcinomas and 2 lung squamous cell carcinomas (Table 1). When whole-genome sequencing was available, the mutational load was sufficient to apply the inference framework to each chromosome separately (Fig. 5 and Supplementary Figs. 1–9). The analysis was restricted to regions of chromosomes with same copy number profile in all samples of a tumour and inferences were normalised by copy-number and genome content. The resolution to infer the distribution of mutational distances from tumours was lower compared to healthy haematopoiesis or brain during development. Nevertheless, in most cases, the reconstructed distributions recover important features of the theoretical distribution (Supplementary Figs. 1–9 and 14). We found that mutation rates per cell division were 4–100 times higher in tumours compared to healthy tissue, ranging from 2.91 × 10^–9 (bp/division) in the colon adenoma to 53 × 10^–9 (bp/division) in one lung squamous cell carcinoma (Table 1). Mutation rates differ significantly between patients but not across chromosomes of the same patient (Supplementary Figs. 11 and 12). Overall this suggests important differences in mutation accumulation at the single cell level between tumours and is in agreement with recent experimental in vitro single cell mutation rate inferences^29,30.

a–c Examples of the mutational distance distribution on single chromosomes for three different colorectal carcinomas for which 6, 7 and 9 multi-region bulk samples were sequenced at whole-genome resolution (dots = data, dashed line = theoretical prediction based on MCMC parameter estimates—see insets). The distribution of mutational distances differs between patients, with Patient 04 (MSI Microsatellite Instability) showing one order of magnitude larger mutational distances. d–f Per-cell mutation rate per chromosome separated by trinucleotide mutational signature. Results are consistent across chromosomes, as expected (Methods). g–i The mean overall mutation rates are (μ02=(1.0+0.46−0.07)×10−8,μ03=(2.4+0.41−0.19)×10−8andμ04=(3.1+0.35−0.12)×10−8bp/division$({\mu }_{02}=({1.0}_{-0.07}^{+0.46})\times {10}^{-8},{\mu }_{03}=({2.4}_{-0.19}^{+0.41})\times {10}^{-8}\mathrm{a}\mathrm{n}\mathrm{d}{\mu }_{04}=\left({3.1}_{-0.12}^{+0.35}\right)\times {10}^{-8}\mathrm{b}\mathrm{p}/\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}$, dashed lines), 20–60 times higher compared to healthy somatic cells. Patient 04 is MSI+ highlighted by signature 6. j–l Estimates of per-cell survival rates per chromosome are consistent across chromosomes of the same patient (Median: β02=0.51+0.05−0.05,β03=0.65+0.02−0.02,β04=0.34+0.01−0.01${\beta }_{02}={0.51}_{-0.05}^{+0.05},{\beta }_{03}={0.65}_{-0.02}^{+0.02},{\beta }_{04}={0.34}_{-0.01}^{+0.01}$), but vary considerably between patients (Supplementary Fig. 12).

To further unravel the underlying differences in mutation accumulation during tumour growth, we decomposed somatic mutations into the most prevalent trinucleotide mutational signatures³¹ for three whole-genome sequenced colorectal carcinomas and inferred per-cell mutation and per-cell survival rates per signature in each chromosome (Fig. 5). Again, we find significant differences between patients (Supplementary Fig. 13), further supporting inter-tumour differences of mutation accumulation at the single cell level.

The inter-patient variation of the cell survival rate was evident. Whereas in healthy tissue almost all cells survive during development, in tumours cell survival rates vary between 0.34 in one MSI+ colon carcinoma and 0.86 in one renal cell carcinoma (Table 1). Again, per-cell survival rates were overall consistent if inferred from chromosomes of individuals, but varied significantly between patients (Fig. 6 and Supplementary Fig. 12). The underlying reasons for this inter-patient variation may be cell intrinsic and/or extrinsic, e.g. high cell death due to genomic instability, high mutational burden or immune surveillance. It will be of high interest to further unravel these differences on a patient specific basis in future studies. It should be noted that the inferred cell survival rates are high compared to previous estimates^32,33. This is a direct consequence of the joint inference of mutation and cell survival rates that was not possible in earlier work.

For each of the 16 tumours analysed we plot the per-cell mutation rate versus the per-cell survival rate. Median estimates and 95% credibility intervals for the mutation and cell survival rate are derived from the MCMC inferences as described in the main text. Dashed lines correspond to values of healthy tissue (μ_h = 1 × 10⁻⁹, β_h = 1/3). White background corresponds to β values that allow for growing cell populations as β = 1/3 corresponds to stable (homeostatic) populations. Shaded area describes values of β that would lead to population extinction. Most cancers scatter across the map, indicating extensive inter-patient heterogeneity.