Deep learning models simultaneously trained on multiple datasets improve base-editing activity prediction

Cell culture and doxycycline-inducible ABE/CBE stable cell lines generation

Wildtype HEK293T cells from ATCC (catalog number CRL-3216) were cultured in DMEM media supplemented with 10% FBS and 1% penicillin-streptomycin (D10 medium) at 37 °C with 5% CO2. Routine mycoplasma testing using PCR (cat no. PM008) yielded negative results the cells used for this study. TRE-ABE7.10 and TRE-CBE (BE4) stable cell lines were generated using PiggyBac transposon systems. HEK293T cells were transfected with pPB-ABE7.10-hygromycin (modified from the Addgene#102919), pPB-CBE-hygromycin (modified from the Addgene#100806) vectors, and pCMV-hybase at a 9:1 ratio. HEK293T-ABE7.10 and HEK293T-BE4 cells were cultured in selection medium containing 50 μg/mL hygromycin.

Lentivirus transduction with the 12 K library

HEK293T-ABE7.10 and -BE4 cells were cultured in D10 medium supplemented with 50 μg/mL hygromycin throughout the experiment. The 12 K library was generated previously¹⁷. For the lentiviral 12 K CRISPR gRNA-target pair library transduction, the following steps were carried out: 1) Day −1: 2.5 × 10^6 cells per 10 cm dish were seeded in 12 dishes. One dish per group was allocated for cell number determination pre-transduction, one for drug-resistance (puromycin) control, and the remaining 10 dishes for lentivirus library transduction; 2) Day 0: Cell count was performed to estimate cell number per dish, determining the volume of lentivirus for transduction with a MOI of 0.3. Low MOI (0.3) ensured most infected cells received only one copy of lentivirus. Infected cells were cultured in D10 medium at 37 °C; 3) Day 1: Transduced cells were split into three dishes per original dish; 4) Day 2: Sub-group 1 (10 dishes) was harvested and labeled as Day 2 post-12K library transduction, with cells pooled for genomic DNA extraction. Sub-group 2 (10 dishes) received fresh D10 medium with 50 μg/mL hygromycin + 1 μg/mL puromycin (Dox-free). Sub-group 3 (10 dishes) received D10 medium with 50 μg/mL hygromycin + 1 μg/mL puromycin + 1 μg/mL doxycycline (Dox induction of base editor expression). For WT HEK293T cells (Group 3) screening, hygromycin but not puromycin was used in the culture medium; 5) Transduced cells were split every 2–3 days when confluence reached 90%. 6) The transduced cells were harvested for further genomic DNA extraction and targeted deep sequencing eight days after growing in the Dox selection medium.

Targeted deep sequencing

Genomic DNA from transduced cells was extracted using the phenol-chloroform method and subsequently purified. PCR amplification of the surrogate target sites was carried out using the TRAP-NGS-F1 and TRAP-NGS-R1 primers (sequences listed in Supplementary Table 2). Each PCR reaction contained 5 μg of genomic DNA as a template, yielding approximately 7.6 × 10^5 copies of the gRNA-target pair construct, providing about 63-fold coverage of the 12 K library. To achieve comprehensive coverage, 32 parallel PCR reactions were conducted, resulting in approximately 2016-fold coverage of each construct. The PCR products were then purified using a 1.5% gel and pooled in equal amounts for subsequent deep amplicon sequencing. The amplicon libraries were subjected to deep sequencing on the MGISEQ-2000 (MGI of BGI in China) platform. All the samples were subjected to paired-end 150 bp deep-sequencing on the MGISEQ-2000 platform.

Raw data processing

