How DNA packaging and surrounding sequence affect prime editing success

The cell line HEK293T was obtained from the American Type Culture Collection (ATCC) and was maintained in DMEM supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin at 37°C with 5% CO₂. HEK293T cells were transfected in 24-well plates with Lipofectamine 2000 (Invitrogen) following the manufacturer’s instructions. The amount of DNA used for lipofection was 1 μg per well. Transfection efficiency was always higher than 90% as determined by fluorescent microscopy or flow cytometry, following the delivery of a control GFP expression plasmid.

The plasmids encoding PE2 (Plasmid#132775), as well as the plasmid expressing the pegRNA (Plasmid #132777), were obtained from Addgene. PegRNAs were cloned into the plasmid backbone with paired oligonucleotides (IDT, Supplementary Table S1), as previously described (Anzalone et al., 2019). Briefly, the oligonucleotides used to create the guide sequences were hybridized, phosphorylated, and cloned into the sgRNA vector using BsaI-HFv2 (NEB), in a reaction that included T4-PNK (NEB) and T4 DNA Ligase (NEB).

DNA amplicons for NGS were generated by PCR using KAPA HiFi HotStart (Roche), according to the manufacturer’s instructions, using primers with overhangs compatible with Nextera XT indexing (IDT,). Following validation of the quality of PCR products by gel electrophoresis, the PCR products were isolated using AMPure XP PCR purification beads (Beckman Coulter). Indexed amplicons were then generated using a Nextera XT DNA Library Prep Kit (Illumina), quantitated, and pooled. Libraries were sequenced with a MiSeq Nano flow cell for 251 cycles from each end of the fragment using a MiSeq Reagent Kit v2 (500 cycles). FASTQ files were created and demultiplexed using bcl2fastq v2.17.1.14 Conversion Software (Illumina). Deep sequencing was performed by the Roy J. Carver Biotechnology Center at the University of Illinois at Urbana-Champaign. Supplementary Table S2

The in-house experiments of prime editing for 25 target sites () were conducted in three biological replicates on different days by the same person. Supplementary Tables S1, S2

The editing efficiency data for PE2 at endogenously edited sites in Anzalone et al. (2019) were obtained from the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/↗) under the accession code PRJNA565979. In-house experimental sequencing data for 25 target sites can be found in the SRA under the accession code PRJNA949853. Raw sequences from Anzalone et al. and in-house experiments were aligned to the human genome (GRCh38) using the Bowtie package (version 2.4.1), with a maximum fragment length of 500 base pairs (bp) for paired-end sequences (Langmead and Salzberg, 2012). Poorly aligned sequences were filtered using the Samtools package (version 1.7) with options -h -F 4 -q 10 (Li et al., 2009). The editing efficiency for each biological replicate was calculated by dividing the number of reads that contain the edit of interest, but were otherwise perfectly aligned, by the number of perfectly aligned reads containing either the wild-type or edited sequence. For paired-end sequencing data, half of the number of reads found in strand pairs was excluded from both the mutant and wild-type sequences to prevent double-counting of reads. The final editing percentage for each target locus was determined by averaging over the editing percentages in all biological replicates. The processed editing efficiencies in Kim et al. (2021) were obtained from their Supplementary Tables 3, 4.

Genomic locations of H3K9me3 and H3K27me3 histone modifications were determined from the aggregate of chromatin immuno-precipitation sequencing (ChIP-seq) datasets in HEK293 cells from the Encyclopedia of DNA Elements (ENCODE) consortium (ENCODE Project Consortium, 2012; Luo et al., 2020), with the accession codes having the prefix “ENC,” and the Gene Expression Omnibus (GEO) (Barrett et al., 2013), with the accession codes having the prefix “GSM.” For H3K9me3, the IP sequencing data used were GSM4301086 (Broche et al., 2021), the combination of ENCFF002AAX and ENCFF002AAZ as replicates, and the combination of GSM3452796 and GSM3452797 also as replicates (Tarjan et al., 2019). The corresponding inputs were GSM4301092, ENCFF000WXY, and GSM4445881 (Lamb et al., 2019). Raw reads were aligned to GRCh38 using Bowtie (version 2.4.1). For H3K27me3, the IP sequencing data were GSM3907592 (Lamb et al., 2019), GSM4301076 (Broche et al., 2021), GSM4586041 (Fan et al., 2020), and GSM4859391 (Fan et al., 2021); the corresponding inputs were GSM3907592, GSM4301081, GSM4586041, and GSM4859385. To determine the consensus genomic locations of H3K9me3 enrichment relative to the input, we first constructed six row vectors from the l1-normalized density of reads, partitioned into 1-kb genomic bins, with the three IP data and three input data treated as separate samples. We performed singular value decomposition (SVD) on the resulting matrix of the six row vectors (Supplementary Method S2). A similar analysis was performed to detect the consensus genomic locations of H3K27me3 enrichment relative to the input. SVD and further analysis of variance (ANOVA) were performed using the decomposition and classification of epigenomic tensors (DeCET) package (Leistico et al., 2021).

