Clonal dynamics limits detection of selection in tumour xenograft CRISPR/Cas9 screens – Cancer Gene Therapy

Selection of suitable cell lines for in vivo CRISPR/Cas9 screens

We hypothesised that a cell line with a high proportion of tumour-initiating cells would be better suited for in vivo CRISPR/Cas9 screens as it would undergo a reduced population bottleneck due to a larger founder population and therefore be less sensitive to stochastic effects. We evaluated 22 UT-SCC head and neck squamous cell carcinoma cell lines and an HCT116 colon carcinoma subline (HCT116/54C) by inoculating immunodeficient mice subcutaneously with different numbers of cells to identify which cell lines could initiate tumours with few cells. In initial experiments, only 11 of 23 cell lines grew as tumour xenografts within 100 days of inoculation of 5 × 10⁶ cells into NIH-III mice (Fig. 1A). Four of these cell lines (UT-SCC-1B, UT-SCC-54B, UT-SCC-74A, and UT-SCC-76A) had low tumour take rates and slow growth with ≤50% of tumours established by day 90, while two others (UT-SCC-1A and UT-SCC-42B) had frequent ulceration of skin over the tumours.

The five cell lines with suitable growth kinetics (UT-SCC-16A, -74B, -110B, -126A and HCT116/54C) were selected to conduct limiting dilution assays in NIH-III and NSG mice to estimate the TD₅₀. We initially screened all cell lines for growth as tumour xenograft models in NIH-III mice by inoculating 5 × 10⁶ cells per flank into 3–5 mice, followed by 10⁶ and 10⁵ cells for those with a high take rate (>75%). For tumour models that grew as xenografts at 10⁵ or 10⁶ cells, we proceeded to test in NSG mice, which were expected to be more conducive to tumour growth, at 10⁵, 10⁴ and 10³ cells until fewer than 50% of tumours formed. Tumours were formed in 100% of NIH-III mice at ≥10⁶ for all cell lines, but the establishment rate fell with cell inocula ≤ 10⁴ in NSG mice (Fig. 1B). As expected, tumour growth took longer on average with smaller inocula (Fig. 1C, D) and tumour establishment rate was higher and faster in NSG mice than NIH-III mice (Fig. 1E). The lowest TD₅₀ of approximately 1100 cells was estimated for HCT116/54C and UT-SCC-126A, followed by UT-SCC-16A, UT-SCC-110B and UT-SCC-74B (Table 1). Despite its low TD₅₀, UT-SCC-126A was not investigated further because it caused considerable weight loss in some animals (e.g. three of four NSG mice inoculated bilaterally with 10⁴ UT-SCC-126A required culling early due to ~20% bodyweight loss).

Pilot xenograft experiments with HCT116/54C and UT-SCC-74B GeCKOv2 libraries

We next compared tumour xenografts of GeCKOv2 libraries of HCT116/54C (joint lowest TD₅₀) and UT-SCC-74B (a second model with a higher TD₅₀) cell lines to ascertain whether the TD₅₀ was consistent with model suitability for in vivo screens. Two mouse strains—NSG and NIH-III mice—were compared by inoculating 10⁷ cells and harvesting tumours when the smallest in each group reached ~250 mm³. HCT116/54C GeCKO tumours grew faster in NSG mice than NIH-III mice and faster than UT-SCC-74B GeCKO tumours, which grew equally in both mouse strains (Fig. 2A).

A Growth curves for HCT116/54C and UT-SCC-74B tumours in NSG and NIH-III mice. Bars represent the mean and SEM of 3–6 tumours. B Percentage of all or NTC gRNAs detected for plasmid, HCT116/54C GeCKO cell inoculum, HCT116/54C GeCKO tumours in NSG mice, HCT116/54C GeCKO tumours in NIH-III mice and UT-SCC-74B GeCKO cell inoculum and tumours in NSG mice. Replicates for plasmid and cell samples are PCR replicates from the same gDNA sample (n = 2 each); tumour replicates are individual tumours in separate hosts (n = 3 per group). The bar indicates mean ± SEM. C Lorenz curves to show the distribution of all gRNA read counts. The curve for summed PCR replicates (plasmid and cell samples) or median tumours (by Hoover index) in a group is highlighted. The black line is the line of equality (plas: plasmid; HCT: HCT116/54C; 74: UT-SCC-74B; NIH: NIH-III). D Hoover index of each sample as a measure of inequality of the read counts for all and NTC gRNAs. E Letter-value plot of log normalised counts per million (log-ncpm) for all gRNAs in plasmid, cell inocula and GeCKO tumour samples (zeroes given pseudocount of 0.5). The white line inside the black boxes indicates the median (M), and the black boxes indicate the upper/lower quartiles (F; fourths). The next smallest boxes indicate the upper/lower eighths (E) and so on. Counts for plasmid and cell samples were summed from two PCR replicates prior to normalisation. The remaining outliers are shown as open circles.

