Sequence features associated with sgRNA efficiency in CRISPR/Cas9 knocking out

To explore sequence features that contribute to sgRNA efficiency, we computed the log odds ratio of nucleotide frequency between DNA sequences targeted by efficient and inefficient sgRNAs (Fig. 2A–C). These sequences are 40 bp in length, including the 19-bp or 20-bp spacer targets as well as their 3′ and 5′ flanking DNAs. All sequences were aligned at the PAM. The log odds ratios are correlated between the ribosomal and nonribosomal sgRNA sets (Fig. 2E) as well as between the two libraries developed independently (Fig. 2F). The correlations indicate that many sequence features are robust against the variation of sgRNA designs, species, and the spacer length. Meanwhile, some features are discordant among the three different sgRNA training sets. For example, guanine is preferred at the 5′ end of the spacer in the “ribosomal” and “nonribosomal” sets, but not in the “mESC” set. This might be ascribed to subtle differences in sgRNA design. If the 5′ end of the spacer is not “G,” Wang et al. prepended a “G” for the expression of sgRNAs from a U6 promoter, in the form of G(A/C/T)X₁₉NGG; otherwise, a form of GX₁₉NGG was taken without any prepended “G” (Wang et al. 2014). Alternatively, Koike-Yusa et al. chose the form of GX₁₉NGG regardless of nucleotide composition at the 5′ end (Koike-Yusa et al. 2014). These different designs of sgRNAs can result in the variation of the nucleotide pairings and 5′ overhang in sgRNA-DNA binding structure, which could further impact sgRNA efficiency.

Open in a separate window

Figure 2.

Preference of nucleotide sequences that impact sgRNA efficiency. (A–C) Logos showing the sequence preference of the three sgRNA sets defined in Figure 1. The height of the nucleotides represents the log odds ratio of nucleotide frequency between efficient and inefficient sgRNAs. (D) A logo showing the selected features that reproducibly impact sgRNA efficiency in the three sgRNA sets. The height of the nucleotides represents the coefficients computed from the Elastic-Net. (E,F) Scatter plots showing the correlation of sequence preference for sgRNAs targeting ribosomal versus nonribosomal genes in Wang data (E) and sgRNAs in Wang data versus Koike-Yusa data (F). Each dot represents a nucleotide in a 40-bp region centered by the spacer. The sequence preference is measured as the log₂ odds ratio of nucleotide frequency between efficient and inefficient sgRNAs.

To determine the reproducible features, we first selected nucleotides that satisfy the following criteria across the three training sets: (1) Signs of the odds ratios are concordant; and (2) magnitudes of the odds ratios are above a threshold in all three sgRNA sets, where the threshold was computed from a statistical significance analysis (see Methods). To identify the dominating nucleotide features, we next applied Elastic-Net, a penalized linear regression model with proven utility in feature selection (Zou and Hastie 2005). In our analysis, the predictor variable of Elastic-Net was a binary vector representing the presence or absence of the nucleotides selected based on the odds ratio, and the response variable was 1 or −1 corresponding to the efficient and inefficient sgRNAs, respectively. The union of the three sgRNA training sets was used for Elastic-Net feature selection and parameter estimation.

With this two-stage feature selection approach, we identified 28 sequence features, as shown in Figure 2D. Most features are located within the spacer region. Our result confirmed several features reported previously. First, guanines are strongly preferred at the −1 and the −2 positions proximal to the PAM sequence, which are associated with the sequence preference in Cas9 loading (Wang et al. 2014). Second, thymines are disfavored at the four positions closest to the PAM, consistent with the fact that multiple uracils in the spacer cause low sgRNA expression (Wu et al. 2014). Third, in line with recent findings, the nucleotides downstream from the PAM contribute to sgRNA efficiency, whereas the sequences upstream of the spacer have no significant effect (Doench et al. 2014). We also identified novel features that reproducibly impact the sgRNA efficiency. For example, cytosine is preferred at the −3 position, which is the DNA cleavage site by the CRISPR/Cas9 complex (Cong et al. 2013). This feature might contribute to the efficiency of cleavage or the introduction of mutation upon DNA repair. Moreover, adenines are preferred from position −5 to −12, and guanines are preferred at positions −14 to −17.

Experimental validation of the sequence model in predicting mutation rates and protein knockout efficiency