RNA-DNA hybridization energies were obtained by computing the difference in length-normalized Gibbs free energy (ΔGo) at 37°C between a paired RNA-DNA oligomer and two unpaired oligonucleotides. The ΔGo values of paired oligonucleotides were computed by adding up all the ΔGo values of the dinucleotide components of the oligonucleotides, in addition to a helix initiation term that accounts for forming the first base pair in the double helix. The ΔGo values of all paired dinucleotides and initiation terms were obtained from Sugimoto et al. (1995).

Nucleosome occupancy in the lymphoblastoid cell line GM12878 was determined from MNase-sequencing bigwig tracks from the ENCODE portal under the accession code ENCFFOOOVME (Luo et al., 2020). Average nucleosome occupancy in a particular genomic range was calculated by determining the sum of MNase-seq signals at every genomic coordinate in the protospacer region, normalized by the length of the protospacer.

PE off-targets were computed by aligning the protospacer sequence to regions in the hg38 genome that are also upstream of an NGG PAM site. A genomic locus was considered to be an off-target if its aligned sequence matched the reference protospacer sequence up to three mismatches outside the PAM site and if there were no mismatches in the GG dinucleotide of the PAM. These off-targets were determined by using the Cas-OFFInder package (Bae et al., 2014).

Linear regression with combined l1 and l2 penalties was computed in R (version 4.2.1) using the function glmnet. cv in the glmnet (version 4.1–4) library (Zou and Hastie, 2005). Given a set of N sequences, sequence Si of length L and edit percentage EPSi was modeled according to the equation

whereis the indicator variable for the presence of A, C, G, or T at theth position in theth sequence,are the regression coefficients for the indicator variables, andis the constant intercept. The vector of regression coefficients is estimated as I i j k j i β j k β 0

whereis the vector of edit percentages andandare parameters to be tuned via cross-validation. After performing 10-fold cross-validation with a default value of, the value ofwas tuned to minimize the mean cross-validated error across all 10 folds of the observed edit percentage compared to the predicted edit percentage of the withheld validation dataset. E P S λ α α = 0.5 λ

It was previously observed that having G as the last templated nucleotide at RTTL = i tended to decrease the editing efficiency compared to the shorter pegRNA design with RTTL = i−1 (Anzalone et al., 2019). In the absence of G, we noticed that the nucleotide A had a similar effect as G. We thus hypothesized that the stepwise differences in editing efficiency between two consecutive RTTLs can be predicted by the last templated nucleotides at those RTTLs. To confirm our hypothesis, we built a linear regression model represented by the equation

whereis the edit percentage at RTTL =in the range [10 nt, 20 nt] for each target locus,is the indicator variable for nucleotidebeing theth templated nucleotide, thes are the regression coefficients, andis an error term. E P i i I i N , N i a ϵ

The data provided by Kim et al. contained editing efficiencies for PBSL = 13 at only nonconsecutive RTTLs (RTTL = 10, 12, 15, and 20); thus, the OLS was applied on the sum of stepwise differences in three ranges of RTTL: 1) ∑10≤i<12ΔEPi, 2) ∑12≤i<15ΔEPi, and 3) ∑15≤i<20ΔEPi. To assess the performance of the OLS model, 8-fold or 10-fold cross-validation was performed using data of Anzalone et al. or Kim et al., respectively. The target loci were partitioned into eight (or 10) sets, and OLS was trained on the ΔEPi (or ∑iΔEPi) at the target loci in the union of seven (or nine) sets (training set), where 10≤i<20. The OLS was validated on the remaining held-out set (test set). This process was repeated eight (or 10) times permuting the test set.

