Measuring the distribution of fitness effects in somatic evolution by combining clonal dynamics with dN/dS ratios

<jats:p>The distribution of fitness effects (DFE) defines how new mutations spread through an evolving population. The ratio of non-synonymous to synonymous mutations (dN/dS) has become a popular method to detect selection in somatic cells. However the link, in somatic evolution, between dN/dS values and fitness coefficients is missing. Here we present a quantitative model of somatic evolutionary dynamics that determines the selective coefficients of individual driver mutations from dN/dS estimates. We then measure the DFE for somatic mutant clones in ostensibly normal oesophagus and skin. We reveal a broad distribution of fitness effects, with the largest fitness increases found for TP53 and NOTCH1 mutants (proliferative bias 1–5%). This study provides the theoretical link between dN/dS values and selective coefficients in somatic evolution, and measures the DFE of mutations in human tissues.</jats:p>

section_content: samples from healthy individuals over time) and in rare cases of patients under active surveillance. An alternative approach is to infer selective coefficients directly from somatic genome sequencing data. Methods to identify positively-selected (driver) mutations rely on finding genes that have significantly more mutational 'hits' (typically hits are non-synonymous mutations) than would be expected by chance, after correction for factors known to influence the mutation rate across the genome 8. Conversely, negatively selected genes are expected to show a paucity of mutations 9, 10. This idea is formalised in the calculation of the dN/dS ratio -a method originally developed in molecular species evolution -that has recently been adapted for use to study somatic evolution (both cancer and normal tissue) 3, [9] [10] [11] [12] [13] [14] [15]. The intuitive idea behind dN/dS is to measure the rate of non-synonymous (dN) mutations (possibly under selection) and compare that to the rate of synonymous (dS) mutations (presumed neutral). The ratio of these two numbers, each normalised for the local sequence-specific biases in the mutation rate, putatively identifies a signature of selection: dN/dS > 1 indicating positive selection, dN/dS = 1 indicating neutral evolution and dN/dS < 1 indicating negative selection. Transforming dN/dS values to selective coefficients in somatic evolution is an unaddressed problem. dN/dS was originally developed in the context of species evolution using the Wright-Fisher process, a classical population genetics model that assumes that evolution occurs over very long timescales, which permits new mutations to fix within lineages, and also that the population size is constant, with all individuals having equal potency and nonoverlapping generations. Under the Wright-Fisher model, the dN/dS of a locus is related to its selective coefficient by the relation 16 :

= 2 1 - *+,- Where is the effective population size and the selection coefficient. However, in somatic evolution the assumptions of the Fisher-Wright model are violated. Somatic evolution is rapid and new mutations are infrequently fixed in the population 17, clonal dynamics are complex and population sizes unlikely to be constant 18. Further, the lack of recombination in somatic evolution can result in strong hitchhiking effects. In addition, since in somatic evolution the ancestral genome is known it circumvents the need to measure dN/dS across a phylogeny (a necessary step for dN/dS analysis in species evolution). Violations of some of these assumptions was previously recognised to make the interpretation of dN/dS problematic 19, 20, and consequently the relationship between selective coefficients and dN/dS values is uncertain. The size distribution of clones (called the site frequency spectrum in population genetics nomenclature) also contains information on the selective coefficients of newly arising mutations. Mathematical descriptions of the dynamics of populations of cells can make predictions on the shape of the clone size distribution under different demographic and evolutionary models 21, 22, and this approach has been used to quantify the dynamics and cell fate properties of stem cells across many tissues [23] [24] [25]. We and others have also used similar approaches to infer the evolutionary dynamics of tumours in deep sequencing data [26] [27] [28] [29]. 

To date, dN/dS analysis and the analysis of the clone size distribution have been performed independently, with conflictual results 30, 31. Here we develop the mathematical population genetics theory necessary to combine these approaches and explore how the interindividual measure of selection at a locus as provided by dN/dS values is related to the underlying cell population dynamics that generate intra-individual clone size distributions. This approach naturally accounts for the nuances in somatic evolution that can make the interpretation of dN/dS difficult. We show how this unified approach allows for greater insight into patterns of selection than either method in isolation, and importantly reveal the precise mathematical relationship between dN/dS values and selective coefficients in somatic evolution. We use this approach to infer the selective advantage of mutations in normal tissue and examine the evolutionary dynamics of cancer subclones.

