Predicting CRISPR/Cas9 Off-Target Effects with Mismatches and Insertions Using Stacked BiGRU Models

We first assess whether our Crispr-SGRU can accurately predict off-target activity. We compare it with four existing hybrid CNN and BiLSTM methods, namely, CRISPR-Net [21], CrisprDNT [26], CRISPR-IP [23], and CRISPR-M [24], on eight public imbalanced datasets: CHANGE-seq, I1, K562, HEK293t, BE3, II5, I2, and II6. Prior to this, we briefly comment on some comparisons of these methods. First, CrisprDNT used hybrid one dimensional CNN and BiLSTM, whereas others applied two-dimensional CNN and BiLSTM. Second, besides CNN and BiLSTM, CRISPR-IP and CrisprDNT used attention mechanism and transformer, respectively. Third, CRISPR-M operated by incorporating both sgRNA–DNA sequence pairs and epigenetic data, whereas others considered only sgRNA–DNA sequence pair information. Thus, we compared these methods by considering sgRNA–DNA sequence pairs along to ensure fairness of the comparisons. Fivefold cross-validation is applied in the training and testing phases to achieve the averaged results. The quantitative results for each of the predictors are represented as heatmaps in Figure 1. We observe that models trained and tested on balanced datasets achieve superior performance than on imbalanced datasets. Our model achieves comparable performance to others in terms of recall (0.693) and outperforms them in terms of precision, F1 score, and MCC, with mean values of 0.607, 0.641, and 0.644, respectively. Crispr-SGRU shows excellent performance on the HEK293t dataset (IR = 247), with precision of 0.703, recall of 0.758, F1 score of 0.729, and MCC of 0.729. We also observe that Crispr-SGRU exhibits remarkable performance on the imbalanced II5 dataset (IR = 1773.6) and II6 dataset (IR = 6846.6) in terms of F1 score and MCC.

As shown in Figure 2 and Figure 3, our model outperforms the comparison models on both the AUROC and PRAUC evaluations. It shows the best performance with a PRAUC of 0.521 and 0.442 on imbalanced datasets II5 and II6, with a significant improvement of 11% and 4.9% compared with the second-best method. On dataset I1, our model outperforms other methods in terms of AUROC and PRAUC, with values of 0.986 and 0.757, respectively. Similarly, we observe that Crispr-SGRU remarkably suppresses other methods on datasets HEK293t, I2, and II6 in terms of these two evaluation indicators. Furthermore, our model outperforms CrisprDNT on datasets I1 and I2 on the measurement of PRAUC, with values increasing by 0.4% and 5.8% on average, respectively. For more details, see Supplementary Table S1. We also assessed the performance of our Crispr-SGRU on balanced datasets. For the sake of clarity, we compared it with two better methods, CRISPR-Net and CRISPR-IP as in Dhanjal et al., and II1 datasets under five-fold cross-validation, and the averages of the individual performances are summarized in Supplementary Table S2. Our Crispr-SGRU clearly outperformed the compared methods on balanced datasets. Collectively, these observations indicate that Crispr-SGRU outperforms competing methods.

To assess the robustness of our model for off-target prediction, we compared it with the above four methods on datasets K562, HEK293t, BE3, and II5 under a leave-one-sgRNA-out procedure. Methods assessed by each sgRNA separately can reduce the risk of biasing the results towards sgRNAs with dominant evaluation calculation. For a fair comparison, we adopted the same dataset under a leave-one-sgRNA-out cross-validation procedure for model training and testing to avoid the impact of sgRNA selection. Specifically, for each dataset, we randomly selected the samples of an sgRNA as an independent testing set and the remaining sgRNAs and the corresponding sequence pairs served as a training set in each iteration. To assess the performance on independent datasets, for each sgRNA, all methods were trained on the rest of this dataset, excluding the dataset containing this sgRNA. This approach can guarantee that the training set and testing set have no sgRNA overlap.shows the details of the selected sgRNAs and the number of positive and negative samples in each dataset. Supplementary Table S3