We trained a DNN model to predict the edit percentages of all pegRNA designs from Kim et al. (2021) and thereby learn the salient sequence features influencing prime editor efficiency. The inputs to the DNN were 47 × 5 matrices, where the dimension 47 was chosen to accommodate the length of the “wide-target sequences” in Supplementary Table 4 of Kim et al. The first four columns of the input matrix were used for one-hot encoding of the edited strand sequence at locations from −21 to +26 relative to the nick site. To account for the fact that a given site could be targeted by different pegRNA designs of varying lengths, the fifth column of the input matrix was used to one-hot encode for the coverage of the corresponding locations by each input pegRNA.

To capture position-specific effects of sequence features relative to the nick site, we used independent filters for different positions in the target sequence, instead of using convolutional filters with shared weights. The input sequences were thus divided into all possible overlapping 8-mers, resulting in 40 8-mers for a sequence of length 47 nt. Each 8-mer was passed through 10 filters of kernel size, and the resulting output was then flattened to a vector. The flattened output was processed through a fully connected layer of 10 neurons. The output of the fully connected layer was then passed through a single neuron whose output was the prediction of the edit percentage of the pegRNA design at the corresponding target site. All layers of the DNN used a rectified linear unit (ReLu) activation function. To prevent overfitting, one dropout layer, where 20% of layer outputs were set to 0, was applied after the input layer, and another dropout layer was applied after the filter layer. The DNN was constructed and trained using the Python package Keras version 2.2.0. 8 5 ×

The full dataset from Kim et al. was divided into 80% training, 10% validation, and 10% test sets, grouping together pegRNA designs with the same target site but with different PBSLs or RTTLs into a common set. Since the edit percentage of approximately 60% of the pegRNA designs was less than 10%, we assigned training weights to each pegRNA design in the training set in order to balance the dataset and prevent the DNN from learning biased features from over-represented poorly edited sites. For this purpose, the training set was partitioned according to their observed edit percentages, with a bin size of 1%. A pegRNA design, indexed by i, was given a weight

where NEP is the number of pegRNAs in the partition corresponding to the range [intEP%; intEP+1%). The aforementioned weights were applied to the loss function during DNN training. Training was terminated when there was no improvement in the weighted mean-squared error loss between the predicted logEP and the observed logEP for the validation set after 10 consecutive epochs. The DNN model at the epoch with the lowest validation loss was considered for further analysis. The Pearson correlation was determined between the predicted and observed edit percentages in the test set. The DNN was trained using the stochastic gradient descent optimizer with Nesterov momentum of 0.9, initial learning rate of 10^–3, and decay factor of 10^–6.

To extract sequence features learned by the DNN, we used a modified version of the simulated annealing (SA) algorithm to determine the optimal sequences that maximize the output of the DNN (Tsallis and Stariolo, 1996; Finnegan et al., 2020); we then used the MaxEnt algorithm to identify salient positions and nucleotide preferences in the optimal sequences (Finnegan and Song, 2017) (Supplementary Methods S3, S4).

In order for a PE to be able to edit its target site, the spacer region of pegRNA must first access and hybridize to the complementary DNA sequence at the targeted locus. We thus hypothesize that target sites residing in heterochromatin, a condensed state of highly packaged DNA, may be blocked from access by PEs and thereby have low editing efficiency. Consistent with our reasoning, it has been previously observed that chromatin structure can interfere with the CRISPR-Cas9 endonuclease activity (Wu et al., 2014; Daer et al., 2017; Yarrington et al., 2018; Liu et al., 2019). We thus investigated the effect of heterochromatin on PE efficiency using genome editing data from Anzalone et al. (2019), Kim et al. (2021), and additional in-house validation experiments. To determine heterochromatin locations in HEK293 cells, which were used in all three studies, we integrated several publicly available H3K9me3 and H3K27me3 ChIP-seq datasets as a proxy for indicating closed chromatin regions (Methods). After obtaining consensus regions of enrichment for each of these histone modifications using singular value decomposition of the joint IP and input data matrices (Methods; Supplementary Method S2; Supplementary Figure S1), we computed the genomic distances from each target site to the nearest H3K9me3- and H3K27me3-modified regions and used these distances as features for predicting editing efficiency. When considering edit percentage as a function of the distance from the protospacer to the nearest H3K9me3-modified region, we observed that the edit percentages at sites close to H3K9me3 were approximately 0%, whereas higher edit percentages were possible farther away (Supplementary Figure S2A).