section_name: Results

section_content: We present a general mathematical framework for the interpretation of frequencydependent dN/dS values in somatic evolution. First, we construct null models of the evolutionary dynamics in the absence of selection, and then augment these models to incorporate the consequences of selection. Evolutionary dynamics differ between normal tissues and cancer cells: in normal tissues maintained by stem cells, the long-term population dynamics is controlled by an approximately fixed-size set of equipotent stem cells undergoing a process of neutral competition 32, whereas in tumour growth the overall population increases over time. In each scenario, we develop a null model to predict the expected genetic diversity in the population in the absence of selection. Positive selection causes selected variants to rise to higher frequency than expected under neutral evolution (Figure 1a ), and negative selection has the opposite effect. This insight guides how we model the effects of selection (i. e diversity of non-synonymous mutations). Specifically, we defined the function (,,, ) as the expected distribution of mutations with selective (dis)advantage found at a frequency, for a given evolutionary dynamics scenario, where mutations accumulate at a rate. For the remainder of the paper we use passenger mutations to refer to those mutations that have no functional effect (s=0) and driver mutations those that have s>0. When comparing to data, driver mutations are taken as equivalent to non-synonymous mutations and passengers equivalent to synonymous mutations. The functional form of (,,, ) encapsulates the population dynamics of the system with parameter vector, which may, for example, include the growth rate of a tumour, or loss replacement rate of stem cells in normal tissue. The direct interpretation of depends on the system under question. Following the logic of the effect of selection above, for 6 > we have that:

(,, 6, )> (,,, ). 

A(B,> @,-,C)

A(B,> ?,-DE,C)

[1] We discuss the general properties of this model. Firstly, when = 0 (neutral evolution), the numerator and denominator are equal resulting in <, <= = 1, as expected. Secondly, dN/dS increases as a function of frequency (clone size) for positive selection, and decreases as a function of for negative selection (Figure 1b ), for all (,,, ) that we consider. Thirdly, the shape of the curves predicted by the underlying population model encodes the value of the selection coefficient; for example the steepness of the increase is proportional to the selection coefficient (Figure 1C ). These observations are a natural consequence of positive selection driving selected mutations to higher frequency (Figure 1a ). Unfortunately, directly using equation [1] to measure selective coefficients from the slope of the dN/dS curve as function of frequency is often impractical. Real sequencing data often suffers from a limited number of mutations detected at any particular frequency and measurement uncertainties in these frequencies. To circumvent these issues, we introduce "interval dN/dS" (i-dN/dS) that aggregates over a frequency range to reduce the influence of these sources of noise. Interval dN/dS is defined as:

[2] Fixing the integration range [ PQR, PST ] allows for robust inference of in potentially sparse and noisy sequencing data using maximum likelihood methods (see Methods).
section_name: A general approach to integrate dN/dS and clone size distributions

section_content: In healthy tissue, only mutations that are acquired in the stem cells will persist over long times, and so we restrict our attention to these cells. Quantitative analysis of lineage tracing data has shown that the stem cell dynamics of many tissues conform to a process of population asymmetry 32. In this paradigm, under homeostasis, the loss of stem cells through differentiation is compensated by the replication of a neighbouring stem cell, thus maintaining an approximately constant number of stem cells. These dynamics are represented by the rate equations:

where SC refers to a single stem cell which divides symmetrically to produce either two stem cells or two differentiated cells (denoted as D above), is the rate of cell division per unit time, and is the probability of a symmetric divisions. The product is referred to as the loss/replacement rate. Differentiated cells will ultimately be lost from the population over long time scales. Under homeostasis, these processes should be exactly balanced with ∆= 0. With ∆≠ 0, the fate of a stem cell is 'biased', introducing positive or negative selection into the model. Previous mathematical analysis shows that this model is a good description of the clonal dynamics in the oesophagus and skin 23, 33, 34. Using the previous analytical results describing the temporal evolution of the clone distribution (see supplementary methods for detailed discussion) we derive the frequency distribution (,,, ) for oesophagus and skin as 21, 23, 35 : [4] Where E is the starting population size and T the mutation rate, which may be different for drivers ( ≠ 0) and passenger mutations ( = 0). () is a scaling factor that depends on ∆, the bias toward self-renewal, which we interpret as our selection coefficient in this system. Specifically:

+f [6] We note that at long times (large ()) equation [4] converges to a 1/ distribution for the site frequency spectrum of a fixed size population 36. () can be interpreted as the average size of a labelled clone after time, which even under homeostasis grows over time and compensates for some clones being lost due to drift. From these expressions, we can then write down a closed-form expression for i-dN/dS as a function of clone frequency (see methods) that allows for maximum likelihood estimation of parameter values (∆, ). We confirmed the accuracy of our derivation using simulations (Figure 2a ), and performed power calculations to determine the minimum number of mutations required to correctly infer the underlying population dynamics. We determined that 8 mutations per gene was sufficient to accurately recover ∆ (Figure 2b ) with accuracy increasing for higher mutation burdens (Figure 2c ). Selection advantages in histopathologically normal human oesophagus We inferred the selective advantage of driver mutations in human oesophagus using published deep sequencing data from Martincorena and colleagues 14, 37 that documents the clonal expansion of a panel of putative driver mutations in histopathologically-normal oesophageal biopsies. We used the dndscv bioinformatics tool 3 to calculate frequency-dependent dN/dS values from these data (clone size measured in fraction of mutant reads multiplied by 2mm 2 -the area of the biopsy -and assuming 5,000 stem cells per mm 2 tissue). dN/dS values varied considerably as a function of mutation frequency (Figure S1 ). We considered the average frequency-dependent dN/dS values across all genes in the panel, on a patient-by-patient basis. Our theoretical model of i-dN/dS calculated from these data fitted strikingly well (Figure S2 ). Estimates of the loss/replacement rate of the stem cell population were in the range 1. 2-5. 0 per year (Figure S2&S3 ). Inference of the selective advantage (measured in terms of the bias towards self renewal ∆) revealed an average bias of 0. 004 (0. 002 -0. 005 95% CI) per missense mutation (Figure S2 ). Nonsense mutations caused a five-fold greater bias towards self-renewal of 0. 021 (0. 008 -0. 032 95% CI) (Figure S3 ). After removal of all genes that are strongly selected, global dN/dS values on the remaining 48 genes show dN/dS of approximately 1 across the frequency range (Figure 2d ), and i-dN/dS analysis revealed somatic mutation does not associate with a proliferative bias (D=0). We then fitted the data on a gene-by-gene and patient-by-patient basis for cases where sufficient mutations were available to perform the fit (Figure 2e -g; Figure S4 ). A broad range of selective advantages were inferred (Figure S4&S5 ). Mutations in TP53 showed large biases across all patients for both missense, D=0. 057 (0. 05-0. 068 95% CI) and nonsense mutations, D=0. 094 (0. 091-0. 097 95% CI) (Figure 3a-b ). This was also true for mutations in NOTCH1 with D=0. 029 (0. 019-0. 036 95% CI) for missense and D=0. 072 (0. 034-0. 089 95% CI) for nonsense mutations. NOTCH2, PIK3CA, CREBBP and FAT1 also showed a bias toward selfproliferation in multiple patients (Figures 3a-b ), though most had a small effect on fitness (range 0. 003 -0. 029 for missense mutations and 0. 030 -0. 041 for nonsense mutations). Together these data suggest a distribution of fitness effects (DFE) characterized by many small effect mutations with few large effect mutations (Figures 3c-d ), as in seen in organismal evolution 2.
section_name: Frequency-dependent dN/dS values in stem cell populations

section_content: Martincorena and colleagues have also published data on the expansion of driver mutations in ostensibly normal human skin 18. Analyses of these data with interval dN/dS revealed a per-patient average selective advantage per mutation (again measured in terms of the bias towards self renewal ∆) of ∆=0. 001 for missense mutations and four-fold higher for ∆=0. 004 for nonsense mutations (Figures S6a-c ). Performing the analysis on a gene-by-gene basis was limited by the low detected number of mutations, and the limited frequency range (clone size range). Good fits to the data were obtainable for NOTCH1 missense mutations in patient PD18003 with fitness estimated to be ∆=0. 0149 (0. 0148-0. 0150 95% CI), and TP53 missense mutations also in patient PD18003, ∆=0. 0054 (0. 0051-0. 0058 95% CI) Figure S6. These fitness coefficients were similar to the oesophagus data. For missense mutations we were also able to produce the distribution of fitness effects across the skin cohort, which showed similar characteristics to the oesophagus data of a small number of high effect mutations and a larger number of smaller effect mutations, Figure S6f.
section_name: Driver mutation selective advantage in normal skin