To experimentally test the model in Figure 2D, we first assessed the mutation rates mediated by a small set of sgRNAs targeting the AAVS1 safe harbor genomic locus. We selected 10 targets with low efficiency scores and 10 with high scores in the AAVS1 region (Supplemental Table 3) and designed CRISPR/Cas9 lentivirus based on the selected targets. Upon examining the transduced 293T cells, we found that sgRNAs with higher scores showed much higher indel mutation rates than those with lower scores (P = 1.9 × 10⁻⁵; t-test) (Fig. 3A). This confirmed that the sequence features we selected have significant effect on the rate of on-target mutation in a CRISPR/Cas9 system.

Open in a separate window

Figure 3.

Experimental validation of the sequence model in predicting sgRNA efficiency. (A) A SURVEYOR gel picture (top) and a bar chart (bottom) showing the indel rates of the sgRNAs predicted to be inefficient (low sequence score) or efficient (high sequence score). The sgRNAs were selected to target the AAVS1 locus. The experiment was conducted in 293T cells. (B) A scatter plot showing the correlation of the predicted sequence scores and the protein knockout efficiency for sgRNAs targeting AR and FOXA1 in LNCaP-abl cells. The knockout efficiency is measured as the percentage of reduction in protein level upon sgRNA infection.

We next tested if the sequence model is informative to predict knockout efficiency on the protein expression level. We designed sgRNAs against AR and FOXA1 and tested them in prostate cancer LNCaP-abl cells using Western blot to assay protein expression (Supplemental Fig. 2; Supplemental Table 4). The predicted sgRNA efficiency score showed strong positive correlation with the protein knockout efficiency (P = 9.0 × 10⁻⁶; Pearson's correlation) (Fig. 3B), where the efficient and inefficient sgRNAs can be clearly separated. Our result suggested a cutoff of zero as a reasonable threshold for sgRNA classification, in line with the configurations in model training, in which the values of 1 and −1 were assigned to efficient and inefficient sgRNAs, respectively.

Predictive power of the sequence model in negative selection screens

With the satisfactory results in the low-throughput validation mentioned above, we next assessed the predictive power of the sequence model using published high-throughput CRISPR/Cas9 knockout screening data. We designed four schemes for this in silico validation. First, randomized threefold cross validation was carried out on the categorized sgRNAs in the ribosomal set to assess the predictive power under identical experimental settings; second, the Elastic-Net model was trained on the ribosomal set and tested on the nonribosomal set, of which the result was expected to reflect the within-library performance under different configurations of sgRNA sampling; third, to evaluate the inter-library and inter-species performance, the model was trained on the union of ribosomal and nonribosomal sets in human and tested on the mESC set; last, the model was trained on the union of all three sets and tested on an independent library of sgRNAs (Shalem et al. 2014). The sgRNA list used in the fourth validation was compiled by Doench et al. (2014) to assess their model in predicting sgRNA “activity,” a measure analogous to the “efficiency” defined in this paper. We achieved reasonable predictive power in all four validations, in which the Area Under Curve (AUC) scores were above 0.7 in Receiver Operating Characteristic (ROC) (Fig. 4A,B). Our model also outperformed the model proposed by Doench et al. (2014) in predicting sgRNA activity (Fig. 4B). With a cutoff threshold of zero, 50%–60% of the inefficient sgRNAs are predictable, at the cost of 10%–20% of efficient sgRNAs misclassified.

Open in a separate window

Figure 4.

Predicting sgRNA efficiency from sequence context in CRISPR/Cas9 knockout screens. (A) ROC curves showing the predictive power of the proposed model. (Red) Threefold cross-validation on the sgRNAs targeting ribosomal genes in Wang data; (blue) trained on ribosomal genes, and tested on nonribosomal genes in Wang data; (green) trained on Wang data, and tested on Koike-Yusa data. The black error bars on the red curve represent standard deviations computed from 10 iterations of random sampling in cross-validation. (B) ROC curves comparing the performance of the proposed model and the Doench et al. (2014) model in predicting sgRNA efficiency in Shalem data. (C) Scatter plot showing the correlation between the predicted sequence score and the relative sgRNA abundance for ABL1 and BCR in KBM-7 cells. The P-values were computed based on the Pearson correlation test. (D) Box plot showing the distributions of correlations between sequence scores and relative sgRNA abundances for essential and nonessential genes in KBM-7. The distribution of random background was computed by permuting the sequence scores within each gene in the data set. (E) Distributions of relative sgRNA abundances in KBM-7 cells, where the sgRNAs were categorized based on the predicted efficiency and the essentiality of their targeted genes.

Since the sequence model was trained based on the sgRNAs targeting common essential genes, we further tested if the model is applicable to the sgRNAs targeting cell-type-specific essential genes. We examined the correlations between sequence scores and the relative sgRNA abundances for ABL1 and BCR, two genes coding for an oncogenic fusion protein in KBM-7 cells (Andersson et al. 1995). For both genes, significant correlations were observed between predicted sgRNA efficiency and the changes of sgRNA abundance in the screen (Fig. 4C). To evaluate whether sequence-specific sgRNA efficiency also influences sgRNAs targeting nonessential genes, we compared distributions of the correlations for essential and nonessential genes identified in KBM-7 (Fig. 4D) or HL-60 (Supplemental Fig. 3). As expected, sequence scores and relative sgRNA abundances are in strong negative correlation for essential genes. Meanwhile, there are weak but statistically significant negative correlations for nonessential genes, which might be attributed to (1) the false negatives in identifying essential genes; (2) decreased cell growth due to intrinsic toxicity of DNA cleavage, similar to the toxicity reported on Cre recombinase (Silver and Livingston 2001); and/or (3) potential off-target effects of sgRNAs. Notably, no correlation was observed for those control sgRNAs lacking genomic targets (Supplemental Fig. 4), consistent with the requirement of chromatin binding for an sgRNA to exert its cell functions.

We next checked if a proper selection of sgRNAs based on the prediction result could refine the sgRNA pool. As shown in Figure 4E, the signals of essential and nonessential genes are not well separable with the predicted inefficient sgRNAs. In contrast, the efficient sgRNAs produced highly informative signals that distinguish essential genes from others. This indicates that the removal of inefficient sgRNAs can reduce the number of required sgRNAs without compromising the sensitivity of the CRISPR/Cas9 knockout screens.

Generality of the sequence model in positive selection screens

Our model was trained based on the sgRNAs targeting essential genes that were negatively selected in the experiments. We next asked if the model is applicable to positive selection screens. To answer this, we curated four lists of genes involved in the resistance to drug or toxin treatments, including the BRAF inhibitor Vemurafenib (PLX) in melanoma A375 cells (Shalem et al. 2014), alpha-toxin and 6-thioguanine (6TG) in mouse embryonic stem cells (Koike-Yusa et al. 2014), and anthrax toxin in HeLa cells (Supplemental Table 5; Zhou et al. 2014). The sgRNAs targeting these genes were categorized based on the reproducibility of positive selection, i.e., the increase of the sgRNA abundance in multiple biological replicates upon treatments. As shown in Figure 5A–E, the predicted efficient sgRNAs expressed higher reproducibility of positive selection in all the experiments, indicating our model's applicability to positive selection screens. Next, we assessed the predictive power of our model and the model by Doench et al. (2014) on the curated sgRNA sets, in which the sgRNAs “selected in all replicates” were defined as positives and the sgRNAs “not selected in any replicates” were defined as negatives. Our model achieved an AUC score of 0.711 and significantly outperformed the Doench et al. (2014) model (Fig. 5F). Notably, 42.3% of the efficient sgRNAs contain a cytosine at the cleavage site, in comparison with 25.9% for inefficient sgRNAs, suggesting this novel feature also contributes to sgRNA efficiency in positive selections.

Open in a separate window

Figure 5.

Assessment of the sequence models in predicting sgRNA efficiency in positive selection experiments. (A–E) Bar charts showing the capability of selection and the experimental reproducibility for predicted efficient and inefficient sgRNAs. The tested sgRNAs target the genes known to be involved in the resistance to different drug treatment or external stimulus. (F) ROC curves comparing the performance of the proposed model and the Doench et al. (2014) model in predicting sgRNA efficiency in positive selection experiments. In the evaluation, the positive test set consists of the sgRNAs selected in all replicates in B–E; and the negative test set consists of those not selected in B–E.

Our prediction is therefore consistent with the phenotype in both positive and negative selections, where the abundance of efficient sgRNAs is expected to change in opposite directions. This line of evidence, together with our validation results on mutation rate and protein knockout efficiency, confirmed that the selected sequence features impact cell phenotype mainly through the loss-of-function of the target genes, but not through their intrinsic toxicity on cell growth.

Selecting phenotype-associated genes

For the CRISPR/Cas9 knockout data sets, we ran the beta-0.1 version of RIGER (http://www.broadinstitute.org/cancer/software/rnai/riger/) to call the negatively selected genes (Luo et al. 2008). In the configurations of RIGER, the measures of relative sgRNA abundance were used to rank the sgRNAs, and a Kolmogorov-Smirnov test was chosen to determine the significant genes. The threshold of P-values was set to be 0.001. The positively selected genes in Figure 5 were collected based on experimental validations in the original studies (Koike-Yusa et al. 2014; Shalem et al. 2014; Zhou et al. 2014), and a cutoff threshold of two folds in relative sgRNA abundance were used to determine efficient and inefficient sgRNAs. The phenotype-associated genes in the CRISPRi/a data set were retrieved from the original paper (Gilbert et al. 2013).

Thresholds for identifying efficient and inefficient sgRNAs

As shown in Figure 1, we used different thresholds for identifying efficient and inefficient sgRNAs in Wang data and in Koike-Yusa data. This is because the sgRNAs targeting essential genes in Wang data and Koike-Yusa data have different distributions of log ratios due to the variance of cell growth rate, Cas9 activity, and the depth of selection in the experiments. In Wang data, the sgRNAs with a log ratio smaller than −1.0 in both HL-60 and KBM-7 were identified to be efficient, and the thresholds for inefficient sgRNAs are in the fifth percentile of the log ratios of control sgRNAs in HL-60 and KBM-7. The thresholds in Koike-Yusa data were determined based on the bimodal distribution of signals (Supplemental Fig. 1). Notably, when the thresholds in Koike-Yusa data were set to be identical to those used in Wang data, we observed little variation of nucleotide frequency (r = 0.96). Therefore, our sequence model is insensitive to the threshold chosen.

Extracting sequence features

For Wang, Koike-Yusa, and Shalem data, sequences of spacer targets and their flanking regions were extracted from hg19 or mm9 genome assembly. Since the length of sgRNA spacers varies from 18 bps to 22 bps in the Zhou et al. (2014) platform, we mapped the spacer sequences to hg19 assembly to find the genomic loci of the targets and extracted the 40-bp sequence (spacer + 3′ + 5′) based on the aligned PAM loci. The sequences in CRISPRi/a data were retrieved from the original paper (Gilbert et al. 2014; Konermann et al. 2015).

Determining reproducible features in CRISPR/Cas9 knockout screens

To measure the sequence preference, we computed the log₂ odds ratio of nucleotide frequency between efficient and inefficient sgRNAs, for the “ribosomal,” “nonribosomal,” and “mESC” sets. Let $s_i^1$, $s_i^2$, and $s_i^3$ represent the log₂ odds ratios for the ith nucleotide computed from the 3 sgRNA sets, we define a score s_i as

where sgn(s) is the sign of s; and $\left|s\right|$ is the magnitude of s.

To determine the background distribution of s_i, we randomly permuted the sgRNAs in the efficient and inefficient categories within each set and computed s_i based on the permuted data. This results in an empirical null distribution of s_i. A threshold T was then set to be at the 95% confidence interval of the null distribution, and the nucleotides with $\left|s_i\right|>T$ were selected as reproducible features.

The Elastic-Net model

Suppose $\mathit{X}=[X_1,X_2,\dots ,X_N{]}^{\mathit{T}}$ is the set of encoded sequence vectors, and $\mathit{Y}=[y_1,y_2,\dots ,y_N{]}^{\mathit{T}}$ is the set of outputs representing the efficiency of sgRNAs, where N is the number of sgRNA samples for training. Let M be the length of the input vectors, the Elastic-Net model computes the parameters $\beta =[{\beta }_1,{\beta }_2,\dots ,{\beta }_M{]}^{\mathit{T}}$ that minimize an objective function E

where α and λ are parameters estimated using cross validation, ${||\beta ||}^1=\sum _i\left|{\beta }_i\right|$, and $||\beta |{|}^2=\sum _i{\beta }_i^2$. We used the glmnet R package to implement the Elastic-Net.

When training the Elastic-Nets for sequence feature selection, we found the optimal α to be 1.0 for CRISPR/Cas9 knockouts and 0.0 for CRISPRi/a (Supplemental Figs. 5, 6). We note that the Elastic-Net is equivalent to the Least Absolute Shrinkage and Selection Operator (LASSO) when α = 1.0 (Tibshirani 1996) and is equivalent to the Ridge regression model when α = 0.0 (Tikhonov and Arsenin 1977). In practice, the LASSO model is preferred to select a small number of features from a large amount of candidates, whereas the Ridge regression is more relevant in situations when many features cumulatively contribute to the response. Therefore, our result on the selection of optimal α is explainable since the sgRNA efficiency depends on fewer sequence features in Cas9 knockout in comparison with that in CRISPRi/a.

Application of previous model

To apply the Doench et al. (2014) model, we downloaded the python script at http://www.broadinstitute.org/rnai/public/analysis-tools/sgrna-design. For fair comparison, the tested sgRNA targets and our regression model (Fig. 2D) were truncated into a 30-bp sequence harboring the spacer target, following the requirement in the Doench et al. (2014) model.

Cell lines and cell culture

The LNCaP-abl (abl) cell line was provided by Zoran Culig (Innsbruck Medical University, Austria). The abl cells were cultured in the RPMI 1640 phenol red-free medium supplemented with 10% charcoal/dextran- treated fetal bovine serum, 2 mM glutamine, 100 µg/mL penicillin, and 100 units/mL streptomycin for the experiments. The 293T cells obtained from the American Type Culture Collection were maintained in DMEM media supplemented with 10% fetal bovine serum, glutamine, and penicillin-streptomycin.

Plasmid construction, lentivirus production, and transduction into human cells

Twenty sgRNAs targeting the AAVS1 locus, nine sgRNAs targeting the AR gene, and six sgRNAs targeting the FOXA1 gene were designed and selected according to our prediction of sgRNA efficient scores. The sgRNA oligos were cloned into lentiCRISPRv2 plasmid (Addgene) as previously described (Shalem et al. 2014). Each plasmid containing inserted sgRNA sequence was verified using Sanger sequencing. To make lentivirus, the constructed lentiCRISPRv2 plasmids were cotransfected into 293T cells with packaging plasmids pMD2.G and psPAX2 (Addgene) using X-tremeGENE HP DNA Transfection Reagent (Roche) in 12-well plates according to the manufacturer's instructions. These 35 types of packaged lentivirus were then transduced into 293T and LNCaP-abl cells using 24-well plates, followed by puromycin selection for 3 d. The surviving cells were maintained for another 3 d before isolation of the DNA and proteins.

SURVEYOR nuclease assay and Western blot

The DNA fragments of the AASV1 locus were amplified from extracted genomic DNA by PCR using the Q5 high-fidelity DNA polymerases (NEB). SURVEYOR nuclease assays were subsequently performed according to the manufacturer's instructions (Trans genomic). The indel mutation rates were calculated as previously described (Ran et al. 2013). Western Blot for AR and FOXA1 knockout in LNCaP-abl cells by different sgRNAs was carried out using whole cell lysates as described previously (Xu et al. 2012), The antibody used for immunoprecipitation was anti-AR (N-20, Santa Cruz Biotechnology) and anti-FOXA1 (ab23738, Abcam), or anti-GAPDH (sc-25778, Santa Cruz Biotechnology). The images were analyzed using ImageJ (Schneider et al. 2012).

We thank Drs. Carl Novina, Tim Wang, David Sabatini, Eric Lander, Ophir Shalem, and Neville Sanjana for providing the data sets and helpful discussions. The project was supported by the National Institute for Health Research (NIHR) grants U01 CA180980 (to X.S.L.), R01 GM113242-01 (to J.S.L.), National Science Foundation (NSF) grant DMS-1120368 (to J.S.L.), and the Claudia Adams Barr Award in Innovative Basic Cancer Research from the Dana-Farber Cancer Institute.

Author contributions: H.X. and X.S.L. designed the study. C.A.M., J.S.L., and M.B. provided conceptual advice. H.X., C.H.C., and D.W. developed the algorithm and performed the analyses. T.X. conducted experimental validations. W.L., T.X., L.C., F.Z., and M.B. helped with technical clarifications and result interpretations. H.X., T.X., and X.S.L. wrote the manuscript. All the authors participated in discussions and manuscript revision.