The raw data were sequenced using the MGISEQ-2000 platform (G400), generating PE150 paired-end reads. FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, version: 0.11.3) and fastp (https://github.com/OpenGene/fastp, version: 0.19.6)²⁴ were used for data quality control and filtering with the default parameters. The paired-end data were assembled using flash (http://www.cbcb.umd.edu/software/flash, version: 1.2.11)²⁵. BWA-MEM (version: 0.7.17)²⁶ was used to map the assembled data to the designed oligo sequences to preliminarily distinguish the data of each synthetic gRNA-target pair surrogate site. The oligo sequences have been reported in the previous CRISPRon study and are available in Supplementary Data 1¹⁷.

Data filtering

The pysam module of Python-3.8 was used to split the aligned data according to the site number of the chip, and the reads of different sites were obtained. Then, we used three steps of strictly controlling parameters to filter the data of each site. Firstly, according to the structure of the chip, g + gRNA (20 bp) + scaffold (82 bp) + GTTT should remain unchanged at the beginning and end of each site. Then, to remove the chip synthesis errors, the pseudo-editing sequences found in WT HEK293T cells transduced with the 12 K SURRO-seq library were removed from both ABE and CBE groups. The above three filtering steps were completed with Julia-1.5.3 language.

Calculation of base editing efficiency

To improve the suitability for machine/deep learning, the SURRO-seq dataset was further processed and base editing efficiency and outcome editing efficiency were calculated separately. Base editing efficiency is relative to the number of sites edited, and for each site, the frequency of mutating (number of mutations/total number of reads) was calculated, including the frequency of occurrence of indels, although at a lower percentage. To mitigate the impact of enrichment and depletion, we kept this subset of the dataset with a total number of reads greater than or equal to 100, resulting in 11,484 gRNAs for ABE and 11,406 gRNAs for CBE. The outcomes for all edited results were calculated by masking indel loci as wild type. The frequency (number of mutations per outcome/total number of reads) and proportion (number of mutations per outcome/number of reads in which mutations occurred) were then calculated separately.

Sequence context analysis for base editing efficiency

To explore the impact of proximal sites on the editing efficiency of base editors, we constructed sequence motifs for both ABE and CBE using our SURRO-seq dataset. We employed a logistic regression model to predict gRNA editing efficiency and generated sequence logos based on the obtained weights. Initially, we selected sequences containing the target nucleotides (“A” for ABE and “C” for CBE) along with their flanking nucleotides ( ± 3nt). These sequences were then encoded using a one-hot encoding scheme for efficiency prediction. The gRNA editing efficiency was coded as 0 or 1 to facilitate binary classification modeling. If the editing efficiency exceeded the mean value at certain position, the target value was 1; otherwise, it was 0. We subsequently utilized a logistic regression model for efficiency prediction, employing 5-fold cross-validation by keeping highly similar target sequences (≤ 4 mismatches for a 30nt target sequence) in a single fold. Then model performance was assessed on the 6th partition with equal size reserved for independent testing using Pearson correlation coefficient between actual and predicted efficiency values. Additionally, we extracted weights from the logistic regression model to visualize the effect of sequence context on base editor efficiency. Given the high efficiency observed from positions 4 to 8 within the spacer, we focused our motif analysis on target sequences with target nucleotides located within these positions.

Data collection and pre-processing for machine learning

In addition to our SURRO-seq dataset, we collected datasets for both ABE and CBE: Song Train/Test from Song et al., ⁷ Arbab dataset from Arbab et al. ⁶ (although there is no information about which part of the data was used for their training and testing) and Marquart Test from Marquart et al. ⁸. We also collected two ABE datasets (Kissling ABE7.10 and ABE8e) from Kissling et al. ¹². Note that no detailed information was provided regarding which data were used for training and testing in the original study. In total, three base editors are involved in this study: ABE7.10 (including SURRO-seq, Song, Arbab, Marquart, Kissling ABE7.10 dataset), ABE8e (Kissling ABE8e dataset) and BE4 (SURRO-seq, Song, Arbab, Marquart dataset). In addition, we not only collected the data evaluated in HEK293T cells but also the Arbab dataset from other cell-lines, including mESC (Arbab mESC Test) and U2OS (Arbab U2OS Test) to explore the model performance across cell lines. Further details about the datasets are given in Supplementary Table 3.

The datasets were processed as follows: 1) Removal of gRNAs without “NGG” PAM at the DNA target sites; 2) Removal of gRNAs for which the corresponding target sites are supported by less than 100 total reads; 3) Removal of gRNAs with over 100,000 total reads, considering that their total reads are higher than other gRNAs; 4) Removal of gRNAs not containing the nucleotide A (ABE) or C (CBE) within the editing window (positions 3 to 10); 5) Exclusive retainment of outcomes with intended target nucleotide transitions (A > G for ABE, C > T for CBE) within the editing window, considering that only few reads supported transitions to unintended nucleotides, such as A > Y for ABE and C > R for CBE and the transitions occur mainly in the editing window (Fig. 1C, D). The outcomes from one target DNA sequence were further grouped according to the positions of the transitioned nucleotides in the editing window and the frequencies within each group were summed; 6) Generation of outcomes at target sites with the potential to be edited but without any reads from the experimental support, by setting their outcome frequency to 0; 7) Removal of gRNAs for which the corresponding target site is supported by fewer than 100 edited reads, as in previous studies^6,7. As we observed in the Song datasets that the sum of the number of reads for the outcomes edited by a specific gRNA can be higher than the total number of reads reported for this same gRNA, we recalculated the outcome frequencies accordingly. See Supplementary Table 4 for more details about dataset processing and Supplementary Table 5 for information on the general composition of these datasets.