Almost all gRNAs were detected in the GeCKOv2 plasmid library (99.98%), the HCT116/54C GeCKO library used for inoculation (99.7%), and cells cultured in vitro for a further 12 d (99.3%; Fig. 2B; Table S1). We observed high representation in HCT116/54C GeCKO tumours in NSG hosts with 92.2–96.9% gRNAs detected, with slightly lower representation in NIH-III hosts at 88.3–93.5% (Fig. 2B). The vast majority of gRNAs were also detected in UT-SCC-74B GeCKO cells (98.5%) but UT-SCC-74B GeCKO tumours in NSG mice had considerably lower representation of 66.0–68.5%, with many detected with only a single read (Fig. 2B; Table S1). In considering only NTC gRNAs, >98% were detected in all HCT116/54C GeCKO tumours but only ~86% in UT-SCC-74B GeCKO tumours (Fig. 2B; Table S2). This suggests that cells receiving non-targeting gRNAs had a growth or survival advantage within the tumours compared to cells receiving most gRNAs targeting a gene.

The decreases in representation in the tumour samples were accompanied by disproportionate increases in the read counts of the top-ranked gRNAs (i.e. those with the highest counts; Table S1). The unequal read count distributions in the tumour samples were also evident in the Lorenz curves of the gRNA count data (Fig. 2C), with the UT-SCC-74B GeCKO/NSG tumours showing extremely unequal read counts (top 1% of gRNAs accounting for 76–78% of reads). The NTC gRNAs showed only slightly smaller read count inequality than gene-targeting gRNAs, suggesting that the major driver is stochastic effects rather than changes in fitness (Table S2, Fig. S1A). As a metric of inequality, we calculated the Hoover index, which can be interpreted as the proportion of reads needed to be redistributed to form a uniform distribution (Fig. 2D; Tables S1 and S2). The Hoover indices followed the inverse pattern of gRNA representation (compare Fig. 2D to B). The inverse relationship suggests that the apparent loss of gRNA representation is at least partly driven by difficulty sampling low-frequency gRNAs in highly unequal gRNA distributions.

The skewed distribution of read counts, especially in UT-SCC-74B tumours, meant that the data were zero-inflated at one end and contained very high counts for a few gRNAs at the other end. We found that most count normalisation methods, including the commonly used TMM [40] and RLE [41] methods, as well as the GMPR method [42] developed for zero-inflated microbiome data, performed poorly with these highly unequal count distributions containing many sampling zeroes. Therefore, we devised a new normalisation procedure, which we call MPM, for our data to account for its unique characteristics (Supplemental Methods and Results; Figs. S4–S7). MPM is a generalisation of the GMPR method, but a key difference is the treatment of zeroes—these are excluded completely in GMPR but are substituted in MPM for pseudocounts when a zero is present in one sample during a pairwise comparison (double zeroes are excluded). This treatment is more consistent with all zero counts being modelled as sampling zeroes, as suggested in previous work [43]. The exclusion of all zeroes in GMPR results in the loss of information that a count is lower in one sample than another sample and can bias library size estimates upwards; MPM using an appropriate pseudocount value performs better in this regard (see Supplemental Results).

Following MPM normalisation, each sample had similar median normalised gRNA counts (Fig. 2E; Fig. S1B; Fig. S9; Tables S1 and S2). This data representation again reveals that the loss of representation from the plasmid to the cell libraries and the tumours is accompanied by wider, more skewed gRNA distribution, with the upper tails becoming progressively more extreme.

Increased read count inequality following subsequent growth of HCT116/54C GeCKO tumours