section_content: We next investigated the selective advantage of driver mutations in cancer. We first investigated whether or not differences existed between dN/dS values for clonal mutations (ie truncal, present in all cells in a cancer) and subclonal mutations (present in a subset of cells in a cancer) were apparent. Using sequencing data from 2,619 cancers from TCGA that had sufficient cellularity and depth (see Methods) we calculated the mutation copy number (MCN) for each mutation and grouped mutations into subclonal, clonal and amplified across the cohort, where mutations with MCN < 1 were subclonal, MCN == 1 were clonal and MCN > 1 were amplified (Figure 4a ). We than calculated global dN/dS ratios for a panel of 198 high confidence driver genes (Methods). Across all cancers, the signal of positive selection was more pronounced for clonal mutations (Figures 4b-e ), with the highest dN/dS values found in amplified mutations 38. Subclonal mutations on the other hand demonstrated much lower dN/dS values. The same pattern was also evident in individual cancer types (Figure 4e,d & S7 ). In many cancer types (colorectal, ovarian, glioblastoma) subclonal mutations showed no evidence of subclonal selection (neutral evolution; dN/dS = 1), Figure 4e, d & Figure S7.
section_name: Clonal mutations have greater dN/dS than subclonal mutations in cancers

section_content: We applied our mathematical approach above to calculate i-dN/dS in cancer evolution. In cancer evolution (,,, ) must account for tumour growth dynamics and subclonal mutations which may rise and fall in frequency due to selection and drift. The well-studied Luria-Delbrück distribution and its extensions describes these dynamics 39. Specifically, the Luria-Delbrück distribution describes the expected number of mutational lineages at a particular frequency assuming an underlying birth-death process for individuals in the population. For neutral mutations the site frequency spectrum has a characteristic

dependence, where is the frequency of the mutations 35, 40. Hence:

where ; is the passenger mutation rate and ; is the survival probability of a lineage at division. We previously showed that in many cancers across types (approx. 30% of cases), subclonal mutations closely follow the prediction of this neutral model 26. Extensions to the classic Luria-Delbruck distribution describe the differential fitness of mutants. We defined the relative fitness advantage as the ratio of net growth rates between wildtype 'passenger' mutations ( ; ) and driver mutations ( < ) : = X @ X ?

-1 [8] s > 0 indicated positive selection while s < 0 indicated negative selection. We also defined the birth and death rates of the respective wildtype (passengers) and mutants (drivers) as ;, ;, < and <. Here, the site-frequency distribution again follows a power law but with exponent dependent on the relative fitness advantage of the mutant 35, 40 :

Here, N is the tumour population size at the time of sampling. Using these expressions (equations 7&9), we derive i-dN/dS (see Methods). The equation exhibits the same qualitative behaviour as for the stem cell model, in that dN/dS increases as a function of frequency for positive selection and decreases for negative selection (Figure 5a ). Using a simulation-based model to generate synthetic data, we confirmed the accuracy of the model by accurately recovering the inputted selection coefficient by application of the theoretical model and maximum likelihood inference (Figure 5a ). Subclonal dN/dS is strongly influenced by the ability to resolve low frequency variants. We generated synthetic tumour cohorts that modelled subclonal selection, and simulated 'perfect sensitivity' for mutation detection. In these cases, where all mutations were resolved, we measured dN/dS»1 (and hence infer a selection coefficient of 0), despite some lineages being positively selected (Figure S9 ). If only higher frequency variants were analysed, then the measured dN/dS > 1 and the correct selective coefficient is inferred (Figure S9 ). We note that at very low frequencies the detected mutations are newly arisen in the population, and so are as yet 'unfiltered' by selection. Consequently the ratio of nonsynonymous to synonymous mutations is expected to be proportional to the respective mutation rates of the two mutation types. The abundance of low frequency mutations also increases exponentially with decreasing clone frequency, and so including very lowfrequency variants 'drowns out' the effects of selection (Figure S9C ). We note that the limited sequencing data of the majority of currently available cancer genomic data means that typically only high frequency variants are detected.
section_name: Interval dN/dS for cancer