Each target sequence is 30nt long, encompassing 4nt upstream, 20nt target DNA protospacer (corresponding to the spacer), 3nt PAM and 3nt downstream from the PAM. To develop our model, we combined our SURRO-seq dataset, Song Train, Arbab dataset, Kissling ABE7.10 and Kissling ABE8e datasets (train and/or test set), while Song Test, Marquart Test and Arbab mESC and U2OS Test were only used for model evaluation (test set). The combined dataset consists of 17,941 unique 30nt target DNA sequences for ABE and 19,010 unique sequences for CBE. The identical 30nt sequences between our SURRO-seq dataset and Song Train were modest (20 gRNAs from ABE and 22 gRNAs from CBE) and were retained in both dataset sources. While the gRNAs from Kissling datasets are generated from the same library, 932 gRNAs exist in both datasets after above filtering steps and all the gRNAs were retained. See Supplementary Fig. 13 for details.

Generation of target DNA and outcome features

We processed the 30nt target DNA sequences and their corresponding edited outcome sequences by the following steps: 1) the 30nt target DNA sequences were one-hot encoded; 2) outcome sequences in the editing window were one-hot encoded/labeled (edited or not) based on the positions of the target nucleotide transitions in the editing window. Therefore, each outcome sequence was represented by one vector with length equal to that of the editing window. The value in the vector was set to 1 if the nucleotide was edited, and 0 otherwise. 3) The spacer RNA-target DNA binding energy ΔG_B was calculated using the CRISPRoff pipeline 1.1.2¹⁹ supplied with RNAfold 2.5.1²⁷. 4) CRISPR/Cas9 indel frequencies were predicted by CRISPRon¹⁷ (see details in Supplementary Note 2). 5) The origin of each training set was encoded in a “labeling” vector of size 5 for ABE (SURRO-seq dataset, Song Train, Arbab dataset, Kissling ABE7.10 dataset, Kissling ABE8e dataset) and 3 for CBE (SURRO-seq dataset, Song Train, Arbab dataset), as detailed in Supplementary Table 1.

Generation of dataset partitions

The dataset was split into six partitions with approximately equal sizes (within a count of one gRNA) for model development as previously done in CRISPRon¹⁷, resulting in five partitions for cross validation and the 6th partition for test. The procedure took the sequence similarity between 30nt target DNA sequences into account by assigning highly similar ones to the same partition while remaining dissimilar ones were randomly assigned into different partitions. In addition, we considered that part of the targets were included in the training of CRISPRon¹⁷, which is applied for Cas9-gRNA efficiency prediction, and that none of the data points in the 6th partition (test set) should belong to the CRISPRon training set to prevent overestimation (Supplementary Note 2). This was done as follows (see detailed workflow in Supplementary Fig. 14): 1) The pairwise Hamming distance between all 30nt target DNA sequences was calculated using the scikit-learn pairwise_distances function²⁸, and extracting highly similar sequences (≤4 mismatches for a 30nt target sequence), as for CRISPRon¹⁷; 2) The 30nt target DNA sequences employed in CRISPRon model development across six validation sets were assigned to the same partitions as in CRISPRon¹⁷; 3) 30nt target sequences highly similar to those assigned in step 2 were put in the same partition as their similar ones; 4) Remaining similar DNA targets were grouped and assigned to one partition randomly; 5) The remaining dissimilar targets from different sources were randomly distributed into different partitions. The test dataset from Song et al. ⁷ provided by the authors, was assigned to the 6^th partition for subsequent benchmark analysis, while all the 30nt sequences from Song Train were divided into the five partitions used for training. Sequences from Song Train which were assigned to the test partition of CRISPRon were not allocated to any partitions in step 2 but were randomly distributed into five partitions in step 5. Additionally, the DNA targets involving same variant loci in Kissling dataset were also grouped as highly similar ones as suggested by the authors¹², which were assigned to different partitions in step 3 and 4. After these assigning steps, 30nt target DNA sequences and their potential edited outcomes from diverse sources were distributed across six partitions without highly similar sequences among partitions and keeping all the data not used for CRISPRon training in the 6th partition, to maintain the test set fully independent for further benchmark analysis (Supplementary Table 6).