This observation motivated us to search for binary classification of target sites as being either weakly or strongly editable, choosing an efficiency threshold of 1% to binarize the data. We then trained a multivariate logistic regression model to classify all 32 endogenous target sites in Kim et al., using the distances to H3K9me3- and H3K27me3-modified regions as features (Methods; Supplementary Method S2). The probability threshold separating the binary editing categories was set to maximize the finite positive likelihood ratio of the true positive rate to the false positive rate. Performing the t-test for each regression coefficient and determining the bivariate and univariate logistic regression decision boundary showed the H3K27me3 variable to be insignificant (p=2.52×10−1;Supplementary Table S3; Supplementary Figure S2B). After removing the insignificant feature and retraining the model, the probability threshold that maximized the finite positive likelihood ratio corresponded to a distance threshold of 26 kb for the nearest H3K9me3 peak from the protospacer (Figure 1A). We observed a substantial increase in PE efficiency at the target sites located at least ∼10 kb away from H3K9me3 modification, suggesting that chromatin assumed an open conformation around this distance threshold. Similarly, it was reported that transcription of genes started to increase ∼10 kb away from H3K9me3 peaks (Barski et al., 2007). Using the 26 kb distance threshold, the accuracy of the classifier on the test set consisting of the data from Anzalone et al. and additional in-house validation experiments was 75%, with an area under the receiver operator characteristic curve (AUROC) of 0.90 (Figures 1B, C). Our findings thus support that heterochromatin may be sufficient to block PEs from targeting; however, heterochromatin is likely not necessary to hinder PE targeting efficiency, and other factors may also contribute to low edit rate, as discussed as follows.

Having confirmed that high-order chromatin accessibility affects PE efficiency, we next investigated the effect of local chromatin structure on prime editing. It has been previously reported both in vitro and in vivo that nucleosomes inhibit the efficiency of CRISPR-Cas9 endonuclease, suggesting that histone proteins may block the target DNA access to Cas9 (Hinz et al., 2015; Yarrington et al., 2018). Given the shared components between the CRISPR-Cas9 and PE constructs, we thus sought to test whether nucleosome positioning at the target sites might also inhibit PE editing efficiencies. Since there were no publicly available MNase-seq datasets for HEK293 and most nucleosome positioning was shown to possess some degree of consistency across different cell lines (Gaffney et al., 2012) and partially exhibit intrinsic DNA sequence preferences (Ioshikhes et al., 1996; Segal et al., 2006; Gupta et al., 2008; Jin et al., 2016), we used the available nucleosome occupancy data in the lymphoblastoid cell line GM12878 (Methods). We found a significant negative correlation between PE edit percentage and nucleosome occupancy in the endogenous target sites of Kim et al. and Anzalone et al. (Pearson r=−0.46,p=1.36×10−2; Figure 1D; Methods). While target sites in nucleosomal DNA tended to have low edit rates, target sites in nucleosome-free regions did not necessarily have high edit rates. Similar to the effect of heterochromatin, our findings thus support that nucleosomes are sufficient, but not necessary, to partially block PEs from targeting, as other factors may also contribute to low edit rate, as shown in the following sections.

Given that the PBS and spacer regions of pegRNA must hybridize to the target DNA sequence in order for editing to occur, with the spacer being the longer of the two sequences, we examined the off-target effect of having several genomic loci complementary to the spacer region on editing efficiency. We found only a weak but statistically significant negative correlation between the number of spacer off-targets and the edit rate averaged over all pegRNAs with the same spacer (Spearman ρ=−0.11; p=7×10−7; Figure 1E; Methods). This result indicates that the presence of off-targets does not substantially affect the editing efficiency at the intended target and is consistent with the previous finding for single-guide RNAs (sgRNA) of CRISPR-Cas9 (Moreno-Mateos et al., 2015).

