We set out to develop a computational strategy to correct for the effect of off-targeting on gRNA depletion without the need to filter out unspecific gRNAs. To validate our strategy, we focused our analysis on loss-of-function screens from the 19Q4 release of the Project Achilles Avana dataset from the Cancer Dependency Map initiative2,12, performed across 26 distinct cellular lineages (Fig. 1a). This dataset represents a useful proof-of-principal scenario as previous work described the presence of unspecific gRNAs in Avana8,15 and showed that these unspecific guides confound the analysis of essentiality and contribute to false-positives hits in negative selection screens8. As important, the Avana genome-wide library targets gold-standard sets of curated essential and non-essential genes14. These gene sets are commonly used to evaluate the performance of CRISPR tools because gRNAs that target them can be considered true-positives and true-negatives in the context of viability screens, since gRNAs that disrupt essential genes are expected to drop out from the population of infected cells over time, while the abundance of gRNAs that disrupt non-essential genes is expected to remain unchanged.
a Number of screened cell lines in each of the 26 lineages represented in the dataset. b Left, number and percentage of gRNAs in the Avana library that have 0, 1, or more than 1 perfect targets (H0) in the human genome (hg38 assembly) or that have specificities (s) lower or equal/higher than 0.16. Right, number of perfect target sites for gRNAs with H0 > 1. The gRNA in the Avana library with the highest number of perfect targets sites is highlighted. c–d Summary of the number (c) and percentage (d) of genes targeted by increasing number of promiscuous gRNAs (1-6+). Promiscuity is defined at three distinct thresholds (H0 > 1, s < 0.10, s < 0.16). e Specificity scores of all gRNAs (n = 108,345) with varying numbers of perfect target sites in the genome. X-axis is truncated at 31 for simplicity. Dashed line highlights the position of the 0.16 threshold. f Specificity scores of gRNAs with varying number of perfect (left; n = 2423) or single mismatched (right; n = 2374) targets in the genome in the absence of any additional Hamming neighbors. gRNAs with no perfect target in the hg38 genome were excluded from the analysis in e and f. All boxplots show minimum, maximum, median, first, and third quartiles.
To systematically evaluate the effect of specificity on gRNA depletion, we enumerated all potential off-targets—between Hamming distances of 0 and 3—for each gRNA in the Avana library using GuideScan16. GuideScan is a retrieval-tree-based algorithm that outperforms Bowtie-based tools in the identification of off-target loci16, providing an exact and direct enumeration of the potential target space of a gRNA within a user-specified number of mismatches to the guide sequence (Supplementary Note 1 and accompanying Supplementary Table 1, Supplementary Data 1, Supplementary Data 2). Our catalog of potential off-target loci for Avana using this approach surpassed the original off-target estimates reported for this library15 as well as the estimates used by Project Achilles in the DepMap data processing pipeline2 (Supplementary Data 2, Fig. 1b–d, Supplementary Note 1). To summarize the specificity of each gRNA in this library, we further computed GuideScan’s specificity score. This score aggregates Cutting Frequency Determination values (or CFD, describing the likelihood of an off-target being cut by Cas9 based on the number, position, and identity of mismatches to a 20 nucleotide (nt)-long gRNA)15 for all potential target sites enumerated by Guidescan16, so that the most specific targeting gRNAs receive a score of 1 and the most unspecific a score of 0 (Fig. 1e, f). In agreement with previous studies8,9, gRNAs with low specificities were on average more depleted from the population during viability screens, often beyond the levels observed for gRNAs targeting known essential genes (Fig. 2a, left, Supplementary Fig. 1a–d). This observation held true even for guides that had a single perfect target site in the genome (H0 = 1) but increasing numbers of off-targets with mismatches (Fig. 2b, Supplementary Fig. 1c). Of note, when we looked at gRNAs targeting known non-essential genes14—whose representation in the library should remain unchanged over the course of the screen—we found that gRNAs with specificity scores below 0.16 were significantly depleted compared to highly specific guides (specificity score = 1; Kolmogorov–Smirnov test, adjusted for multiple testing). Fold-change distributions of gRNA with a specificity score equal or above 0.16, however, were indistinguishable from those of highly specific guides suggesting that above this threshold the effect of off-target cutting on the guide’s representation in the library is statistically minimal (Supplementary Fig. 1a, Supplemental Note 2).
a Boxplots of z-scores for gRNA log2-fold changes for Avana Project Score 19Q4 screens (n = 689) across multiple specificity bins before (left) and after (right) correction. Specificity values correspond to highest value of each bin. Dashed lines indicate median depletion of specific gRNAs (GuideScan score = 1) targeting known non-essential (top) or essential (bottom) genes. One-sided Pearson correlation values between specificity and mean depletion of gRNAs in each bin are shown, along with their significance. Boxplots show minimum, maximum, median, first, and third quartiles. b As in a but plotting only gRNAs that have a single perfect target site in the genome (H0 = 1; H0, hamming distance of 0). One-sided Pearson correlation values between specificity and mean depletion of gRNAs in each bin are shown, along with their significance. Boxplots show minimum, maximum, median, first, and third quartiles. c Schematic representation of CSC. Briefly, the correction takes as input the gRNA sequences and corresponding depletion values. For each of the gRNAs, CSC enumerates all off-targets up to Hamming distance of 3 (H3) and computes its specificity score. CSC then uses the off-target information and depletion values to model the effect of off-targeting on gRNA abundance. Finally, it outputs both the corrected depletion values and the specificity metrics for each of the gRNAs. logFC, log fold change.
To determine the extent to which off-target mediated gRNA depletion acted as a confounder in the Achilles dataset, we calculated Bayes Factors (BF) for each gene in individual screens17. In this context, BF are an assessment of gene essentiality, with positive values indicating a gene is essential and negative values indicating a gene is non-essential. Gene Set Enrichment Analysis (GSEA) showed that genes targeted by unspecific gRNAs were significantly enriched in high BF values, particularly as the number of unspecific gRNAs per gene increased or as the specificity of the gRNAs that target each gene decreased (see Supplementary Fig. 1e–g for an example cell line). This suggests that, in agreement with previous reports8,9, off-targeting may contribute to false-positive dependencies even when multiple independent gRNAs per gene are present in a library. Taken together, these data also validate the Avana dataset from Project Achilles as a suitable model to test our CRISPR Specificity Correction (CSC) algorithm for its ability to correct the confounding effect of off-targeting on gRNA depletion.
CSC takes as inputs the sequence and depletion values of all gRNAs in a screen (Fig. 2c). As a first step, it uses the sequence information to retrieve for each guide the number of potential target sites it has at zero (H0), one (H1), two (H2), or three (H3) Hamming distance to the gRNA sequence as well as their GuideScan specificity score. Although cleavage of sites with a Levenshtein distance of 1 to the gRNA—often referred to as ‘bulged’ sites—has also been reported18,19, these types of mismatches rarely preserve Cas9 activity15 and therefore were not considered. To retrieve off-target information for each guide in a library, CSC accesses hash tables whose keys contain every Cas9 gRNA that can target the human (hg38) or mouse (mm10) genomes. In the tables, each of these keys is mapped to the gRNA’s specificity metrics which have been previously computed using GuideScan16. CSC then uses the five specificity metrics as covariates in a model to assess the contribution of off-target parameters to gRNA depletion via a multivariate adaptive regression spline (EARTH, see Methods). Piecewise linear spline functions, have been previously used to model the effect of genome copy-number alterations on gRNA depletion as they can model nonlinearities and are well suited to deal with the saturating effect of high numbers of cut sites12. The EARTH model utilized by CSC has the advantage of automating both the partitioning of the data through its hinge functions as well as the variable selection step, thus providing a flexible approach for off-target correction. In addition, the backward pass improves generalizability of the model and minimizes overfitting of the training data by pruning the least effective terms and selecting the best performing model through generalized cross validation (see Methods for more details on model selection).
Specificity-corrected depletion values for each guide are outputted along with the enumeration of off-target sites and the gRNA’s GuideScan specificity score (Fig. 2c). Importantly, specificity metrics for all inputted gRNA sequences are provided even in the absence of accompanying depletion values. Thus, CSC not only allows users to correct off-target influence in their screening data but also provides direct access to GuideScan’s specificity metrics for any Cas9 gRNA that targets the human or mouse genome. This allows users to determine the potential target space of a gRNA or set of gRNAs even prior to a screen and may further aid CRISPR library design.
To validate CSC, we applied it to all screens from the DepMap 19Q4 Achilles dataset. As predicted, CSC removed the correlation between gRNA specificity and gRNA depletion (Fig. 2a, b). Additionally, when inferring gene essentiality in each cell line of the Achilles dataset, we found that correction of off-target mediated gRNA depletion by CSC significantly increased both the recall of constitutive essential genes at 5% False Discovery Rate (FDR) of constitutive non-essential genes (Fig. 3) as well as the Area Under the Curve (AUC) of precision-recall curves (Supplementary Fig. 2a, b). Both these metrics evaluate the success of predicting true gene essentialities. This trend remained true when looking at aggregate data for the entire Achilles data set (Fig. 3a), as well as at the level of individual lineages and cell lines (Fig. 3b, c, Supplementary Fig. 3). Of note, CSC markedly outperformed the filtering strategy implemented by Project Achilles to deal with unspecific gRNAs2 (Fig. 3, Supplementary Fig. 3b).
a Boxplot showing recall values at 5% FDR for the 19Q4 Project Achilles dataset (n = 689) before correction (grey), with the Achilles filter (blue), or with the CSC correction (yellow). Each dot represents the median recall value of a lineage (n = 26). Minimum, maximum, median, first, and third quartiles are shown. p-values were calculated using a two-sided Wilcoxon test. b Median recall values at 5% FDR for each lineage. c Top, example plot for a single lineage (skin). Bottom, example precision-recall plot for a single cell line (A375 melanoma cells). Lollipop graphs are plotted by increasing values of the uncorrected pipeline. The maximum recall value at 5% FDR was used for each comparison.
With the increased recall at 5% FDR, we observed a concomitant increase in the number of genes identified as essential. In total, 12,444 genes scored as a dependency in at least one screen when CSC was implemented, compared with 5831 and 6018 genes for data not corrected for off-targeting or when unspecific gRNAs were removed with the Achilles filter, respectively. To determine if these results constituted true gene essentialities, we first looked at the occurrence of known false positives. Because the number of false positives is influenced by the total number of genes identified as hits, we varied Bayes Factor thresholds for each screen to yield the same number of hits across each pipeline as previously described20. Then, for each screen we counted the number of gold-standard non-essential genes14 that had scored as hits. We found that CSC led to a significant reduction in the number false positives compared to uncorrected data, again outperforming the filtering strategy implemented by project Achilles (Fig. 4a). At this fixed number of positive hits, CSC also led to a significant increase in both precision and recall (Supplementary Fig. 2c). We then examined the expression levels of genes inferred as essential by each of the three analysis pipelines. We found that those genes identified as essential after off-target correction by CSC tended to be well expressed in the cell line in which they scored as hits (Fig. 4b, Supplementary Fig. 4a). By contrast, for each screen, genes scoring as essential in uncorrected data but not in data corrected with CSC tended to have significantly lower expression levels in the respective cell lines. In fact, a substantial subset of these genes was below bulk RNA-seq detection threshold, suggesting that they may represent false-positive hits. Similarly, genes identified as essential when data was corrected with the Achilles filter but not when data was corrected with CSC tended to be lowly expressed in the cell lines they scored as hits or below detection threshold, suggesting these may also represent false positives (Fig. 4b, Supplementary Fig. 4a).
a Boxplot showing number of false positive hits for the 19Q4 Project Achilles dataset (n = 689) before correction (grey), with the Achilles filter (blue), or with the CSC correction (yellow). For comparison purposes, Bayes Factor thresholds were varied to return the same number of hits across all three pipelines. p-values were calculated using a two-sided Wilcoxon test. b Violin plots showing the expression levels (log2(TPM+1)) of genes in the cell lines in which they were identified as dependencies. Left graph shows genes identified both before and after correction (common), only after correction (CSC only) or only before correction (uncorrected only). Right graph shows genes identified as dependencies both in pipelines that implement CSC and Achilles (common), or only in one of the pipelines (CSC only, Achilles only). Dot represents the mean value. c Boxplots showing the LOEUF scores for genes identified as dependencies by only one of the pipelines. Exclusive dependencies are defined as genes that score in more than 15 cell lines in one pipeline and none in the other (uncorrected versus CSC: n = 832; Achilles versus CSC plot n = 835). LOEUF scores for constitutive essential and non-essential genes are shown on the left as a reference (n = 1169), and their median values are highlighted across all plots with a dashed line. d Recall values as in Fig. 2a but calculated using known constitutive essential and non-essential genes targeted by at least one unspecific gRNA (H0 > 1). Filter corresponds to a pipeline where gRNAs with more than one perfect target site are removed before the analysis. Boxplots show minimum, maximum, median, first, and third quartiles.
Next, we looked for evidence of functional essentiality for genes identified as putative dependencies. Essential genes are under evolutionary pressure to maintain sequence integrity, and mutations that inactivate their functions are expected to be depleted from natural human populations. In contrast, non-essential genes whose disruption does not affect organism health or function, are expected to tolerate the accumulation of inactivating mutations. As such, we took advantage of the Genome Aggregation Database (gnomAD)21, which catalogs high-confidence predicted loss-of-function variants and uses these to classify human genes according to the mutational constraint they are under. Specifically, the LOEUF score places genes along a spectrum of tolerance to inactivating mutations, where genes that play essential cellular roles, and therefore are under high mutational constraint, receive low LOEUF scores, while genes whose disruption has no impact on cell viability or organismal health and are therefore under low mutational constraint in the human population receive high scores21. We confirmed that these scores can reflect gene essentiality by applying them to the curated sets of essential and non-essential genes14 (Fig. 4c, left). We then retrieved LOEUF scores for genes consistently identified as essential by one analysis pipeline (defined as genes that score as hits in more than 15 distinct cell lines) but not the other. We found that genes consistently identified as essential in data corrected using CSC but that did not score as essential in any screen when using uncorrected data or data corrected using the Achilles filter tended to have low LOEUF values. On average these scores were well below those attributed to constitutive non-essential genes (Fig. 4c). This suggests that genes exclusively identified as dependencies after correcting off-targeting with CSC are under mutational constraint in the human population and therefore presumably play essential roles in human cells. In contrast, genes that scored in more than 15 distinct screens only before data correction or only after removing promiscuous gRNAs through the Achilles filter—but that did not score as hits after data was corrected for off-targeting using CSC—tended to have higher LOEUF values, often above those of constitutive non-essentials. These results suggest that their inactivation may be well tolerated in humans. Taken together, these data are consistent with the notion that genes identified after computational correction by CSC reflect true essentialities, and that CSC implementation minimizes the occurrence of false-positive hits. It also suggests that CSC outperforms the current filtering approach implemented by Project Achilles to deal with unspecific gRNAs (Fig. 4a–c).
To test CSC’s performance in the context of highly unspecific libraries, we selected only genes that are targeted by at least one gRNA containing multiple perfect target sites in the human genome (H0 > 1) with the Avana library. The resulting subset of gRNAs approximates better the low specificity of libraries designed to target non-coding regulatory elements in the genome9 (Supplementary Fig. 4b, c), while still targeting true-positive and true-negative genes that can be used to benchmark CSC. The low specificity of gRNA libraries cutting within non-coding motifs stems from the fact that predicted cis-regulatory elements are typically small, with lengths several fold below those of average exons. As a result, the traditional approach used for screens of protein-coding genes, where known unspecific gRNAs are excluded at the step of library design, is not viable for a large fraction of regulatory sequences. Indeed, filtering out gRNAs suspected of off-target activity would render many of these regulatory sequences untargetable9 (Supplementary Fig. 4d).
We used this pool of unspecific Avana gRNAs to reanalyze all screens from the DepMap 19Q4 dataset and re-calculate precision and recall curves. As expected, in the context of this highly unspecific set of gRNAs, the simple filtering of known promiscuous guides was an ineffective strategy to correct off-targeting (Fig. 4d). By contrast, we found that correction of off-target mediated gRNA depletion by CSC substantially increased recall (5% FDR) across all lineages of the dataset (Fig. 4d). This suggests CSC can help retrieve known gene dependencies even from highly promiscuous gRNA pools.
To test the generalizability of CSC, we next applied it to the genome-wide screens released by Project Score1 (Fig. 5). This dataset was generated using the Sanger genome-wide library22, which was designed to contain gRNAs with 19 nt-long complementarity to the genome and therefore follows distinct off-target rules than those established for the more common 20-nt design. Like Avana, the Sanger library contains a non-negligible fraction of gRNAs with perfect (H0 > 1) or near perfect (H1 > 0) off-targets (Fig. 5a, Supplementary Data 3). In addition, as described for Project Achilles screens, gRNAs with increasing numbers of perfect target sites were increasingly depleted from essentiality screens released by Project Score (Fig. 5b, left). This was also true for gRNAs with only one perfect target site in the genome (H0 = 1) but increasing numbers of off-targets with single-mismatches (Fig. 5b, right). This suggests that even in the context of a 19-nt gRNA design, cleavage of mismatched off-targets can cause measurable toxicity to the cells. In contrast to the Achilles dataset, we found no correlation between GuideScan’s Specificity score and gRNA depletion in Project Score screens (Fig. 5c). This is expected as this score is computed using the CFD metric which was empirically developed using 20 nt-long gRNAs. Thus, the Project Score dataset provides a unique challenge for CSC and an opportunity to test its flexibility in correcting off-target mediated gRNA depletion in CRISPR viability screens based on distinct gRNA designs. Applying CSC to Project Score screens resulted in correction models that included both H0 (number of Hamming 0 targets) and H1 (number of Hamming 1 targets) as the main covariates, in stark contrast to models generated on Avana data where GuideScan’s specificity score was consistently the most important covariate (compare example metric output files for each dataset, provided as Supplementary Data 4 and Supplementary Data 5 for an Avana and Sanger screen, respectively). This reflects the weakness of GuideScan’s specificity score in summarizing off-target cleavage by 19-nt gRNAs (Fig. 5c) and highlights the ability of CSC to discriminate the most informative covariates for its correction models. As before, we evaluated the performance of CSC in the Project Score dataset by calculating the recall of known essential genes at 5% FDR (Fig. 5d). We also calculated the AUC for precision-recall curves before and after correction. We found that CSC led to a significant increase in both metrics (recall p-value <2.2e-16, AUC p-value = 5.58e-8; two-sided, Wilcoxon test). Together, these data demonstrate the generalizability of CSC and its ability to increase the sensitivity of essentiality screens performed with distinct gRNA library designs.
a Number and percentage of gRNAs in the Sanger library that have 0, 1, or more than 1 perfect targets (H0) in the human genome (hg38 assembly) or that have 0, 1, or more than 1 targets with a single hamming mismatch (H1). b Left, boxplots of z-scores of gRNA log2FC for all Project Score screens across multiple H0 bins (i.e., increasing numbers of perfect targets). Right, boxplots of z-scores for gRNA log2FC for all Project Score screens across multiple H1 bins (i.e., increasing numbers of targets with a single hamming mismatch) for gRNAs that have a single perfect targeting the genome (H0 = 1). Dashed lines indicate median depletion of specific gRNAs (H0 = 1, H1 = 0) targeting known non-essential (top) or essential (bottom) genes. One-sided Pearson correlation values between mean depletion and number of off-targets, as well as the significance of the correlation are shown below. n = 324. c As in b but plotting z-scores of gRNA log2FC against GuideScan’s specificity bins. Specificity values correspond to highest value of each bin. Note how for 19-mer gRNAs there is no correlation between GuideScan’s score and depletion. One-sided Pearson correlation values between mean depletion and number of off-targets, as well as the significance of the correlation are shown below. n = 324. d Left, change in recall (5% FDR) for each Project Score screen after being corrected with CSC. On the right we show examples of precision-recall curves for screens showing increased recall (MCAS, second most improved) or decreased recall (OAW42, second most decreased) after correction. Dashed line highlights 0.95 precision value (which corresponds to a false discovery rate of 5%). Note how the improved recall in the MCAS screen is accompanied by an increase in the area under the curve (AUC). In contrast, the decreased recall in the OAW42 screen is accompanied by no change in the AUC. All boxplots show minimum, maximum, median, first, and third quartiles.
We note that for a small number of Avana and Sanger screens (38 and 26 screens, respectively) correction of off-targeting by CSC led to an apparent drop in recall at 5% FDR (Supplementary Fig. 3 and Fig. 5d). This effect does not appear to be caused by distortions introduced by CSC, as the mean distance between the fold-changes of known essential and non-essential genes is minimally impacted by the correction and is statistically identical between screens where CSC increases or decreases recall (Supplementary Fig. 2d; two-sided Kolmogorov–Smirnov test). In addition, while increased recall is accompanied by a significant increase in the AUC of the corresponding curve (Project Achilles p-value <2.2e-16, Project Score p-value = 9.404e-09; two-sided, Wilcoxon test), lower recall is not accompanied by statistically lower AUC values (Project Achilles p-value = 0.1159, Project Score p-value = 0.4186; two-sided, Wilcoxon test). This suggests that the drop in recall at 5% FDR we have documented in a minor fraction of screens represents only a local feature of the curve and not a general trend across all decision thresholds (See Supplementary Fig. 2b and Fig. 5d for representative curves). To further ensure the robustness of the corrections and provide users with full information of how they were generated, CSC outputs for each screen a file describing the model used along with its performance metrics. These include the RMSE of the model (i.e., the square root of the variance of the residuals), a measure of how accurately it predicts the impact of off-targets on gRNA depletion (see Supplementary Data 4 and Supplementary Data 5 for example files). Additionally, CSC allows users to define RMSE thresholds above which no corrections are performed and only off-target descriptions for each gRNA in the library are outputted.
Finally, BAGEL2 is a recently published hit-calling algorithm that includes some degree of multi-targeting correction20. BAGEL2 is efficient at gene essentiality calling but it does not segregate the multi-targeting correction and hit-calling functions preventing us from directly comparing it to CSC. However, BAGEL2 uses Bowtie to generate its off-target predictions using settings that can miss perfect and near perfect target sites (see Supplementary Note 1, Supplementary Data 1, and Methods). Additionally, it discards from its analysis gRNAs for which it identifies more than 10 perfect sites or more than 10 sites with a single mismatch, and therefore will have difficulty in handling libraries with particularly unspecific gRNAs, as those that may be used for the screening of non-coding regulatory elements. Additional advantages of CSC that may be particularly useful to users include the full description of the off-target information of every gRNA submitted as an input, as well as information of how off-target metrics were incorporated in the correction model, how well the model performed, and how the correction modified original depletion values. Finally, CSC provides users with means to correct for the confounding effect of off-targets as part of their standard hit-calling pipelines, therefore providing increased flexibility in its employment.
In summary, we present here a flexible computational correction that minimizes the confounding effect of unspecific gRNAs in CRISPR-Cas9 essentiality screens leading to improved sensitivity and reduced false-positive hits in genome-wide screens. We believe CSC will be a powerful aid to ongoing efforts to catalog genomic loci required for cellular fitness, particularly in the context of screens targeting highly repetitive genomic regions—such as non-coding regulatory elements—where the design of specific libraries and gRNA filtering approaches are not feasible9. To facilitate its incorporation into existing CRISPR analysis pipelines including those that correct for genomic amplifications (Supplementary Fig. 2e, see Methods), we make the software freely available as a Python package at https://bitbucket.org/arp2012/csc_public/src/master/.