Test sets for the evaluation and comparison of models

Test sets from various dataset sources were collected: SURRO-seq Test (the 6th partition from SURRO-seq dataset), Song Test, Arbab Test (the 6th partition from Arbab dataset), Kissling ABE7.10 Test (the 6th partition from Kissling ABE7.10 dataset, ABE only), Kissling ABE8e Test (the 6th partition from Kissling ABE8e dataset, ABE only), Marquart Test and Arbab mESC and U2OS Test (see Generation of dataset partitions). To create fully independent test sets, each test set was processed as follows: 1) Pairwise Hamming distance between test sets and training sets used in any method was calculated on 20nt gRNAs by the cdist function in SciPy²⁹; 2) The gRNAs with Hamming distance ≤ 3 (number of mismatches) on 20nt gRNAs from any training set were removed as was done in CRISPRon¹⁷; 3) The 30nt target DNA sequences used to train the final version of CRISPRon were removed and the different test sets were reduced to the following number of gRNAs: SURRO-seq Test: 868 (ABE), 1297 (CBE); Song Test: 430 (ABE), 474 (CBE); Arbab Test: 1,096 (ABE), 1180 (CBE); Marquart Test: 1482 (ABE), 1,489 (CBE); Kissling ABE7.10 Test: 144 (ABE), Kissling ABE8e Test: 152 (ABE); Arbab mESC Test: 2554 (ABE), 1481 (CBE); Arbab U2OS Test: 361 (ABE), 289 (CBE). In addition, we extracted the outcomes with a frequency greater than 5% in each test set for benchmark analysis considering these outcomes are more critical for practical application, which are called “SURRO-seq Test*”, “Song Test*”, “Arbab Test*”, “Marquart Test*”, “Kissling ABE7.10 Test*” and “Kissling ABE8e Test*” respectively. To evaluate the performance on a dataset including the vast number of low efficient gRNAs appearing the unprocessed datasets we made a test containing the firstly discarded gRNAs (edited reads <100) resulting in “SURRO-seq Test with low efficient gRNAs”: 2326 (ABE), 2165 (CBE). Correspondingly, by removing low efficient gRNAs ( < 100 edited reads) from the Song Test, we obtained “Song Test without low efficient gRNAs”: 383 (ABE), 338 (CBE).

Development of deep learning models designed for base editors

Deep learning models were developed for ABE and CBE using Keras/Tensorflow³⁰ 2.10.0 and Python 3.10. The models were designed to simultaneously predict gRNA editing efficiency and outcome frequency from a 30nt target DNA sequence and an edited outcome. The model requires four inputs: the 30nt target DNA sequence, the 8nt labeling of editable positions, and two molecular features: ΔG_B, a powerful feature for Cas9-gRNA efficiency prediction¹⁷, and the predicted Cas9-gRNA efficiency, which is positively correlated with base editor efficiency⁷ (Supplementary Fig. 1). The one-hot encoded 30nt target DNA sequences are fed into three convolutional neural networks with filter sizes of 3, 5, and 7. The outputs from these three filters are flattened and then further passed through one fully connected layer. The one-hot encoded outcome vector representing the edits of target nucleotides in the editing window is concatenated with the molecular features and outputs from the three convolutions. Subsequently, this combined input is fed into two sequential fully connected layers before computing a 2D output: the gRNA editing efficiency and the outcome frequency (Supplementary Note 1 and Supplementary Fig. 2).