Using the best combination of cell line and mouse host, we evaluated the feasibility of conducting an in vivo screen with this model. HCT116/54C GeCKO tumours (n = 36) in female NSG mice were allowed to grow until the median tumour size reached ~250 mm³ then the mice were randomised into four groups: one (A) in which tumours were harvested immediately, a mock treatment group (B) and two drug treatment groups (C with 6-thioguanine, 6-TG and D with evofosfamide). Tumours in B-D were allowed to grow until the largest tumour in each group reached ethical size limits, at which point all tumours in the group were collected (day 38, 42 and 43 post inoculation for groups B, C and D, respectively). There was a small non-significant decrease in tumour growth for the two drug treatment groups compared to controls (Fig. 3A). Representation of gRNAs in the group A tumours was high (90.1 ± 3.5%; mean ± SD) but was lower than in the similarly-sized tumours in the pilot study (94.8 ± 2.8%), most likely reflecting that the tumours grew slightly slower than in the pilot study, taking 14 days to reach an average size of 230 mm³ compared to 11 days to reach an average of 340 mm³. Moreover, there was a substantial decrease in representation of all gRNA and NTC gRNAs in tumours in groups B–D relative to the group A tumours at day 14 (Fig. 3B; Table S3). A much smaller loss of representation was observed relative to the cell inoculum by keeping the cells used for inoculation in culture for the duration of the in vivo screen (95.0% gRNAs detected after 38 additional days).

A HCT116/54C GeCKO tumour growth after the commencement of drug treatment at day 14 (arrow) (mean ± SEM; n = 8–10). B Percentage of all or NTC gRNAs detected for HCT116/54C GeCKO cells (inoculum or cultured an additional 14 or 38 days after inoculation) and HCT116/54C GeCKO tumours collected at 14 days, 38 days with no drug, 42 days with 6-thioguanine (6-TG) or 43 days with evofosfamide (evo). Replicate cell samples are PCR replicates from the same gDNA sample (n = 2 each); tumour replicates are individual tumours grown bilaterally (n = 8 for 14 d and evo; n = 10 for 38 d and 6-TG). The bar indicates mean ± SEM. C Lorenz curves to show the distribution of all gRNA read counts. The curve for summed PCR replicates (plasmid and cell samples) or median tumours (by Hoover index) in a group is highlighted. The black line is the line of equality. D Hoover index of each sample as a measure of inequality of the read counts for all and NTC gRNAs. E Letter-value plot of log normalised counts per million (log-ncpm) for all gRNAs in cell and tumour samples (zeroes given pseudocount of 0.5). Counts for plasmid and cell samples were summed from two PCR replicates prior to normalisation. Sample name refers to group, animal number and flank position (e.g. A1L represents a tumour on the left flank of mouse 1 from group A).

Similar to the small tumours in the pilot study, the decreases in gRNA representation were mirrored by increases in inequality and skew in the gRNA counts. Lorenz curves and Hoover indices show inequality was much greater in the 14-day tumours (Hoover indexes 0.54–0.62) compared to the cell inoculum (0.29) or cells cultured for an additional 14 days after inoculation (0.32) and somewhat greater than that of cells cultured for an additional 38 days after inoculation (0.50; Fig. 3C, D; Table S3). All large (38–43 day) tumours had extremely unequal gRNA distributions as shown by the Lorenz curves, Hoover indices of 0.89–0.96, and 74–96% of reads from only the top 1% of gRNAs. Letter-value plots of normalised gRNA counts similarly showed distribution becoming progressively wider going from the cell culture samples to the 14-day tumours and then large (38–43 day) tumours (Fig. 3E; Fig. S10). The latter gRNA distributions all had very long upper tails indicating a strong skew towards relatively few gRNAs having very high counts, while most others had zero or very few counts. In the most skewed large tumour sample (C4R), 17.6% of reads were from a single gRNA (Table S3). The median gRNA count was zero for all large tumours. As with the pilot study, similar but slightly smaller increases in count inequality and skew could be seen in an analysis of NTC gRNAs in the tumour samples (Fig. 3B, D; Fig. S1C, D; Table S4).

Random variation in clonal dynamics precludes gRNA-induced selection in HCT116/54C GeCKO tumours