Our aforementioned logistic regression model has revealed that target sites close to heterochromatin are unlikely to be editable. Once a target site away from heterochromatin is selected, the next important step is to optimize the pegRNA design by adjusting the lengths of PBS and RT templates. It was previously reported that altering the PBS length (PBSL) and RT template length (RTTL) could have a drastic effect on PE2 efficiency, even when targeting the same site for the same edit, and that having G as the last templated nucleotide tended to decrease the PE2 efficiency (Anzalone et al., 2019). We further observed in Anzalone et al.’s data that 1) at the EMX1 locus containing no Gs in the RT region (Anzalone et al.'s Figure 2B), having A as the last templated nucleotide consistently resulted in lower edit rates, suggesting that other nucleotides apart from G might also modulate the editing efficiency; 2) at the FANCF locus (Anzalone et al.'s Figure 2B), even though the last templated nucleotide was G at RTTLs of both 10 and 18 nts, the edit rate for 18 nts was much higher than that for 10 nts, suggesting that the presence of G alone could not explain the variability of edit rates as a function of RT template sequence composition. These findings motivated us to use the sequence content in the target region as features in predicting the PE2 efficiency as PBSL and RTTL were varied. We thus developed an ordinary least squares (OLS) linear regression model to predict stepwise differences in edit percentages between two adjacent RTTLs, using the last templated nucleotides as features (Methods).

Anzalone et al. generated edit percentages (defined as the fractions of total reads with correct edits) of PE2 at eight different gene loci using pegRNAs with varying ranges of RTTL for each locus, with the largest overlap of lengths among the designs being between 10 and 20 nts. We thus chose this range [10 nt, 20 nt] to train an OLS linear regression model for predicting the stepwise differences, ΔEPi=EPi+1−EPi, in edit percentages between i+1 and i RTTL (Figure 2A; Pearson r = 0.79; Methods). Eight-fold cross-validation of holding out and testing on each gene locus supported its robustness, with the root mean square error (RMSE) of the test set being mostly similar to that of the training set (∼2.3%) (Supplementary Table S4). Our model learned the effect of the last templated nucleotide on ΔEPi to be in increasing order of G, A, C, and T for the (i+1)^th templated nucleotide and decreasing order of G, A, C, and T for the i^th templated nucleotide (Supplementary Table S5). The relative ordering of G < A < C < T resembles the order of ionization energy of the nucleotides (Burrows and Muller, 1998), but clarifying the connection would require further experimental and theoretical investigation.

Having confirmed that our OLS linear regression model performed well on the eight target sites from Anzalone et al., we further investigated whether a similar approach could be generalized to Kim et al.’s data (Kim et al., 2021), which measured editing efficiencies at thousands of target sites by high-throughput sequencing. The PBSL was set to 13 nts, to be consistent with the range used in Anzalone et al. Since Kim et al. measured the edit percentages only for RTTLs = 10, 12, 15, and 20, we could not directly apply our stepwise approach on these data. We thus divided the data for different RTTLs into three ranges and trained an OLS linear regression model to predict consecutive stepwise differences in edit percentages for each range separately (Methods): 1) RTTLs from 10 to 12; 2) RTTLs from 12 to 15; and 3) RTTLs from 15 to 20 (Supplementary Figure S3). Ten-fold cross-validation verified the robustness of our model within each range (5.1% RMSE; average of Pearson correlation across folds = 0.33 for RTTLs in the range [10,12] and [15,20], and 0.25 for RTTLs of [12,15]; Supplementary Table S6). However, the effect sizes of the four nucleotides differed somewhat across the three ranges (Supplementary Tables S7–S9). When the RTTL is in the range [15,20], the effect on ΔERi was in the increasing order of G, A, T, and C for the (i+1)^th templated nucleotide and decreasing order of G, A, C, and T for the i^th templated nucleotide (Supplementary Table S9), similar to the aforementioned results for Anzalone et al.’s data. RTTL in the [12,15] range also showed the pattern that G as the (i+1)^th templated nucleotide was associated with the smallest ΔEPi and G as the i^th templated nucleotide was associated with the largest ΔEPi compared to other nucleotides at these respective positions (Supplementary Table S8). This pattern was reversed for the RTTL in the [10,12] range, where G as the i^th templated nucleotide was associated with the smallest ΔEPi compared to other nucleotides at this position (Supplementary Table S7).