To enhance the generalization performance on diverse datasets, our SURRO-seq dataset was combined with the published datasets Song Train, Arbab dataset, Kissling ABE7.10 (ABE only) and Kissling ABE8e datasets (ABE only) for model development. One fused model was proposed for training on a combined dataset while keeping the original values of gRNA editing efficiency and outcome frequency. Due to the variations between diverse datasets (i.e., the same gRNA may exhibit different efficiencies across datasets), we encoded the dataset of origin as a feature (Supplementary Table 1) and introduced this feature alongside others in the model. The complete model architecture is outlined in Fig. 3A with details in Supplementary Fig. 7.

Models were trained using a 5-fold cross-validation on five partitioned datasets with the 6th partition held out for testing. All the models were trained and evaluated using the mean squared error (MSE). During training, the model underwent 2000 epochs in which weights were updated using batches of 16 samples, with a learning rate of 0.0001 using ADAM optimization. The training was stopped when the MSE computed on the validation set did not improve for 100 consecutive epochs and the best model by MSE on the validation data was saved. To prevent overfitting during training, a dropout rate of 0.1 was applied in each layer. The hyper-parameters were tuned to optimize the deep learning models. Inspired by the network from DeepBE¹⁰, the number of filters in convolutional neural networks and number of nodes in fully connected layers were increased (see Supplementary Data 2 for detailed hyper-parameter optimization). During hyper-parameter optimization, each training in the 5-fold cross-validation was repeated four times with different random seeds and the best model of the four was chosen from each validation. The final MSE for model selection was calculated as the average MSE from the five best models, weighed by number of outcomes in each validation. In the final model, each training in the 5-fold cross-validation was repeated 10 times with different random seeds for robustness. The final output was the average of the output from the best five models chosen for each fold.

Benchmark analysis of base editor prediction methods

CRISPRon-ABE and CRISPRon-CBE were benchmarked against existing base editor methods, including DeepABE/CBE⁷, BE-HIVE⁶, BE-DICT⁸, BE_Endo¹¹, BEDICT2.0¹² (ABE only) and DeepBE¹⁰ (ABE8e only) for the prediction of the gRNA editing efficiency and outcome frequency on independent test sets. For BE-DICT⁸ and BEDICT2.0¹², the outcome frequency with target nucleotide transitions in the editing window (from positions 3 to 10 in 20nt gRNA spacer) was predicted as suggested by the authors. Considering all the test sets are genome-integrated target sequences, we only compared the performance with BE_Endo¹¹ trained using sequence information. All the above models were downloaded from their GitHub repositories.

Jointly evaluation of gRNA editing efficiency and outcome frequency

To evaluate the performance of the base editor efficiency predictions, we calculated an extended K-dimensional Pearson correlation coefficient R_K (K = 2)¹⁸ and an extended K-dimensional Spearman rank correlation coefficient, ρ_K (K = 2), between ground truth and prediction results. These metrics allow us to combinedly evaluate the performance in predicting gRNA editing efficiency and outcome frequency. The formulas of R_K and ρ_K are explained as follows (see Supplementary Note 4 for details):

We consider two tables X (e.g. ground truth) and Y (e.g. prediction) with N rows (number of outcomes) and K columns (here K = 2, gRNA editing efficiency and outcome frequency). The covariance between X and Y was defined in ref.¹⁸ as:

where \({\bar{X}}_{k}=\frac{1}{N}\mathop{\sum }_{n=1}^{N}{X}_{nk}\), \({\bar{Y}}_{k}=\frac{1}{N}\mathop{\sum }_{n=1}^{N}{Y}_{nk}\) and \({w}_{k}=\frac{1}{K}=\frac{1}{2}\).

The correlation coefficient is dessfined as

For R₂ calculation, X and Y are the values of gRNA editing efficiency and outcome frequency. For ρ₂ calculation, X and Y are the ranks of gRNA editing efficiency and outcome frequency. It can be shown that ρ_K takes a similar form as Spearman rank correlation when K = 1 (see Supplementary Note 4).

Statistics and reproducibility

No statistical method was used to predetermine sample size. The detailed data processing steps were described in “Methods” section. The datasets for model development are provided with the software package for reproducibility. All statistical analyses were performed using Python (version 3.10.6) SciPy library (version 1.10.1).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.