As Figure 4 and Figure 5 depict, Crispr-SGRU shows excellent performance on the K562 dataset, with an average F1 score of 0.828 and MCC of 0.826. It achieves comparable performance to CrisprDNT, with averaged value of 0.880 and 0.910 in terms of all evaluation metrics. Our model significantly surpasses CrisprDNT on the HEK293t dataset based on precision, recall, F1 score, MCC, and PRAUC measurements, with improvements of 0.9%, 23%, 4.1%, 6.9%, and 20%, respectively. In addition, Crispr-SGRU has the second-highest AUROC (0.945), about 4% lower than CRISPR-Net on this dataset. Crispr-SGRU achieves the best averaged AUROC (0.999), PRAUC (0.895), and MCC (0.695) on the BE3 dataset, with a 4% improvement in PRAUC and a 1.4% improvement in MCC over the second-best method. Moreover, it achieves the highest average F1 score (0.721) and MCC (0.740) on the II5 dataset. This corresponds to a 5.2% improvement in F1 score and a 5.4% enhancement in MCC over the second-best method. Together, these results suggest that our method surpasses others in terms of robustness.

We evaluated the effectiveness of different network architectures for Crispr-SGRU on off-target prediction. To this end, we carried out a comparative analysis of two variants derived from our final model by modifying its components. First, we verified the efficacy of the Inception by proposing a variant architecture (Crispr without Inception), getting rid of the Inception layer from the full model. We directly used the output of encoded sequences as inputs for the stacked BiGRU instead. Second, we confirmed the efficacy of the stacked BiGRU by constructing a variant architecture (Crispr without BiGRU) by removing the stacked BiGRU. We used the flattened output of the Inception layer as the input of the dense layer. Then, we compared these variants with Crispr-SGRU under five-fold cross-validation on the K562 and HEK293t datasets.

The experimental results are shown in Table 1. We observed that these two variants performed significantly worse than the final model, suggesting that both Inception and stacked BiGRU are essential. Furthermore, we discovered that removing Inception layers can obviously degrade the performance metrics. This proves the effectiveness of multi-scale convolutional operations in detecting the abstract sequence features. Notably, we observed that eliminating the stacked BiGRU module resulted in a significant decrease in all performance metrics across these two datasets, leading to 19.7% and 4.8% decreases on average on the K562 and HEK293t datasets, respectively. Thus, the stacked BiGRU is indispensable in our model for its ability to capture the rich sequence dependencies and contributes most to off-target prediction.

Next, we compared the computational efficiency of Crispr-SGRU and four deep learning-based methods—CRISPR-Net, CrisprDNT, CRISPR-IP, and CRISPR-M—on the HEK293t dataset. The batch size and epoch were set to 256 and 30, respectively. The running time results can be found in. Overall, our method demonstrates a certain level of advantage in efficiency among all methods. We also observe that CrisprDNT requires the most time to execute. This is expected, because CrisprDNT utilizes a transformer module, which inherently possesses high computational complexity. In contrast, CRISPR-IP introduces a single-layer attention mechanism, which runs seven times faster than CrisprDNT. These results illustrate that though deep networks can better capture the characteristics of data, the computational cost of the model increases exponentially as the number of hyperparameters increases. Thus, there is a trade-off between performance and computational cost. Supplementary Table S4

We were also interested in the interpretability of our model. We applied Deep SHAP [37] to estimate the position-dependent nucleotide contribution to our model on six mismatch-only datasets: CHANGE-seq, K562, HEK293t, BE3, II5 and II6. For each dataset, we used all positive samples and randomly selected an equal number of negative samples for analysis. To make the contribution of nucleotides comparable among various datasets, we rescaled the SHAP values by standardization. As depicted in Figure 6, nucleotides at positions 10–20 exhibit a strong positive contribution to off-target prediction on all datasets. This result aligns with previous studies, which corroborate that nucleotides at 8–12 positions adjacent to the PAM determine Cas9 specificity [38,39,40]. In addition, we observe the nucleotide at position 18 has a positive contribution on all datasets, which is in line with previous findings that nucleotides at positions 2–3 distal to the PAM are less tolerant of mismatches [41,42]. Our model consistently pays close attention to positions 1–3. This observation is in accordance with previous findings, suggesting that mismatches at positions 1–3 are less likely to be tolerated [42]. We also find that positions 1–5 have a positive contribution on CHANGE-seq, K562, HEK293t, and II5 datasets. These observations are in agreement with Vora et al. [40], which demonstrates that the PAM-distal-most positions provide valuable information.