The representation and count distribution data in HCT116/54C GeCKO tumours suggested a greater loss of targeting gRNAs compared to non-targeting controls. To more directly test this, we compared the ratio of the 90th percentile read counts of all gRNAs to the 90th percentile read counts of the NTC gRNAs. The 90th percentile was used to obtain a reasonable number of counts since the median and upper quartile were zero for several large tumour samples. If the distribution of read counts for all gRNAs and NTC gRNAs was changing uniformly during tumour growth, this ratio was expected to stay constant for all groups. Instead, there was a progressive decrease in this ratio from the cell samples to the small tumours (14 d) to the large (38–43 d) tumours (Fig. 4A). These data indicate a greater depletion of cells harbouring targeting gRNAs compared to NTC gRNAs in the tumours. A likely explanation of this effect is the overall negative selection of cells harbouring knockouts compared to control cell lines due to a general loss of fitness in these cells. Thus, although there were large random differences in growth and survival among clones, some effects of selection could be detected when averaging across a large number of gRNAs.

A Ratio of counts for the gRNA at the 90th count percentile among all gRNAs to counts for the NTC gRNA at the 90th count percentile. The bar indicates mean ± SEM. B Distribution of log₂ normalised cpm of the top 16 gRNAs (by p-value control versus 6-thioguanine) targeting genes expected to be involved in 6-thioguanine resistance in the tumour samples. The selected genes were HPRT1, NUDT5 and 23 genes in the KEGG mismatch repair geneset (6-TG = 6-thioguanine, evo = evofosfamide). Lines are the average log cpm calculated with aveLogCPM (prior count of 0.5), and statistics are for the no drug versus 6-thioguanine comparison using the edgeR exactTest function with dispersion estimated from large tumour samples only. C Variation in counts of four gRNAs across samples in the large tumour dataset. The four gRNAs represent those with the highest and 10th highest maximum log normalised cpm across all samples for all or NTC gRNAs. D Variation in counts of the six gRNAs in the library targeting CXCL2 across samples in the large tumour dataset. For C, D, the numbers above each individual targeting gRNA represent previously reported estimates of guide activity [14], with scores approaching 1 representing highly effective gRNAs and those approaching 0 representing poorly performing gRNAs.

Given the extremely high variability in the large tumour counts, we reasoned that detection of any gRNA-driven selection among the treatment groups would be rare. Unexpectedly, after filtering for low abundance gRNAs (73,421 targeting gRNAs retained), an initial analysis of the entire in vivo screen data set (including small 14 d tumours) using the edgeR exact test pipeline [37] resulted in the detection of a large number of gRNAs significantly selected either positively or negatively (>10,000 each direction; FDR < 0.1) for the comparisons between untreated vs. 6-TG and untreated vs. evofosfamide. Further investigation revealed that this was caused by the inclusion of the small tumours during the analysis. Including only the large tumours in the analysis resulted in the detection of a more realistic number of gRNAs as being selected (approximately 2000 in each direction for 6-thioguanine and 0 for evofosfamide; FDR < 0.1; Supplemental Data Files 1 and 3). This was due to a change in the dispersion estimates—common dispersion was 4.8 for the dataset including small tumours, and 8.2 for the large tumours only. In the first analysis, the inclusion of small tumours, which had lower variability in counts, effectively biased the estimated dispersion downwards relative to what would be applicable for the large tumours, resulting in failure to control type I error.

Although a number of gRNAs were detected at FDR < 0.1 in the large tumour dataset for 6-thioguanine, gene-level analysis of these gRNA-level data by RRA [38] detected no gene hits at FDR < 0.1 for either 6-TG or evofosfamide compared to untreated (Supplemental Data Files 2 and 4). To further evaluate the performance of the screen, we assessed the enrichment of gRNAs targeting well-characterised 6-TG sensitivity genes HPRT1 and NUDT5 [44], as well as those in the mismatch repair pathway [7], in treated tumours. The top gRNAs (by FDR) targeting this subset of genes were either enriched or depleted in treated tumours vs. controls rather than being predominately enriched (Fig. 4B). For example, one gRNA targeting NUDT5 was depleted (A32983 log₂ FC − 7.9, FDR = 0.078), but another was enriched (A32982 log₂ FC 6.1, FDR = 0.11). In addition, there was large variability in gRNA counts and frequent incidence of zero counts across all tumours, such that any apparent enrichment in one group versus another was largely driven by a small number of tumours having very high counts despite other samples having much lower counts (potentially zero).