Testing the regression model trained on Kim et al.’s data on Anzalone et al.’s data in the corresponding RTTL ranges (Supplementary Tables S7–S9), the predicted ΔEPi still significantly correlated with the observed values, albeit to a lesser extent than the model directly trained on Anzalone et al.’s data (Pearson r = 0.37 ,p=7.13×10−4; Figures 2A, B). Repeating the same analysis in individual RTTL ranges revealed that the predicted ΔEPi was significantly correlated with the observed values only in the 15–20 RTTL range (Pearson r=0.25, p=0.34 in the 10–12 RTTL range; Pearson r=0.05, p=0.81 in the 12–15 RTTL range; Pearson r=0.62, p=2.12×10−5 in the 15–20 RTTL range). We also predicted the trend of absolute edit percentages in Anzalone et al.’s data as a function of the RTTL (Figure 2C). Overall, the model trained on Anzalone et al.’s data was able to reproduce the observed trend that the edit percentage decreased whenever the last templated nucleotide was G (green line in Figure 2C). The models trained on Kim et al.’s data were able to capture the observed trend in longer RTTLs, but not so well in shorter RTTLs (magenta line in Figure 2C), consistent with the fact that our predicted ΔEPi was more accurate in the long RTTL range [15,20].

There were some discrepancies between the two independent datasets regarding the regression coefficients trained on shorter RTTL designs. For example, Anzalone et al. observed a decrease in editing efficiency whenever the last templated nucleotide was G throughout all RTTLs. By contrast, Kim et al. observed that, on average, editing efficiency was highest when the last templated nucleotide was G at RTTLs = 10 and 12, while it was lowest for the last templated G only at RTTL = 20 (Kim et al.'s Figure 2F). Further investigation is needed to understand why G has an opposite effect on editing efficiency at short RTTLs and why this phenomenon is not universal across independent datasets. We shall revisit the position-specific effects of G on editing efficiency in subsequent sections.

The OLS linear regression analysis yielded insight into the preferred RTTLs of pegRNAs for a given target site, but it could not reveal the optimal target, given multiple PAM candidates in the region of interest. We therefore built an elastic net model (Zou and Hastie, 2005), a linear regression model with l1 shrinkage to impose sparsity of features and l2 shrinkage to reduce overfitting, for each combination of PBSL and RTTL, separately (PBSLs = 7, 9, 11, 13, 15, and 17 and RTTLs = 10, 12, 15, and 20) (Kim et al., 2021) and predicted the absolute editing efficiency of a target site using only the sequence information of the PBS, RT template, and its flanking regions (Methods). The predictive variables were the indicator variables encoding the four nucleotides at each position of the 47-bp-long target site sequence, except for the fixed GG dinucleotide in the NGG PAM (Methods). Since this dataset did not contain information about where the target sites were integrated in the genome, features involving relative distance to epigenomic modifications were not included in our model. Our approach was much simpler and more interpretable than the previously reported model pooling together different PBSL and RTTL designs and using 1,766 features (Kim et al., 2021), many of which might have been highly correlated and redundant.

The elastic net predictions significantly correlated with the observed edit percentages in each pair of PBSL and RTTL (Pearson r>0.6; p<10−150; Supplementary Table S10). Figure 3A shows the regression coefficients for PBSL = 13 and RTTL = 15 (other combinations of PBSL and RTTL are given in Supplementary Figure S4). In general, the higher the GC content in the PBS region, the more likely the target site was to have high editing efficiency (Supplementary Figure S5), consistent with the previous report (Kim et al., 2021). However, the extent of this preference depended on the PBSL, that is, the high-performance designs with shorter PBSL required higher GC content in the PBS region than those with longer PBSL (Supplementary Figure S4). Considering that G: C bonds are stronger than A: T bonds, this result suggested that a short PBS lacking a sufficient number of G/C nucleotides may not stably hybridize to the 3’ end of pegRNA, thereby yielding low editing efficiency. To further investigate this idea, we computed the length-normalized RNA-DNA hybridization energy between the pegRNA PBS region and the corresponding DNA and found a significant negative correlation between PBS-DNA hybridization energy and edit percentage for all combinations of PBSLs and RTTLs (Figures 3B, C; Supplementary Figure S6; Methods). Additionally, we confirmed again that the shorter the PBSL, the larger the magnitude of the anti-correlation between editing efficiency and PBS-DNA hybridization energy (Figure 3B vs. Figure 3C; Supplementary Figure S6).

The regression coefficients also showed that the presence of G as the last templated nucleotide tended to increase the editing efficiency when RTTL = 10, regardless of the PBSL, or when RTTL = 12 and PBSL ≤ 9 (Figure 3A; Supplementary Figure S4). This finding recapitulated the observation made by Kim et al. that the editing efficiency was, on average, highest when the last templated nucleotide was G for RTTL = 10 or 12 (Kim et al., 2021); however, our analysis clearly highlighted that this effect depended on the PBSL being sufficiently short for the case RTTL = 12. By contrast, when PBSL ≥ 11 and RTTL ≥ 15, the presence of G as the last templated nucleotide decreased the editing efficiency (Figure 3A; Supplementary Figure S4), consistent with the findings of our aforementioned OLS model and Anzalone et al. (2019); for shorter PBSL, however, the last templated nucleotide had only a minor effect on editing efficiency, and the GC content in the PBS region instead had a pronounced effect.

The large effect size of the GC content in the PBS region for short PBSL designs, together with the accompanying reduction in the negative effect of the last templated nucleotide, suggested that optimal choices of pegRNA for target loci containing high GC content in the PBS region would involve short PBSLs. Supporting this idea, we observed that when the pegRNA contained only G/C in the range [−7,−1] and A/T at −9 and −8, the edit percentage was, on average, highest for PBSL = 7 (Figure 3D; highest median edit percentage of 20.5% when PBSL = 7); here, we considered A/T at the −9 and −8 positions based on the regression coefficients learned by the elastic net for PBSL = 7 (Supplementary Figure S4). When fewer than 6 G/Cs were found in the [−7,−1] range, PBSLs between 11 and 15 had, on average, higher editing efficiency than other PBSLs (Figure 3D; highest median edit percentage of 9.8% when PBSL = 13). Given a target site, these findings thus provided a general guideline for determining the optimal PBSL based solely on the 9-bp-long sequence upstream of the nick site in the edited strand.

Some patterns of regression coefficients were shared among most combinations of PBSL and RTTL. Based on the magnitude of regression coefficients, the most salient nucleotides were those positioned around the nick site (Supplementary Figure S4): C and G nucleotides at locations from −1 to +3, immediately 5’ to the PAM site, led to high editing efficiency, whereas T in the range from −2 to +3 was to be avoided. G/A at the −17 position was positively correlated with editing efficiency, while T at the same position was negatively correlated. In general, C at the +20 position was positively correlated with editing efficiency across designs. For the PAM site, AGG was favored and CGG disfavored for optimizing editing efficiency. Immediately after the PAM site, G at the +7 position was strongly anti-correlated with editing efficiency. Surprisingly, the presence of G at the +10 position was significant and tended to increase editing efficiency, not only for RTTL = 10, in which case G would be the last templated nucleotide as described previously, but also for RTTL = 15 or 20 when PBSL ≥ 11 (Figure 3E). Considering that Cas9 variants similar to SpCas9 in size were known to recognize PAM sequences of varying length between 2 and 8 nts (Anders et al., 2014) and that the +10 position was only 4 bps downstream of the NGG PAM site, the fact that the presence of G at the +10 position consistently increased the editing efficiency across multiple RTTLs suggested that SpCas9 preferentially interacted with this specific nucleotide. Another unexpected finding was that select edited strand nucleotides even outside the regions of the protospacer, PBS, and RT template seemed to modulate editing efficiency. For example, when PBSL = 13 and RTTL = 15, the presence of A at the +17 position, lying outside the range of the RT templated region, tended to decrease the editing efficiency (Figure 3F).

In summary, the elastic net regression models trained on each combination of PBSL and RTTL learned distinct and common sequence features modulating editing efficiency. The effect of the last templated G nucleotide on editing efficiency depended on both PBSL and RTTL of pegRNA designs. Moreover, pegRNAs with short PBSLs highly depended on the GC content of the PBS region, perhaps to help stabilize the hybridization to the complementary DNA. The PBS RNA-DNA hybridization energy was significantly anti-correlated with edit percentage in all combinations of PBSL and RTTL, but pegRNAs with shorter PBSLs showed more pronounced dependence. Finally, the presence of G/C nucleotides was most critical immediately 5’ to the PAM site for all combinations of PBSL and RTTL.

It is challenging to model nonlinear effects of coupled nucleotides at multiple positions using elastic net regression. To extend our linear regression approach, we thus trained a DNN model on Kim et al.’s data (Kim et al., 2021) to predict the edit percentage of a pegRNA design given its length and the sequence of the target site and flanking regions (Figure 4A; Methods). The Pearson correlation coefficient between the observed and predicted edit percentages in the test set was 0.73 (Figure 4B; Methods).

To understand optimal sequences associated with high edit percentages, we applied a simulated annealing (SA) method for maximizing the prediction of a trained DNN over its input space via Markov Chain Monte Carlo (MCMC) sampling (Finnegan et al., 2020). We accelerated the simulation by using a broadened sampling distribution stemming from nonextensive statistical mechanics (Tsallis and Stariolo, 1996) and obtained the target sequences maximizing the predicted edit percentages for each combination of PBSL and RTTL (Supplementary Method S3; Supplementary Tables S11–S14). To identify the DNN-learned salient features in these optimal sequences, we next applied MaxEnt, another MCMC sampling method based on the maximum entropy principle, for generating new input sequences that produce similar DNN predictions as those of the initial input sequence (Finnegan and Song, 2017). Upon initializing the MaxEnt chains at the optimal target sequences, nucleotide preferences important for maximizing PE efficiency were extracted by calculating the position-wise Kullback–Leibler (KL) divergence between the nucleotide distribution in the sampled sequences and the uniform null distribution; in this formalism, the larger the KL divergence, the more relevant the nucleotide (Figure 4C; Supplementary Figure S7; Supplementary Method S4).

The DNN confirmed the sequence features previously detected by the elastic net model. For example, it learned that C at the +1 position, G at the +3 position, and G at the −17 position relative to the nick site were associated with high edit percentage; other similarities with the elastic net results included the preferences for A as the first nucleotide in the PAM sequence, for C or G in the PBS region, and for C at the +20 position. In addition, the DNN also discovered new features. For example, although the elastic net model detected no strong preference for any particular nucleotide at position +9, except for perhaps a weak preference for G (Supplementary Figure S4), the DNN found T at this position to be associated with optimal editing. This discrepancy in the feature importance arose from the difference in model architectures, that is, the elastic net model learned only the independent effects of individual nucleotides, whereas the DNN accounted for aggregate effects of nucleotides captured by the filters. When considering only the marginal effect of single nucleotides at the +9 position, target sequences with T at that position actually had a lower edit percentage on average than those with other nucleotides at that same position (9.00% for T vs. 9.72% for A, C, and G at +9; Figure 4D). However, imposing the +8 position to be T and the +10 position to be G showed that target sequences with T at +9 had a higher edit percentage on average than those with other nucleotides at +9 (14.33% for T vs. 10.72% for A, C, and G at +9; Figure 4E; Supplementary Figures S8A, B). Similarly, the marginal effect of T at the −4 position, residing in the PBS region, was negligible; however, the DNN uncovered a role of T in the context of a GC-rich background: for sequences with 100% GC content in the range [−7,−1], except at the −4 position, the presence of T at −4 substantially increased the edit percentage by 6.90% on average (Supplementary Figures S8C, D). Finally, the DNN found a preference for CAG at [+1,+3] in optimal target sequences, whereas the elastic net model found only a marginal effect of G/C at +2 and a less pronounced trinucleotide effect of C (G/C) G in the same region (Supplementary Figures S8E, F). In summary, our DNN was able to predict the observed edit percentage with high accuracy and yielded interpretable results regarding both marginal and context-dependent sequence features of optimal target sites.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/↗, PRJNA949853.