shows the performance of Crispr-SGRU and two student models: the student model with KD (Student with KD) and student model without KD (Student without KD) on K562 dataset under fivefold cross-validation. Compared to Student without KD, Student with KD boosts the performance in terms of all evaluation metrics, with AUROC, PRAUC, F1 score, and MCC increased by 0.2%, 3.1%, 1.1%, and 4.6%, respectively. This observation illustrates that knowledge is successfully transferred from the teacher model to the student model. Thereby, Student with KD can be applied to interpret the decision basis of Crispr-SGRU. Supplementary Table S5

To improve the reliability and transparency of decision-making in the model and reveal the role of the stacked BiGRU in the student model, we visualize the model weight of Student with KD of each base-pairing of an sgRNA–DNA sequence pair and analyze the important sequence fragments obtained from the model. We apply Deep SHAP to visualize the model weight maps for 16 base pairs for each position in sgRNA–DNA sequence pairs for Student with KD. Figure 7 shows the importance of the sequence context for each position in sequence pairs for Student with KD on K562 dataset. We observe base-pairing at positions 14–20 achieve higher Deep SHAP values, which means these positions have significant influence on prediction among all sites of sequence pairs. Specifically, C-A:18 (representing the mismatch site C-A at position 18) has the highest SHAP value, which means it has the most significant positive impact on off-target prediction. We also find that the mismatch sites A−T:19, C−T:19, C−G:19, and T−A:20 have great effect on predictive outcomes. This phenomenon is expected, since CRISPR/Cas9 specificity depends on the sgRNA seed sequence of the PAM-proximal region [43]. In addition, G−C:17, C−G:18, and T−G:18 are influential, with positive contributions to prediction. This is consistent with previous observations that the cut site-adjacent sequence is an important determining factor for Cas9-induced mutations [44,45,46]. Furthermore, we observe that some mismatch sites of the PAM-distal region, such as G−C:1, C−A:4, and C−A:5, have great influence on off-target prediction, which coincides with a previous finding that mismatches at positions 1–3 are less likely to be tolerated [42]. Together, these results indicate that our model can learn critical sequence information in sequence pairs for off-target prediction.

We used eight public datasets to evaluate the performance of our model. These datasets can be divided into two categories: (i) mismatch-only datasets (e.g., Dhanjal et al. [47], II1 [48], CHANGE-seq [49], K562 [50], HEK293t [50], BE3 [51], II5 [52], and II6 [53]), and (ii) mismatch-and-indel datasets (e.g., I1 [21] and I2 [21]). Each entry in these datasets is composed of an sgRNA sequence, a target DNA sequence, and a binary label, with 0 representing the negative sample and 1 representing the positive sample. All datasets contain only off-target sites with NGG PAMs. Dhanjal et al.’s dataset consists of 9214 positive samples and 9917 negative samples. The II1 dataset consists of 2273 positive samples and 2580 negative samples. These two datasets are balanced, with imbalance ratios of 1.10 and 1.14, respectively. The CHANGE-seq dataset contains 67,476 positive samples and 2,806,151 negative samples. The HEK293t and K562 datasets contain 18 and 12 sgRNAs, respectively. Among all 30 sgRNAs, 656 samples have been identified as off-targets. The HEK293t dataset consists of 536 positive samples and 132,378 negative samples. The K562 dataset contains 120 positive samples and 20,199 negative samples. The BE3 dataset contains 9 sgRNAs accounting for 79 validated off-targets [51]. The II5 dataset has 95,775 negative samples, while the number of positive samples is only 54. The I1 dataset consists of 10 sgRNAs, with a total of 7371 active off-targets, 340 with 1 bp indel and up to 3 bp mismatches and 7031 with up to 6 bp mismatches. Dataset I2 comprises 213,943 sequence pairs from 6 sgRNAs. This dataset is comprised of 60 validated off-targets, of which 13 contain indels. The II6 dataset contains 52 positive samples and 383,407 negative samples. More details are shown in Table 2.