The high count dispersion, together with the pattern of read counts among individual tumours for a particular gRNA, suggest a strong random component of clonal expansion within each tumour. Individual clones (with gRNAs as barcodes) potentially expand to a very high frequency in an apparently random fashion, which masks any potential effects of selection. Indeed, we found that the identities of the top gRNAs were different in each sample (Fig. S2A), and Spearman correlation analysis failed to identify any grouping of large tumours by treatment group or overall (Fig. S2B). To further demonstrate this, we plotted the read counts of the gRNAs with the highest (CXCL12 gRNA #B11840) and 10th highest (SLC18A3 gRNA #A44491) log-ncpm across all samples, as well as the NTC gRNAs with the highest (NTC gRNA #0404) and 10th highest (NTC #0161) log-ncpm across all samples. The read counts of these four gRNAs varied widely across large tumour samples (Fig. 4C). For instance, the log-ncpm for CXCL12 gRNA #B11840 in one 6-TG-treated tumour was 26.5, while zero counts were detected in another 6-TG-treated tumour. This wide variation suggests largely random expansion of clones within these tumours. In smaller tumours, the counts were less variable but were still unusually high or low in certain samples, indicating a degree of random clonal expansion was already present in these smaller tumours. To confirm this variation applied across all gRNA targeting the same gene, we also investigated all gRNAs targeting CXCL12 and observed a similar level of variation between samples as we had seen for CXCL12 gRNA #B11840 (Fig. 4D). Estimated guide activities (in the range 0–1; ref. [14]) for the CXCL12 gRNAs vary from 0.31 to 1.00. Although the CXCL12 gRNA with the most variable counts (#B11840) was the one with the highest estimated activity (1.00) while the least active guide (#A11852; 0.31) had the most dropouts, the other guides with estimated activities in the range 0.33–0.89 showed a similar pattern of count variability. This is consistent with the variation in counts being driven by stochastic effects instead of gRNA-editing effects.

Count dispersion increases with growth time and is much greater in tumour samples

To obtain a statistical measure of the level of random variation in gRNA counts in different sample groups, we used edgeR to estimate the negative binomial dispersion for different subsets of our dataset. To simplify comparisons, we determined the common dispersion, which assumes all gRNAs have the same dispersion, and took the square root of this to estimate the biological coefficient of variation (BCV; square root of common dispersion) for each subset. We performed this analysis on both technical replicates and biological replicates. Among technical replicates, the BCV across sequencing runs for the same sample was low at 0.012 (Fig. 5A), which is expected given this variation should follow a Poisson distribution. BCV of a similar magnitude was obtained for PCR replicates of the GeCKOv2 plasmid sample (0.022), but interestingly the BCV for PCR replicates of various samples of HCT116/54C GeCKO cells was greater (0.091–0.164), suggesting amplification bias during the PCR process used for gRNA readout was contributing to count variation. Several HCT116/54C GeCKO cell samples collected at different times were available, with culture periods from the initial post-transduction sample (‘time zero’) varying from 16–55 d. When all cell samples (n = 5) were analysed as a single group, the BCV was 0.489, while cell samples cultured ~30 d (28–31 d, n = 2) had a BCV of 0.392 and those cultured ~16 d (16–17 d, n = 2), 0.208, showing that dispersion was increasing with growth time. The BCV of even small tumour samples was considerably greater at 0.84–1.18, while for the large tumours, it increased further to 2.81-3.17. These increases in dispersion with growth in both the cell and tumour samples are consistent with the widening of the count distribution that occurred with increased growth time (Fig. 3E). We also estimated the count dispersion for various combinations of sample groups to illustrate the risk of underestimating dispersion when a lower dispersion group is included among higher dispersion groups. The BCV with all large tumour groups was 2.83, but when 14 d tumours were included, this decreased to 2.16, further decreasing to 1.92 if the pilot tumours were also included. Combining cell samples and all tumours led to additional decreases in the BCV estimate (Fig. 5A).

A Biological coefficient of variation (BCV) estimated in edgeR for different subsets of samples. Data points represent different sample groups, and lines represent combinations of sample groups. Seq Reps: Sequencing replicates (n = 2) of the same library (14 d tumour) produced by PCR. PCR Reps: PCR replicates of the plasmid DNA (Pls) or the same genomic DNA samples from cultured HCT116/54C GeCKO cells (Cells), n = 2 each. Cells: Replicate HCT116/54C GeCKO cell samples cultured in parallel from the time zero sample for 16–17 d (~16 d; n = 2), 28–31 d (~30 d; n = 2), or all samples combined (n = 5), including one sample cultured for 55 d; any PCR replicates were collapsed by summing prior to calculating BCV for these biological replicates. Tumours: Pilot study tumours for HCT116/54C GeCKO in NSG (NSG, n = 3) and NIH-III mice (NIH; n = 3) and UT-SCC-74B GeCKO in NSG mice (74B, n = 3). Larger study tumours collected early (14 d, n = 8), control-treated tumours (38 d, n = 10), 6-thioguanine-treated tumours (6-TG; n = 10) and evofosfamide-treated tumours (evo; n = 8). B Estimates of the median number of clones present in each sample from the pilot and larger study, based on the estimated probability of survival of cells with neutral (NTC) gRNAs, and the median number of cells per gRNA. C Estimates of BCV in simulated datasets of reduced library size by random binning and/or subsampling gRNAs and samples from small tumour data (NSG pilot + 14 d). Boxplots show the distribution of the BCV estimates for 100 simulated datasets for n = 3, n = 10 or n = 20 tumours and indicate the number of gRNAs. D Estimates of BCV in simulated datasets of reduced library size from large tumour data (untreated, 6-TG and evo; 38–43 d). The number of clones/gRNA in C and D is estimated based on the median number of clones/gRNA for a full library multiplied by the binning factor. E Power curves determined from the common dispersion estimates of simulated datasets at a log₂ effect size of 2 for small and large tumours. Circles and solid lines indicate the curve obtained from the dispersion estimate from n = 3 tumours; the median across 100 simulated datasets is plotted. The dark grey area represents the lower and upper quartiles, and the light grey area is the range. The squares and dashed line indicate the power curve when using the maximum dispersion estimate across 100 simulated datasets from n = 10 (small tumours) or n = 20 (large tumours); cln/gRNA: clones/gRNA.

Simulations show reduced count dispersion in tumour samples with reduced gRNA library sizes

The GeCKOv2 whole genome gRNA library used in this study consists of 119 461 unique gRNAs. With 10⁷ cells inoculated per tumour, there was an average of 84 cells/gRNA. However, not all these cells will survive to form a clone. To approximate the number of clones/gRNA, we used a binomial model to estimate the overall probability of a cell with neutral (NTC) gRNA surviving in a tumour to form a clone based on the pattern of gRNAs detected (at a given detection level) and the count distribution of the cell inocula. For small NSG tumours, the estimated probability of a cell surviving to form a clone was 0.034–0.10 (1 in 10–30), which implies a median of 2.3–7.0 clones/gRNA (Fig. 5B). In large tumours, the probability of clones continuing to survive decreases to 0.0022–0.0074 (1 in 134–459), implying a number of surviving clones/gRNA of 0.15–0.52 (Fig. 5B). If we assume that the number of clones/gRNA is a more important metric for coverage in in vivo screens compared to cells/gRNA, these estimates suggest that the effective in vivo coverage of our library was several orders of magnitude less than typically recommended for in vitro screens (minimum coverage of at least 100 cells/gene is recommended; ref. [45]).

These estimated numbers of clones/gRNA suggest very low effective in vivo coverage of our library. Improving coverage by increasing the number of cells inoculated per gRNA would result in greater averaging of the random effects on growth and survival, resulting in reduced count dispersion, thereby increasing statistical power to detect gRNA selection. Although there was limited scope to increase the total number of cells inoculated, the number of cells inoculated per gRNA could be increased by reducing the total library size. To assess the effects of using a smaller gRNA library for in vivo screens, we simulated smaller gRNA libraries by randomly combining the gRNA counts in our full-size library into bins of 2–1600 gRNAs, followed by subsampling gRNAs (and samples) to 100,000 (unbinned), 50,000, 25,000, 10,000, 4000, 2000, 1000, 500, 250, 125 or 62 gRNAs; these represented datasets with coverage 84–130,000 cells/gRNA, respectively, but with much lower coverage in terms of clones/gRNA: 3.4–5400 clones/gRNA in small tumours, 0.27–420 clones/gRNA in large tumours, based on our estimates of the number of clones/gRNA in our full library. We reasoned that our binned counts would mimic the random characteristics of clonal expansion in smaller libraries, as clones with different gRNAs in the larger library would be largely equivalent (in terms of random characteristics) to different clones that receive the same gRNA in a smaller library. Any signal due to the gRNAs themselves, which we have already established was small relative to the random clonal expansion effect (Fig. 4), was further averaged out by the random nature of the binning, particularly among bins with a greater number of gRNAs, as well as subsampling to reduce dependence on any particular sample or gRNA.

To assess the effects of a smaller gRNA library on count dispersion and statistical power, we estimated the common dispersion using edgeR following binning/subsampling on our dataset. The dispersion was estimated with a varying number of tumour samples to evaluate the degree to which the dispersion that would apply to a large dataset could be estimated from only a small number of samples; such a situation would occur, for example, when a pilot experiment is conducted to inform power calculations for a later full-sized experiment. For small tumours, our simulations showed that a reduction in BCV to approximately 0.5 could be achieved with 10,000 gRNAs, a further reduction to about 0.3 with 2000 gRNAs, which are of a similar BCV magnitude to samples of HCT116/54C GeCKO cells, and below 0.1 with very small library sizes of <200 gRNAs (Fig. 5C). When these library sizes were considered in terms of the estimated number clones/gRNA, BCV values of 0.15–0.3 were achievable with ~150–700 clones/gRNA in small tumours. For large tumours, a reduction in gRNA library size to 2000 or less reduced the BCV below 1.5, similar to in small tumours with a full library, while a very small library size of 62 gRNA was required to achieve a similar BCV to in cells of approximately 0.5 (Fig. 5D). On a per clone basis, dispersion as estimated by BCV was also higher in large tumour samples compared to small tumour samples, which is consistent with random clonal dynamics being cumulative throughout the tumour growth period. BCV estimates were less precise when determined from a small number of pilot samples, as indicated by the greater range, particularly at the lower end, when estimated from n = 3 compared to using a larger number of samples (n = 10 or n = 20), coupled with a tendency for the median to be lower at n = 3 (Fig. 5C, D).

To estimate the consequence of varying count dispersion on downstream statistical analysis, we conducted power and sample size calculations using these common dispersion estimates (Fig. 5E; Fig. S3; Table S5). When dispersion was estimated from n = 3 tumours, the estimates of the sample size required to achieve 0.8 statistical power had a greater range and lower median, with the minimum across simulations at n = 3 being approximately 1.5- to 2-fold less than that of the most conservative estimates (maximum across simulations) with n = 10/n = 20 tumours (except where the estimated sample sizes required were themselves close to n ≤ 3). This suggests that when the sample size is estimated from a smaller pilot dataset, a conservative sample size inflating factor of 2 is needed to ensure the correct level of statistical power in a full dataset. When using the more precise n = 10 dispersion estimates for small tumours, our simulations showed that at a log₂ effect size of 2, a median estimate of 35 tumours (range 34–37 across simulated datasets) is required to achieve power of 0.8 with 100,000 gRNAs, progressively reducing to 3 tumours (range: 3–4) with 2000 gRNAs. We estimate that there would be 170 clones/gRNA when using a library of 2000 gRNAs. For large tumours (dispersion estimates from n = 20 tumours) and 4000–100,000 gRNAs, >50 tumours per group were required to achieve power of 0.8 (Fig. 5E; Fig. S3A; Table S5) at effect size of 2, falling to 20–50 tumours at 500–2000 gRNAs, and <20 tumours at library sizes of <500 gRNAs. Although our power calculations indicate very small libraries are needed to achieve good power with feasible group sizes (<20–30) at this effect size, these can also be interpreted in the context of the number of clones present: we estimate ~100–400 clones/gRNA will survive in libraries of <500 gRNAs in large tumours. These estimates highlight how the reductions in count dispersion, achieved through the use of a smaller gRNA library, substantially reduce the sample sizes needed to achieve good statistical power.