section_content: Limitations in the quality of currently available sequencing data meant that the theoretically predicted frequency dependence of dN/dS values could not be assessed in cancer genomics data (Figure S8 ). Limited sequencing depth introduces uncertainty into the determination of variant allele frequencies ("sequencing noise") which can result in incorrect classification of mutation clonality. Visual inspection of the mutation copy number histogram for TCGA data (Figure 4a ) showed a very broad dispersion of MCNs, and the resolution at lower (subclonal) frequencies was particularly poor. Issues arising due to sequencing noise are exacerbated in the setting of dN/dS analysis where pooling the data from multiple patients with different sequencing depth and purities is required. Consequently, the range of subclonal frequencies where interval dN/dS could be calculated was severely restricted. We tested whether or not looking at individual genes (rather than individual mutations) allowed for measurement of the DFE. However, the lack of recurrent subclonal mutations on a gene-by-gene basis precluded this approach. Power calculations predicted that a minimum of 30 subclonal mutations in a given gene were required to accurately fit the interval dN/dS model (Figure 5b ). This level of subclonal recurrence of individual mutations was not seen in the data: for example, the average number of subclonal mutations in TP53 per cancer type, as well as the number of subclonal VHL mutations (which has been reported to occur subclonally at an appreciable frequency 41 ) were both well below this cutoff (Figure 5B ). Consequently, large cohorts of tumours sequenced to higher depth are required to apply this approach. Aside, we note that the traditional dN/dS approach, and also our modelling framework, assumes that mutations are independent, and consequently the possibility of hitchhiking of mutations (e. g. nested driver mutations within clones) is neglected. In simulated data, we observed high mutation rates for both driver and passenger mutations led to hitchhiking being common, and subsequent obscuring of the signal of selection (Figure 5c ). In extreme cases this led to dN/dS = 1 (apparent neutral evolution) even in the presence of multiple selected lineages. For most cancers, the number of driver mutations per cancer is thought to be low (<10) 3, but nevertheless in hypermutator cancers the hitchhiking effect is likely to be common. Thus, despite hypermutator tumours tending to have fewer copy-number alterations and hence less problematic estimation of MCNs, the prevalence of hitchhiking precludes analysis of these tumours.
section_name: Currently available cancer sequencing data is insufficient to infer selective advantages

section_content: Here we have shown that the combination of dN/dS values with mutation frequency-based information provides additional quantitative insight into dynamics of somatic evolution than either method alone. Specifically, the combined approach enables direct inference of the selection coefficients of mutations in somatic tissues. 

Using this methodology we have begun the construction of the distribution of fitness effects (DFE) in somatic evolution (Figure 3c,d & Figure S6f ). In histologically normal epithelium, mutations of most genes considered showed minimal effects on fitness (near-neutral evolution), though selection coefficients for some loci, foremost NOTCH1 and TP53 were considerable (>1% and >5% respectively), and consequently the DFE has most mass close to s=0 with a long right-tail of highly-selected variants. We observed that values of selective coefficients of individual genes varies between patients, likely because of inter-patient difference in the precise location of point mutations, but potentially also because of interpatient variation in selective pressure from the microenvironment. Nevertheless, the comparative rank of per-gene fitness coefficients was broadly consistent across patients. This consistency in selective coefficients is in agreement with the observation highly recurrent gene mutations in cancer 42 and evidence of repeatability in cancer evolution 43. 

We have previously measured fitness effects in individual cancers (but were unable to ascribe fitness changes to individual genes) finding increases in growth rate in a selected clone approaching 100% in some cases 27. Care must be taken when comparing selective coefficients between normal and cancer populations, because in the former we quantify selection as tilt away from homeostasis and towards net growth of a lineage, whereas in cancer we infer the relative growth rate of a clone within the tumour as a whole. With this important caveat in mind, nevertheless the fitness increases observed in cancer appear to be much larger than for normal tissues. We hypothesise that this is because the effect of selection is weaker in expanding populations like cancer, wherein the generation of a subclonal expansion requires very large increases in fitness 44. On a cautionary note, our theoretical work shows that the clonality of mutations strongly determine the observed value of dN/dS, and so a misleading picture of the selective forces operating in a tumour (or healthy tissue) will be produced if dN/dS frequency-dependent effects are not corrected for. The accuracy of any estimate of evolutionary dynamics from dN/dS values is of course dependent of the underlying accuracy of the dN/dS measure itself, which is compromised by uncharacterised variability in the mutation rate across the genome 45 and in the uncertain pathogenicity of individual single nucleotide variants (extensions to estimate site level selection coefficients may circumvent some of these issues 46, 47 ). Finally, we note that dN/dS measures cannot elucidate evolutionary pressures in individual samples as insufficient (subclonal) mutations will be found at any individual locus. dN/dS cohort measurements are sensitive to outliers, where a few patients with high selection can drive the results 48. Other approaches, such as using the site frequency spectrum, are likely more powerful for these types of questions. Combining population genetics methods with comparative genomics is a powerful way to infer selection pressures in human somatic evolution, giving new insight into the fundamental parameters that determine evolutionary dynamics in health and disease.
section_name: Discussion