We use the one-hot encoding scheme proposed by Lin et al. [21] to encode sgRNA–DNA sequence pairs, where the DNA or RNA bulge (indel) is denoted using underscores. Each sequence pair is encoded as a seven-row binary matrix, with a five-bit channel (e.g., A, G, C, T, and indel) to encode the sgRNA and DNA nucleotides and a two-bit direction channel to identify the mismatch and indel directions. Specifically, sgRNA and DNA sequences are encoded by one-hot encoding, subsequently being performed by OR operation to get the results of the first four channels (e.g., AA (1, 0, 0, 0), TT (0, 1, 0, 0), GG (0, 0, 1, 0), CC (0, 0, 0, 1)). The fifth channel marks the DNA or RNA bulge, where indel occurs is 1, otherwise it is 0. The sixth and seventh channels are applied to identify the insertion or deletion and the mismatch type. To adapt the indel information in this encoding scheme, we add 1 nt using hyphens beside the target sequence and the PAM sequence to form a 24 nt base sequence or base pair sequence, which is encoded by (0, 0, 0, 0). As such, each sequence pair is encoded by a binary matrix with a size of (length of the sequence pair) ×7. We direct readers to the original publication [21] for more details.

Our Crispr-SGRU model (Figure 8a) has two parts: an Inception-based CNN (Figure 8b) and a stacked BiGRU (Figure 8c). The CNN module is regarded as a feature extractor that learns the local and global features of sequence pairs independently. The stacked BiGRU is designed to recognize sequential patterns. Each sequence pair is encoded into a binary matrix, subsequently being fed into a CNN. The CNN module consists of four two-dimensional convolutional layers with varying filter sizes: 1 × 1, 1 × 2, 1 × 3, and 1 × 5, respectively. The units of these layers are set to 10. Utilizing varying perceptual fields can combine feature information at various scales during parallel processing [55]. This enables the internal layer to choose the suitable filter size for learning the required information of the sequences. Additionally, the parallel pooling operations in Inception can reduce the dimensions of the feature map. The outputs of Inception are input into the stacked BiGRU, which consists of three BiGRU layers of 30, 20, and 10 units, to learn the contextual information of the sequence. This is inspired by previous studies that stacked multiple RNN layers and extracted more generalized sequential features [32,56,57]. The outputs of the stacked BiGRU are flattened and fed into two dense layers with 128 and 64 units. The rectified linear unit (ReLU) [58] is applied as the activation function. We use a 0.35 dropout rate for model regularization to combat overfitting. The output layer consists of two neurons using a sigmoid activation function that quantifies propensity for on- and off-target sites. We use an Adam optimizer [59] with learning rate of 0.0001 to train our model to minimize the loss function.

We regarded the off-target prediction as a binary classification task. This is inspired by a previous study that off-target assessment applying classification outperforms regression [50]. Several class-balanced loss functions have been proposed to solve the problem of training from imbalanced data by introducing a weighting factor inversely proportional to the effective number of samples [36], such as focal loss [60], Tversky loss [61], asymmetric loss [62], dice loss [63], and hinge loss [64]. The formulas of these loss functions are summarized in Supplementary Table S6. Inspired by previous work [26], we focused on the loss function to address the data imbalance problem while training our model. We compared the performance of Crispr-SGRU trained using the abovementioned loss functions and observed that using dice loss achieves superior performance to others (Supplementary Table S7). Thereby, we used dice loss in our model to reduce the impact of favoring the negative samples.

