A machine learning approach for predicting CRISPR-Cas9 cleavage efficiencies and patterns underlying its mechanism of action

Abstract

The adaptation of the CRISPR-Cas9 system as a genome editing technique has generated much excitement in recent years owing to its ability to manipulate targeted genes and genomic regions that are complementary to a programmed single guide RNA (sgRNA). However, the efficacy of a specific sgRNA is not uniquely defined by exact sequence homology to the target site, thus unintended off-targets might additionally be cleaved. Current methods for sgRNA design are mainly concerned with predicting off-targets for a given sgRNA using basic sequence features and employ elementary rules for ranking possible sgRNAs. Here, we introduce CRISTA (CRISPR Target Assessment), a novel algorithm within the machine learning framework that determines the propensity of a genomic site to be cleaved by a given sgRNA. We show that the predictions made with CRISTA are more accurate than other available methodologies. We further demonstrate that the occurrence of bulges is not a rare phenomenon and should be accounted for in the prediction process. Beyond predicting cleavage efficiencies, the learning process provides inferences regarding patterns that underlie the mechanism of action of the CRISPR-Cas9 system. We discover that attributes that describe the spatial structure and rigidity of the entire genomic site as well as those surrounding the PAM region are a major component of the prediction capabilities.

Fig 1. Schematic flow of the cross-validation procedures.

The main components of the learning pipeline for the leave-one-sgRNA-out and leave-study-out cross-validation procedures are presented. ¹ This step was applied to the leave-one-sgRNA-out procedure only. ² In each iteration, the samples of a single sgRNA (in the case of the leave-one-sgRNA-out procedure) or all samples from a single study (in the case of leave-study-out) were excluded from the training data and used as a test set. The algorithm was trained on the rest of the data. ³ Each set of cleaved samples (targets that correspond to a single sgRNA) was oversampled using bootstrapping, thus introducing a subset twice the size of the original one, and an equal-sized set of uncleaved samples was randomly chosen. ⁴ For each original set of cleaved samples in the test set (targets that correspond to a single sgRNA), an equal-sized set of uncleaved samples was randomly chosen.

Fig 2. Comparison of four prediction algorithms on the assembled dataset.

(a-d) Pearson correlation coefficient computed over all the samples in the dataset. The horizontal axis represents the scaled observed values published in the experimental studies, and the vertical axis represents the scores predicted by: (a) CRISTA applied using cross-validation (r² = 0.65), (b) CCTop (r² = 0.23), (c) OptCD (r² = 0.13), (d) CFD score (r² = 0.52). (e) Receiver Operator Characteristics curves computed over all the samples in the test dataset: CRISTA (AUC = 0.96), CCTop (AUC = 0.85), OptCD (AUC = 0.85), CFD score (AUC = 0.91). Positives and negatives represent cleaved and uncleaved sites, respectively. True (and false) positives rate is computed as the true-positives (false-positive) number divided by the number of positives (negatives). (f) Precision-Recall curves computed over all the samples in the dataset: CRISTA (AUC = 0.96), CCTop (AUC = 0.87), OptCD (AUC = 0.88), CFD score (AUC = 0.93). Precision is computed as the true-positive number divided by the sum of true-positives and false-positives. Recall is computed as the true-positives number divided by the positives number. (g) Pearson correlation coefficient computed for each sgRNA: CRISTA (averaged r² = 0.80, sd = 0.13), CCTop (averaged r² = 0.46, sd = 0.22), OptCD (averaged r² = 0.32, sd = 0.28), CFD score (averaged r² = 0.65, sd = 0.28). (h) Receiver Operator Characteristics curves computed for each sgRNA: CRISTA (averaged AUC = 0.99, sd = 0.02), CCTop (averaged AUC = 0.86, sd = 0.13), OptCD (averaged AUC = 0.9, sd = 0.12), CFD score (averaged AUC = 0.9, sd = 0.11). (i) Precision-Recall curves computed for each sgRNA: CRISTA (averaged AUC = 0.99, sd = 0.02), CCTop (averaged AUC = 0.92, sd = 0.09), OptCD (averaged AUC = 0.93, sd = 0.07), CFD score (averaged AUC = 0.94, sd = 0.06). Mean values are marked with horizontal lines. The whiskers reach 1.5 times past the first and third quartiles.

Fig 3. Accuracy across different studies in a leave-study-out cross-validation.

(a) Observed cleavage intensities versus predicted intensities. The top and bottom rows represent the nuclear targets of the ‘unique guides’ and ‘common guides’, respectively. Pearson r² values are shown; "overall" represents the correlation calculated by taking all points, and "mean" is the average correlation calculated for each sgRNA individually. Different colors represent nuclear targets of different sgRNAs. (b, c) ROC and PRC curves. The ‘unique guides’ and ‘common guides’ of each study are represented by different curves. AUC values are denoted in the legend. Each column corresponds to a single experimental platform.

Fig 4. Features importance.

Clustering of top-ranked features and their relative importance. The nodes sizes represent the feature importance as calculated by CRISTA. Edges transparencies represent correlation such that strongly correlated features are connected by darker edges. Yellow and blue edges represent positively and negatively correlated features respectively. Abbreviations: YY- mismatches of type pyrimidine-pyrimidine; RR–mismatches of type purine-purine; MGW–minor groove width; ‘#’ represents counts (for further explanations of the features, see S3 Table). The graph was produced with Cytoscape [55] using the pairwise correlation for every pair of features and their importance scores.

Fig 5. Forward selection results.

The top plot represents the ROC-AUC, PRC-AUC, r², and root mean square error (RMSE) following the addition of every feature from left to right. The bars represent feature importance, i.e., the contribution of every feature to the prediction accuracy as computed by the Random Forest algorithm. The RMSE is divided by two for visualization.