section_content: TCGA Data Processing MAF (Mutation Annotation Format) files from the Mutect2 mutation calling algorithm and copy number segmentation data for 9950 cancers from 26 cancer types were downloaded from the genomic data commons portal using the TCGAbiolinks R package 49. Cellularity and ploidy estimates derived from ASCAT were obtained from COSMIC (https://cancer. sanger. ac. uk/cosmic/download). We then filtered for >2 reads reporting the variant and >9 reads coverage at each locus in both the tumour and normal sample. We removed samples where the effective depth (defined as cellularity times depth) was < 50X and those that had likely undergone genome doubling (ploidy > 2. 5). This left 2619 samples from 17 cancer types which we deemed suitable for analysis. 

Copy number (CN) segmentations together with cellularity estimates were used to correct the variant allele frequency and produce mutation copy number estimates. We assume that the observed CN state ( yyyy ) was a combination of signals from the tumour sample and contamination from normal cells (with two copies) assuming tumour purity c. yyyy = × + 2(1 -) With this, log(R) ratios were transformed into copy number states using the following formula: = 2|2 }~•(€) -1 + • Using these corrected copy number states, mutation copy number (MCN) values were calculated. Given mutation with variant allele frequency Q, copy number Q at the locus and cellularity estimate of the tumour c, the MCN was calculated as follows: 4a ) show a dominant peak at MCN = 1 representing clonal mutations present in a single copy, confirming that the corrections we applied work as intended.
section_name: Methods

section_content: For the oesophagus and skin data we used mutation calls provided by the original studies. In the oesophagus data when a mutation was present in multiple adjacent biopsies we used the sum of the mutation frequency times the area of the biopsies (2mm 2 ) as our readout of clone size and performed the dN/dS analysis on a patient by patient basis. dN/dS calculations For calculating dN/dS ratios the dndscv R package was used which calculates both global dN/dS ratios across the whole exome or a panel of genes as well as per gene dN/dS ratios using a covariate based model to infer dN/dS values with a limited number of mutations 3. In an attempt to enrich for positive selection in some of our analysis we calculated dN/dS for a subset of 198 high confidence driver genes 50. Over or under filtering of possible germline SNPs is known to influence dN/dS values in somatic genomes 3. We previously found that mutation calls provided by TCGA are likely over stringent on filtering germline SNPs resulting in inflated dN/dS values 48. To circumvent this issue, we calculated a baseline dN/dS value by randomly selecting 1,000 genes (excluding drivers) and then running dndscv across the whole TCGA cohort, reasoning that this should on average return dN/dS = 1, and any deviation from this would be due to under/over filtering of SNPs. Repeating this procedure 50 times and then taking the mean value gave us our baseline value which we could then subtract from further dN/dS values we calculate in our analysis. To confirm this procedure produces the expected result of dN/dS = 1 in the absence of selection, we repeated the procedure and again, randomly selected 1,000 genes 100 times and then applied the correction (subtracting the calculated deviation from 1). As would be expected the mean of this distribution was dN/dS = 1, validating our approach, Figure S10. To calculate the interval dN/dS measure we took our corrected mutation frequency data and determined a low cutoff PQR based on the minimum mutation frequency. We then created a vector of frequencies PST that covered the total range of mutation frequencies and calculated dN/dS between PQR and all values of PST. This allowed us to plot dN/dS as a function of PST and fit our interval dN/dS models.
section_name: Oesophagus and skin data

section_content: We used a maximum likelihood approach to fit our models to the data. Defining the observed interval dN/dS as and the model dN/dS as ‡() =. First of all we define the residuals between the data and the model as = - ‡. Assuming that the residuals are normally distributed with mean 0 we can write down the negative log likelihood (NLL) as