KD [65] is applied to transfer the dark knowledge learned by the larger teacher model to a small student model. Some studies have also used this strategy in training to improve model interpretability. For instance, Lu et al. proposed a KD-insight drug–target affinity prediction method, which takes atomic-level pockets as input and learns drug–target atomic-level interaction information, thus giving the model better interpretability [66]. Fan et al. developed a deep learning-based method called ETFC to predict 21 categories of therapeutic peptides [67]. They constructed a student model with KD to investigate the contribution of each AA in the ETFC to each peptide sequence and enhanced the interpretability of their model. We note that sequence pairs are encoded by the encoding scheme proposed by Lin et al. [21] in our model, so the base-pairing information in each sequence pair is inevitably disrupted by an OR operator, which means our model lacks decision-making transparency. Inspired by the above studies, we constructed the KD-insight off-target prediction model, which takes a 16-channel encoded matrix as input and investigates the contribution of each base pair to each sgRNA–DNA sequence pair, thus giving the model better interpretability.

Figure 9a shows the teacher–student workflow for KD. We defined the pretrained Crispr-SGRU as the teacher model, while the student model is based on it with fewer parameters. The details and the encoding scheme and student model are exhibited in Figure 9b,c. We designed a 24 × 16 encoding scheme to encode sequence pairs, where 24 represents the length of the sequence pairs. For an sgRNA–DNA sequence pair, the base pair is represented by 16 channels, i.e., AA-, AT-, AC-, AG-, TA-, TT-, TC-, TG-, GA-, GT-, GC-, GG- CA-, CT-, CG-, CC-, and CG-. Each position in each sequence pair is related to a vector of length 16 with a single non-zero element corresponding to the base pair at that position. For example, AA and CG are represented as vectors [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] and [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], respectively. This encoding scheme can preserve the information of both sequence and mismatch types, thus avoiding loss of important information during model training and enabling accurate interpretation of the base-pair details.

In the KD process, we trained the student model with real sample labels as hard labels, which typically convey information about the predicted classes. The predictive outcomes of the teacher model were regarded as soft labels that could offer more comprehensive insights and salient knowledge derived from the probability distribution across all classes during training the student model. We utilized the KD model to extract more detailed information from the teacher–student model and obtained refined soft labels using the student model. We employed distillation loss (DL) to quantify the loss incurred during model training. We used student loss (SL) to denote the loss incurred during this process. Thus, the loss function of the student model considered the distribution of labels from the teacher model, which is defined as follows:whereis a hyperparameter used to balance the losses ofand. We setin the following experiment. (1) L K D = ∗ + ( − ) ∗ α S L α D L 1 α S L D L α = 0.2

Crispr-SGRU was performed using Python 3.8 and Keras library 2.4.3 with a Tensorflow (2.5.0) backend. All experiments were carried out on a desktop computer with an Intel (R) Xeon (R) Silver 4210R CPU (Intel, Senta Clara, CA, USA) @ 2.40 GHz, Ubuntu 18.04.6 LTS, and one NVIDIA GeForce RTX 3080 Ti with 12 GB of memory (MSI, New Taipei, China). For each dataset, we adopted the standard data split: 85% and 15% for model training and testing, respectively. We trained the models with a batch size of 256 and epoch of 30. Additionally, we applied the bootstrapping sampling algorithm [68] to ensure that each batch had equivalent positive samples and negative samples in the training data, thus preventing gradient update instability during model training.

To assess model accuracy, we used four commonly used evaluation metrics—precision, recall, F1 score [69], and Matthews correlation coefficient (MCC) [70]—which are defined as follows:(2)Precison=TPTP+FP(3)Recall=TPTP+FN(4)F1score=2×Precision×RecallPrecision+Recall(5)MCC=TP×TN−FN×FPTP+FN×TN+FP×TP+FP×TN+FN where TP is the number of correctly predicted off-targets, FN is the number of off-targets predicted to be non-off-targets, TN is the number of correctly predicted non-off-targets, and FP is the number of non-off-targets predicted to be off-targets. The F1 score offers a balanced measurement by considering both precision and recall, thus providing insights into both accuracy and completeness of the model. Additionally, we constructed receiver-operating characteristic (ROC) [71] and precision–recall (PR) curves [72] to visualize model ability and calculated the corresponding area under the ROC curve (AUROC) and area under the PR curve (PRAUC) to quantify the model’s performance. AUROC can better evaluate the prediction performance when the labels are highly imbalanced, which is the case for off-target prediction.