where N denotes the normal probability density function. We can then find the parameters that minimize the NLL and calculate confidence intervals on these estimates using the Fisher information matrix.
section_name: Model fitting

section_content: For the stem cell model, using equations [2]- [6] in the main text, interval dN/dS is given by:

. Given that the data is in terms of area, A we made the transformation =, where is density of stem cells per mm 2, which we set to 5,000 cells /mm 2 for fitting. 

For the cancer model, interval dN/dS is given by:
section_name: Interval dN/dS models

section_content: To confirm our analytical models and investigate the influence of uncertainty in mutation frequencies due to sequencing noise and to challenge some of the underlying assumptions of our theoretical approach, we developed 2 simulation based models. The first one models cancer evolution and the second models stem cell evolution under homeostasis. For the cancer evolution model, we adapted our previously described model 27 so that mutations can be one of two types, neutral passengers or mutations that have an effect on fitness of cells (either positive or negative). We model cancer growth as a continuous time branching process. At each division, daughter cells acquire mutations with a fitness effect s at rate < and passenger mutations (which are neutral) at rate ;. This is implemented by drawing a Poisson random variable with mean given by < or ;. Fitness of passenger mutations is 0, while driver mutations have fitness advantage s, where s is defined by equation [8]. We also implemented a model where fitness was a random exponentially distributed variable with mean s. For the stem cell model we seed a population of -stem cells that then undergo loss/replacement as described by the following rate equations

As only the stem cells are long lived the differentiated cells are not explicitly modelled such that when a stem cell "differentiates" it is effectively lost from the population. As in the cancer model, during division, daughter cells acquire mutations with a fitness effect at rate < and passenger mutations at rate ;. Fitness increases the bias toward self-proliferation Δ of a stem cell lineage. Additional driver mutations do not further increase the fitness of stem cells. To calculate dN/dS across a cohort of simulated tumours or tissue biopsies we count the number of driver mutations < and the number of passenger mutations, ; and then normalize by their respective mutation rates. In our model drivers = non-synonymous and thus every driver has an effect on fitness. Then the ratio of these two numbers gives us the excess or deficit of mutations due to selection -ie the dN/dS ratio. = < / < ; / ; For the interval dN/dS we simply calculate the T between PQR and PST. To introduce uncertainty into mutation frequencies we perform a process of empirically motivated sampling to the true underlying frequency. Firstly, we specify the average depth of sequencing D, then the depth of sequencing for mutation i is given by Q = () The sampled number of read counts is then -= ( = Q, = ) And the sampled variant frequency is then -= -/ Q Code and data availability Code used for the analysis are included as a snakemake pipeline which will reproduce all the analysis and generate all the figures. Julia 51 was used for the majority the simulations and R 52 was used to analyse the data and generate the figures. Some of the analysis rely in bespoke packages written for this which are freely available under and open source licence. Code is available at github. com/marcjwilliams1/dnds-clonesize.
section_name: Simulations

section_content: Global dN/dS values in different frequency bins for patient PD31182 showing that the values depend on the frequency of mutations.
section_name: Figure S1

section_content: Model fits for all patients in the oesophagus data set. Purple points are data and red lines model fits. Fits were performed separately for missense, a and nonsense mutations, b. Each plot is annotated with the inferred bias ∆ and the R 2 value.
section_name: Figure S2

section_content: Model fits per patient and per gene per patient when there were sufficient mutations in the skin dataset. 

Points are data and lines are model fits, a-e. f shows the distributions of fitness effects for missense mutations across the cohort. There were insufficient nonsense mutations in the majority of genes to draw the equivalent plot for nonsense mutations.
section_name: Figure S6

section_content: Interval dN/dS for 192 high confidence driver mutations. We observe no patterns that are predicted by our theoretical model.
section_name: Figure S8

section_content: Generating a synthetic cohort with selection and using all mutations to infer dN/dS values shows that in this case dN/dS~1, while if we restrict our attention to high frequency variants dN/dS>1, a. Inferred selection coefficients are accurate only when using high frequency variants, b. Using our theoretical interval model equation we see that fixing $%& = 1 and taking the limit $%& → 0 results in dN/dS = 1.
section_name: Figure S9