Massively parallel profiling and predictive modeling of the outcomes of CRISPR/Cas9-mediated double-strand break repair

1Genome engineering conventionally involves using a programmable endonuclease (i.e. a zinc finger nuclease (ZFN), transcription activator-like effector nuclease (TALEN) or RNA guided nuclease Cas9 (clustered regularly interspaced short palindromic repeats-CRISPR associated protein CRISPR/Cas9) to introduce a double-strand break (DSB) at a specific location in the genome. In mammalian cells, such DSBs are primarily repaired by one of two pathways—homology directed repair (HDR) and classical non-homologous end joining (c-NHEJ). HDR uses homologous template sequences to repair the DSB, potentially introducing programmed edits via the repair template. In contrast, c-NHEJ directly rejoins the broken ends, often perfectly but occasionally introducing errors, typically in the form of short insertions or deletions (indels). In addition to HDR and cNHEJ, there is evidence for an alternative NHEJ pathway (alt-NHEJ), also termed microhomology mediated end joining (MMEJ), wherein short, homologous sequences in the vicinity of the DSB are used to align the broken ends prior to joining, resulting in deletions or potentially more complex events. Below, we use ‘NHEJ’ to refer to both c-NHEJ and MMEJ/alt-NHEJ, i.e. template-free editing.

2In recent years, CRISPR/Cas9 has emerged as a particularly versatile tool for genome editing. For many if not most applications of CRISPR/Cas9-mediated genome engineering, it is used in conjunction with the cell's endogenous NHEJ machinery to introduce short indels in a targeted fashion, e.g. to disrupt the function of genes or regulatory elements or to introduce irreversible changes that record cell lineage or molecular events. However, despite NHEJ’s central role in this transformative tool, our understanding of the processes that determine the rate and patterns of NHEJ-mediated errors remains incomplete.

3Recent studies have demonstrated that the error outcomes of NHEJ are strongly dependent on sequence context. Other studies show that the characteristics of the broken ends (blunt or staggered end; length of any overhang) also affect end-joining patterns both in vitro and in vivo. However, a systematic profiling of the sequence determinants of NHEJ repair patterns has yet to be undertaken.

4Here, we profiled ∼1.16 million mutational events resulting from Streptococcus pyogenes Cas9 (SpCas9)-mediated cleavage and NHEJ-mediated DSB repair of 6872 synthetic target sequences. From the resulting data, we identify the primary features of sequences adjacent to the sites of DSBs that shape the distribution and relative frequency of NHEJ-mediated mutational outcomes, e.g. nucleotide content and microhomology. We furthermore exploit microhomology to demonstrate the ‘programming’ of deletion patterns. Finally, we develop a logistic regression model to predict insertions and deletions (Lindel) that result from CRISPR/Cas9-mediated cleavage of an arbitrary sequence. A standalone Lindel webtool is freely available (https://lindel.gs.washington.edu), and Lindel predictions have been integrated into the CRISPOR web tool (http://www.crispor.org).

5To generate a library of CRISPR/Cas9 targets that could safely be characterized within human cells, we evaluated ∼1 million random 20mer crRNA sequences, scoring them against the human genome (version hg19) for off-target effects using FlashFry. We excluded guides with an exact match or up to two mismatches against any potential target in the human genome, or those with an off-target score <90, resulting in a modest bias towards targets containing CpG dinucleotides (Supplementary Figure S7), and then selected a final library of 70 000 top scoring guides for synthesis. The resulting sgRNA sequence and their corresponding targets were separated by a common 20 bp spacer sequence and ordered as an Agilent SureGuide Unamplified Custom CRISPR Library array (Figure 1A).

6To analyze the potential impact of programmed microhomology, we selected a subset of 1000 sgRNA-target pairs from the library above and introduced microhomology with different lengths (2, 4 and 6 bp) matching the last 2, 4 and 6 nucleotides upstream of the cleavage site. Each design was assigned a 4 bp barcode, indicating its programmed microhomology pattern (Figure 5A and B). This library of microhomology sequences was ordered as an oligo pool from Twist Biosciences.

7The lentiGuide-Puro (Addgene #52963) vector was modified with two rounds of PCR to remove the existing tracrRNA and filler sequence (primer P1, P2), and to incorporate two BsmBI restriction site for integration of sgRNA-target pairs (primer P3, P4). The modified vector was digested with BsmBI (NEB, Buffer 3.1) at 55°C for 3 h and gel purified with Monarch DNA Gel Extraction Kit (NEB). This digested and purified vector was used for all downstream cloning.

8Oligos with sgRNA-target pairs from Agilent or Twist Bioscience were both resuspended to 10ng/μl. The oligo pool was PCR amplified using KAPA Biosystems HiFi HotStart ReadyMix 2× using primers P5 and P6 and cleaned with the DNA Clean&Concentrator kit (Zymo Research). The purified PCR product was then digested with BsmBI (NEB, buffer 3.1) at 55°C for 1 h to generate compatible sticky ends matching the modified lentiGuide-Puro above, and subsequently cleaned with DNA Clean&Concentrator (Zymo Research). Digested vector and insert were ligated with T4 ligase (NEB) with a molar ratio of 1:3. Ligation products were transformed in to Stable Competent Escherichia coli (NEB C3040H). Transformed cells were cultured at 30°C overnight and plasmid DNA was prepared using a ZymoPURE II Plasmid Kit. The subsampled library with 12 917 targets was bottlenecked by seeding transformed cells on plate. Colonies on plates were transferred to liquid medium to expand them. The precision of the number 12 917 follows from the fact that we can simply count the number of unique guide sequences present in deep sequencing of PCR amplicons from cells.

9We generated a mono-clonal 293T cell line expressing Cas9 by transduction of Cas9-blast lentivirus particles (Addgene plasmid #52962). Cells were cultured in DMEM High glucose (GIBCO) supplemented with 10% fetal bovine serum (Rocky Mountain Biologicals) and 1% penicillin­–streptomycin (GIBCO) and grown with 5% CO2 at 37°C.

10All lentivirus libraries were produced by the Fred Hutchinson Co­operative Center for Excellence in Hematology Vector Production core facility. HEK293T cells were transduced and media was changed to virus free media for 24 h post­-transduction. Cells were passed every 48 h with a split ratio of 1:6. Cells were harvested at day 5 after transduction.

11Genomic DNA was extracted with DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer's protocol. 15 bp unique molecular identifiers (UMIs) were added by one initial round of linear PCR using a primer containing a 5′ sequencing adaptor (P7). For each reaction we used 250 ng of genomic DNA, 0.2 μl 100 mM primer and 25 μl HiFi HotStart ReadyMix 2× (KAPA Biosystems). PCR reaction were performed as follows: 95°C 3 min, 98°C 20 s, 5 cycles of 65°C 1 min and 72°C 2 min, 98°C 20 s, 5 cycles of 65°C 1 min and 72°C 2 min. The subsequent PCR product was cleaned with 1.8× AMPure XP beads (Beckman Coulter) and resuspended in 25 μl of elution buffer. A second round of amplification was performed using primers targeting the 5′ sequencing adaptor (P8) and 50 bp downstream of the cleavage site (P9) for 20 cycles. The resulting PCR product was then size selected using a dual size-selection cleanup of 0.4× and 0.8× AMPure XP beads (Beckman Coulter) to remove genomic DNA and small fragments (<200 bp) respectively. This size-selected product was subsequently re-amplified to add the 3′ sequencing adaptor with primer P8 and P10 for an additional five cycles. The final PCR product was cleaned with 0.75x AMPure XP beads (Beckman Coulter) and was re-amplified to add flow-cell adaptor and sample index for 5 cycles. All PCR reactions used HiFi HotStart ReadyMix 2× (KAPA Biosystems) with the manufacturer's recommended conditions. The library was sequenced on an Illumina NextSeq 500 sequencer using paired-end 150 cycle reads. All primers used are listed in Supplementary Table S2. Sequencing data has been deposited in the Gene Expression Omnibus (GEO) with accession number GSE131421 and the processed data files are deposited in Figshare with a doi link: https://doi.org/10.6084/m9.figshare.7374155.

12Across three replicates, we sequenced a total of 148 million paired-end reads on an Illumina NextSeq 500. We first clustered these paired-end reads by their 15 bp UMI sequence and then filtered out reads with <90% identity within their representative UMI clusters. Sequence identity was identified using edlib. UMIs with fewer than 10 reads were excluded from downstream analysis. This yielded 4 405 379 UMIs (91 325 700 reads), representing ∼61.8% of our sequencing data (Supplementary Table S3). We then selected the most common forward and reverse read sequence for each UMI for further processing. These forward and reverse reads were merged into a single read using PEAR and aligned in a two step process as follows. First, we sought to identify the ‘reference’ sequences for each programmed array sequence. We aligned the merged reads to a backbone sequence where the guides and targets were represented by Ns using EMBO’s needleall software with the following scoring matrix: match = 5, mismatch = -4, gap-open = –20, gap-extension = –0.5. The mismatch penalty for Ns was set to 0. The sequence over the guide region was then extracted and matched against the list of programmed array sequences. Guide sequences with more than two mismatches to the designed guides were excluded, with edit distances assessed with UMI-tools. Second, merged reads were aligned to their discovered reference, in which Ns were replaced by the guide/target sequence identified from the first step, using Biopython.pairwise2 with the following scoring matrix: match = 5, mismatch = –4, gap-open = –13, gap-extension = –0.5. All indels were then right aligned (e.g. Supplementary Figure S3A). Aligned reads with indels within –3/+2 bp of the cleavage site were assigned to their indel class. Aligned reads were excluded for downstream analysis if the sgRNA and target sequence didn’t match, the result from template switch during lentivirus transduction, or unexpected mutations introduced during synthesis, cloning, and PCR. A final library of 1.19 million unique reads (UMIs) were identified. Our library of 1000 microhomology sequences were processed by this same pipeline, yielding a final library of 249 039 UMIs from 31 239 645 paired-end sequencing reads. Scripts and other software are available from our GitHub repository: https://github.com/shendurelab/Lindel.

13Sequence motif analysis was conducted with kpLogo using default settings with a specified k-mer length of 1 or 2. Input sequences were weighted by the frequency of insertion.

14For n from 1 to 10 nucleotides, the last n nucleotides upstream of each deletion were compared to the last n nucleotides of the deleted sequence (as the deletion is right aligned. Supplementary Figure S3A). The length of microhomology was identified as the largest n nucleotides match in sequence.

15We phrased our problem of predicting repair outcomes and their frequencies as that of a classification task with 557 classes. Because large mutation events are rare, we limited our classification effort to deletion events <30 bp, and we grouped insertions ≥3 bp into one class. In total, we defined 557 classes of indels. These classes include 536 deletion alleles, 4 possible single nucleotide insertion, and 16 possible dinucleotide insertion and insertions ≥3 bp. There are a total of 550 potential deletion events that are both <30 bp in length and overlap with the –3/+2 window around the cleavage site. We captured 536 deletion alleles in our data; the missing 14 classes are mainly large deletions. As input to our model, we defined 3033 binary features. These are (i) Sequence features: 384 binary features corresponding to the one-hot encoded target sequence (excluding the PAM region), including 80 for single nucleotide content (4 nucleotides × 20 positions) and 304 for dinucleotide content (16 dinucleotides × 19 positions); (ii) Microhomology features: 2649 binary features corresponding to MH tracts; specifically, for each of the possible deletion event class, we defined five binary features (or 2–4, depending on the size of deletion) corresponding to the length of the MH tract, if any (0-4 bp × 519 + 0–3 bp × 7 + 0–2 bp × 6 + 0–1 bp × 4 deletion event classes = total 2649 binary features. Our 4790 programmed sequences were randomly partitioned into a training set of 3,900 sequences, a validation set of 450 sequences, and a test set of 440 sequences.

16We trained the logistic regression in a standard manner for machine learning models. However, because each target sequence can generate many possible repair outcomes, we trained our models using soft labels that correspond to the probability that each class is observed, rather than hard labels that force each input to correspond exclusively to one class. Each model was trained using the Adam optimizer with a learning rate of 0.001 and a categorical cross-entropy loss. Training proceeded for a maximum of 100 epochs with a ‘patience’ of 1, meaning that training was stopped after two epochs with no improvement in validation set performance. All initializations and the hyperparameters for the Adam optimizer were set to the defaults in Keras v2.1.3 with a backend of Theano v1.0.1. We selected the best model based on performance on the validation set according to the coefficient of determination using grid search over hyperparameters. This search involved separate scans over regularization strengths for L1-regularization and L2-regularization individually with a range of 10−10 to 10−1 (Supplementary Figure S5D–F).

17We compared our model to two other models (ForeCasT and inDelphi) in the same setting. All models used 60 bp sequence centered at the cleavage site as an input, while trying to predict the frequency of 557 classes of indels that we defined above. As both Lindel and ForeCasT are predicting all possible unique repair outcomes, it's straightforward to them compare directly. inDelphi only predicts ∼90 classes deletions with locations, and 1–59 bp deletions regardless their locations. We used the classes that inDelphi predicts that overlap with the ForeCasT and Lindel classes, which included all ∼90 classes using microhomology, 1 bp deletion (assuming it's located the cleavage site), and 1 bp insertions. All classes that inDelphi is not predicting were assigned as 0. The performance were measured using MSE on all ∼450 unique classes for each sequence in our test set (n = 440) and the ForeCasT test set (n = 4298)

18Toward a comprehensive understanding of the sequence determinants of NHEJ-mediated error patterns, we developed a strategy that would allow us to efficiently profile a large number of repair events from each of a large number of sequence contexts (Figure 1). In brief, we designed 70 000 targets balanced in nucleotide content and screened against the human genome for CRISPR/Cas9 single guide RNAs (sgRNAs). We then used array-based oligonucleotide synthesis to encode these targets in cis with their corresponding sgRNAs, separated by a 20 bp spacer. We then amplified and cloned these molecules to a lentiviral vector. In our initial experiments, the complexity of the resulting library of synthetic targets and their cognate sgRNAs was such that we obtained relatively few edited templates per target. Therefore, we re-cloned the library under bottlenecking conditions (Materials and Methods), reducing its complexity to 12 917 targets. We then proceeded with viral packaging and transduction, in triplicate, of a monoclonal human embryonic kidney (HEK) 293T cell line that stably expresses spCas9 (multiplicity of infection of ∼4–8). As such, within any given cell, only one or a few sgRNAs are expressed, and each one directs Cas9-mediated DSBs to a target located immediately adjacent to it. After five days to allow for the introduction of NHEJ-mediated errors at these targets, cells were harvested and genomic DNA isolated. We then PCR amplified the region comprising the targets and corresponding sgRNAs using unique molecular identifiers (UMIs) appended during the first extension cycle to distinguish whether identical edits were derived from the same cell or different cells.

19Summing across the three replicates, we sequenced PCR amplicons to a depth of ∼148 million reads, which were reduced to ∼1.19 million reads after collapsing on the basis of identical sequences and UMIs, and filtering of reads with evidence of lentivirus-mediated template switching or other unexpected sequences (e.g. synthesis or PCR errors). After further filtering of poorly represented targets (those represented by fewer than 10 UMIs), our dataset consisted of ∼1.16 million UMIs corresponding to 6872 unique targets. On average, each target was represented by 168 UMIs and 24 alleles (where ‘allele’ refers to a unique post-editing sequence of a given target). Each allele was aligned to its original sequence, known because the corresponding gRNA sequence is part of the same amplicon, using the Needleman-Wunsch algorithm. Alleles were categorized as wild-type (i.e. unedited), a deletion, or an insertion.

20Overall, targets were highly edited, with only 9.8% of UMIs corresponding to the wild-type allele. Of UMIs containing detectable mutations, 63.6% were deletions and 31.5% were insertions (Figure 2A). The remainder (4.9%) contained some combination of substitutions, insertions and deletions, and are excluded from all of our subsequent analyses. Deletions were dominated by small events; only 1.5% were >25 bp, although we note that deletions >150 bp are not captured by our assay. In contrast, although we believe that our assay should have been able to recover insertions up to ∼500 bp, the overwhelming majority of insertion events were of a single base pair.

21We sought to examine whether repair patterns for any given target were reproducible, as previously shown for a more limited set of templates. For each target, we calculated the frequency of each non-wild-type allele. For any given target, the distribution of frequencies for its alleles were highly reproducible in pairwise comparisons of the three replicates (median Pearson's r = 0.91, 0.93, 0.93, Figure 2B, left). Meanwhile, if we permute the alleles in one replicate on a target-by-target basis and repeat the pairwise comparison, these correlations are greatly reduced (median Pearson's r = 0.20, Figure 2B, right).

22Confirming the observations of, the diversity of mutations strongly varied from target to target. We calculated the Shannon entropy of mutational outcomes for any given target as –∑ pi*log(pi), where pi is the frequency of ith indel of that target (Figure 2C). Entropy values for any given target were highly reproducible between replicates (Figure 2D) and only modestly correlated with sampling depth (Supplementary Figure S1). Of note, some targets consistently exhibited particularly diverse mutational outcomes consequent to NHEJ—that is, high entropy (e.g. Figure 2E, where the most frequently observed mutation occurs in only 10.1% of mutated templates). Other targets were strongly biased towards a more limited set of mutational outcomes—that is, low entropy (e.g. Figure 2F, where the most frequently observed mutation occurs in 80.4% of mutated templates).

23We next sought to investigate the determinants of insertions at the DSB, which were dominated by 1 bp events (Figure 3A). 84% of 1 bp insertions were predicted (and presumably templated) by the nucleotide immediately upstream of the cleavage site (i.e. the 17th nucleotide in target sequence; Figure 3B; Supplementary Figure S2) Although it might have been expected that NHEJ-mediated repair would be symmetric with respect to the site of a DSB, we do not observe templating from the immediately downstream (18th) nucleotide (Figure 3B). Similarly, of 2 bp insertions, a substantially greater than expected proportion (41%) were templated by the sequence immediately upstream of the DSB (i.e. inserted sequence identical to the 16th and 17th nucleotides of the target sequence; Figure 3C). The asymmetric templating of NHEJ-mediated insertions was also described in two other recent studies based on data from yeast and mice.

24Because the ratio of insertions to deletions varied from target to target, we used kpLogo to examine what local sequence features might shape this. We find that the presence of a T or A at the 17th bp of the target was associated with insertion events, while a G or C at this position was associated with deletion events (Figure 3D, left). Additional analyses showed a ‘TG’ dinucleotide flanking the cleavage site to be the most highly biased toward insertion (57% of events with that context are insertions), while a ‘GA’ dinucleotide flanking the cleavage site was the most highly biased towards deletion (17% of events with that context are insertions) (Figure 3D, right).

25We split 2680 targets associated with both insertion and deletion outcomes into training (n = 2000) and test (n = 680) sets, and trained a linear regression model to predict the proportion of insertion events based on position-specific content of the hexamer centered on the DSB (single and dinucleotide k-mers; 104 binary features; Figure 3E). The model performs reasonably well (Pearson's r = 0.70).

26Overall, these analyses confirm that local sequence around the DSB site plays an important role in shaping the outcome(s) of NHEJ-mediated errors. In particular, the asymmetry implied by the high rate of identity between 1 and 2 bp insertions and the nucleotides immediately upstream to the DSB, but not the nucleotides immediately downstream to the DSB (Figure 3B-C), suggests that not all CRISPR/Cas9-mediated cleavages are blunt-ended. Indeed, in vitro studies have shown that the non-complementary strand of the target can sometimes be cleaved by Cas9 at multiple sites upstream of the –3 bp position relative to the protospacer adjacent motif (PAM), while the complementary strand is cut only at that site, instances which would result in a 5′ overhang. The preponderance of 1 bp insertions templated by the 17th rather than 18th base could be explained by fill-in of this overhang followed by blunt-ended ligation (and similarly for the preponderance of 2 bp insertions that are templated by the 16th and 17th bases, rather than the 18th and 19th bases).

27To summarize, we propose a model (Figure 3F) where: (i) some proportion of cleavages of the non-complementary strand by Cas9 occur upstream of the –3 bp PAM cleavage site, while cleavage of the complementary strand always occurs between the 17th and 18th positions, resulting in a 5′ overhang; (ii) 5′ overhangs are preferably repaired by gap-filling and ligation, resulting in the observed bias toward templating by the bases immediately upstream rather than downstream of the DSB; (iii) local sequence context biases the pattern of cleavage on the non-complementary strand, resulting in different frequencies of blunt versus 5′ overhangs for different targets, which in turn biases the ratio of insertions versus deletions. A similar model was recently proposed by Lemos et al. based on asymmetric templating of NHEJ-mediated insertions observed in yeast.

28We next examined patterns of deletion. Microhomology (MH) refers to the use of short regions of identical sequence (1–16 bp) that can mediate the alignment of broken ends (Figure 4A) and is relevant to both c-NHEJ and alt-NHEJ/MMEJ. Here, a deletion event is considered to be MH-mediated if the sequence at the 3′ of a rejoined end is identical to the 3′ end of the deleted sequence, and the size of the MH tract refers to the length of that identical sequence. By that definition, we found that over 75% of deletion events in our dataset are MH-mediated. The length of MH tracts ranged from 1 to 10 bp. Nearly all MH-mediated events (94.6%) involved relatively short tracts of microhomology, i.e. 1–4 bp. Longer MH tracts were observed more rarely (Supplementary Figure S3A), probably simply due to the relative paucity of opportunities in our set of target sequences.

29The frequencies of tracts of various lengths consistent with MH usage were substantially higher than background expectation for all lengths except 1 bp, with that frequency increasing as a function of tract length (Figure 4B). We further investigated the relevance of 1 bp MH by comparing the proportion of 1 bp deletion events in targets with identical versus non-identical nucleotides immediately spanning the cleavage site. We observe a 3-fold greater proportion of 1 bp deletion events when those nucleotides are identical than when they are not (Supplementary Figure S3B), suggesting that 1 bp MH may play a role in aligning, stabilizing and rejoining the broken ends.

30The lengths of MH versus non-MH mediated deletions exhibited distinct distributions (Figure 4C-D). In particular, the distribution of deletion sizes for MH-mediated events peaks at both 1 and 5–6 bp, while an equivalent distribution for non-MH-mediated deletions peaks at both 1–2 and 8 bp. The frequency of longer deletions exhibits an exponential decay for both MH and non-MH mediated events. To investigate this further, we jointly analyzed the frequency of start and end points for deletion events, relative to the position of the canonical cleavage site (Figure 4E and F). Both MH and non-MH mediated deletions exhibited a preference for ‘unidirectional’ events, i.e. either the start or end point is immediately adjacent to the cleavage site, rather than the deletion spanning the cleavage site.

31What explains the excess of deletion events of specific lengths? For MH-mediated events, the excess of 1 bp deletions may simply be attributable to the aforedescribed instances of identical nucleotides spanning the cleavage site (Supplementary Figure S3B). However, the excess of 5–6 bp MH-mediated events is clearly driven by events in the downstream direction (Figure 4E), i.e. deletions between the DSB and the PAM. A potential explanation is that the predilection of PAM-like sequences near the DSB for deletion events (i.e. a G nucleotide at the 17th position or a CG dinucleotide at the 16th/17th position; Figure 3D), coupled with the consistent presence of the CGG PAM sequence at the 21st–23rd position, results in an excess of deletions mediated by CG (5 bp deletion) or G (5-6 bp deletion) microhomology (Figure 4G). Further work would be required to confirm this, as there may be other explanations.

32For non-MH-mediated events, the excess of 8 bp events might be explained by the observation that in the dsDNA-sgRNA-Cas9 complex, the region 1–8 bp downstream of the cleavage site is occupied by Cas9, even after cleavage. Thus, the enrichment of non-MH deletions 8 bp from the cleavage site could simply correspond to the nearest position lacking Cas9 protection from endonucleases during repair (Figure 4H).

33Since MH is widely used in deletion events, we reasoned that we could program a library of targets to generate predictable mutations by introducing MH proximal to the cleavage site. With the same basic experimental scheme (Figure 1), we tested a library of 1000 targets and corresponding guides containing MH tracts of three different lengths (2, 4 or 6 bp) matching the sequence immediately upstream of the expected DSB site, and positioned 6 bp downstream of the cleavage site (Figure 5A-B). The resulting data were processed and analyzed similarly to the previous experiment.

34Intentionally programming MH tracts resulted in a high proportion of events corresponding to the expected deletions (8, 10 and 12 bp deletions for 2, 4 and 6 bp MH tracts, respectively; Figure 5B and C, Supplementary Figure S4). We also observe that the ratio of the programmed deletion increases as a function of length of the MH tract (Figure 5C). However, despite the greater predictability of which MH-mediated outcome would occur, the relative proportion of MH-mediated deletions increased only slightly from 76% to 82% (Figure 5D). Furthermore, we did not observe an excess of ‘imperfect’ MH-mediated events, e.g. an excess of 11 or 13 bp deletions in targets for which a 12 bp deletion was expected (Supplementary Figure S4). Nonetheless, the results show how targets that would result in diverse editing outcomes can be strongly biased towards a specific outcome by the presence of MH tracts (Figure 5E and F).

35The above results above suggest that the NHEJ-mediated repair outcomes for any given target sequence are both reproducible and dependent on sequence context. Accordingly, we next sought to train a machine learning model to predict these outcomes and their relative frequencies. We began by filtering out target sequences that were either poorly reproducible (low correlation between replicates, mainly due to low UMI counts; Supplementary Figure S5A) or poorly edited, resulting in a dataset of ∼1 million UMIs representing 4790 target sequences. On average, each target in this subset of the data used for modeling was represented by 204 UMIs and 28 alleles.

36Because larger events are rare in our data, we focused on predicting deletion events <30 bp in length, as well as all possible 1–2 bp insertion events at the DSB. Across all targets, we identified 557 ‘event classes’. The vast majority of CRISPR/NHEJ-mediated indels arising from any given target sequence should fall into one of these 557 event classes. We therefore framed our machine learning task as one of predicting, for an arbitrary target sequence, the relative frequency of CRISPR/NHEJ-mediated indels falling into each of these 557 event classes. These included 536 deletions (defined solely by their start/end points), all four possible single nucleotide insertions, all 16 possible dinucleotide insertions, and finally, a single event class for insertions greater than 2 bp in length. Of note, the 536 deletion event classes comprise almost all of the 550 possible combinations of start/end positions, with the constraints that deletions must be <30 bp and overlap with the –3/+2 window around the cleavage site. The 14 potential deletions that satisfy these constraints but were not observed in the modeling dataset were mainly large deletions.

37We also defined 3033 binary features to characterize the target sequence for which repair outcomes are being predicted. These are (i) Sequence features: 384 binary features corresponding to one-hot encoded sequence, including 80 for single nucleotide content (4 nucleotides × 20 positions) and 304 for dinucleotide content (16 dinucleotides × 19 positions); (ii) microhomology features: 2649 binary features corresponding to MH tracts; specifically, for each of the possible deletion event class, we defined 2–5 binary features (depending on the size of deletion) corresponding to the length of the MH tract ([0–4 bp × 519] + [0–3 bp × 7] + [0–2 bp × 6] + [0–1 bp × 4] deletion event classes = total of 2649 binary features) (Figure 6A).

38We split the 4790 target sequences in our modeling dataset into subsets of 3900 (for training), 450 (for validation) and 440 (for testing).

39Our model consists of three components: (i) predicting the ratio of insertions to deletions; (ii) predicting the distribution of 536 classes of deletion events; predicting the distribution of 21 classes of insertion events. The overall distribution of predicted outcomes is then determined by intersecting these three components. Of note, because we define deletion classes using the deletion start site and deletion length, two or more classes can effectively represent identical outcomes due to microhomology, but in a sequence specific manner (Supplementary Figure S5C). To address this, while evaluating performance, we simply collapsed identical outcomes.

40We trained predictors for each of the three components independently using logistic regression with a varied number of features and varied strength of L1 or L2 penalties. All of the models were trained on the training set using cross-entropy loss and evaluated on the validation set using the mean squared error (MSE). For the indel ratio predictor, we predicted that microhomology features and sequence context would both be important for prediction. However, including microhomology features did not improve the performance compared to using one-hot encoded sequence alone (MSE = 0.0203 and 0.0201 respectively, Supplementary Figure S5D). For the insertion predictor, it has been shown above that most insertions were templated by the sequence upstream of the cleavage site. We reasoned that sequence context around the cleavage site (±3 bp) should be sufficient to predict insertions. Consistent with this, including the full 20 bp target worsened performance, increasing MSE from 0.00666 to 0.00711 (Supplementary Figure S5E). For the deletion predictor, we compared the performance of models using sequence features only, microhomology features only or all features. The model with all-features performed the best (MSE of 0.000204, as compared with 0.000271 for sequence-only and 0.000208 for microhomology-only) (Supplementary Figure S5F). However, the nearly identical performance of the all-features versus microhomology-only models for predicting deletions is notable. We used the best performing model for each component to build a predictor for the overall distribution of outcomes (Figure 6A).

41Applying this model to the test set of 440 target sequences, which had been entirely held out from the training and validation steps, we compared the observed versus predicted frequencies of indels falling into various ‘event classes’. Observations and predictions were well matched for most targets, with a MSE of 0.000172 (Figure 6B). As a baseline, we also generated a set of predictions based simply on the aggregate frequencies of event classes in the training and validation datasets; as expected, these predictions performed more poorly (MSE of 0.000359; Figure 6B), confirming the improvement conferred by the model. Poorly predicted targets tended to be those with relatively shallower sampling of editing events, i.e. where our observed frequencies are noisier (Supplementary Figure S5B).

42While this manuscript was in preparation, several similar studies were published. Together with this manuscript, all four studies profiled repair outcomes of Cas9-induced DSBs at large numbers of endogenous or synthetic targets (Supplementary Table S1). The primary conclusions, e.g. that sequence context around the DSB, together with MH, are the major determinants of repair outcomes, are consistent between these studies as well as with earlier studies. In addition, Shen et al. built inDelphi and Allen et al. built ForeCasT, as models that predict NHEJ repair outcomes, analogous to the Lindel model described here. The ForeCasT model predicts deletions similarly to Lindel, as well as all possible 1–2 bp insertions (Supplementary Table S1). In contrast, the inDelphi model predicts the frequency of three classes of indels independently, using a neural network model for MH-mediated deletions (90 classes) and non-MH deletions (59 classes corresponding to 1–59 bp deletions, without prediction of location), and k-nearest neighbors model for 1 bp insertions (four classes) (Supplementary Table S1).

43We compared the three models by measuring the MSE on 440 targets in our test set as well as 4298 targets in ForeCasT test set (of note, the predicted probabilities of event classes not predicted by inDelphi were simply set to 0). Our model performed the best on our test set (MSE = 0.000172, 0.000225, 0.000212 for Lindel, ForeCasT and inDelphi, respectively), while ForeCasT performed the best on its test set (MSE = 0.000173, 0.000152, 0.000182, for Lindel, ForeCasT and inDelphi, respectively). We then combined our training set (3900 sequences) with ForeCast training set (10 725 sequences), resulting in a total of 14 625 sequences. We trained on this aggregated training set using the Lindel modeling approach, which resulted in the best overall performance (MSE = 0.000165 on our test set and 0.000125 on ForeCasT test set; Figure 6C, D). This final Lindel model is more accurate at predicting the ratio of insertions to deletions than ForeCasT (MSE = 0.01 and 0.02 for final Lindel model and ForeCasT, respectively; Supplementary Figure S6A,B). We further investigated the source of the errors for these models. Despite the fact that they implement different modeling approaches, both Lindel and ForeCasT’s mispredictions primarily lie with small deletions and 1 bp insertions (Supplementary Figure S6C–F).

44As a common use of CRISPR/Cas9 in conjunction with NHEJ is to introduce frameshifting mutations, we also assessed the observed versus predicted ratios of frameshifting indels for each of the 440 targets in our test set and 4298 targets in ForeCasT test set, and found them to be reasonably correlated (Pearson's r = 0.707, MSE = 0.0122 on our test set and Pearson's r = 0.676, MSE = 0.0098 on ForeCasT test set; Figure 6E; Supplementary Figure S6G). This result compares very favorably with the predictions of a previously published tool that we tested on this same task (Pearson's r = 0.283, MSE = 0.0431 on our test set and Pearson's r = 0.455, MSE = 0.0315 on the ForeCasT test set; Figure 6F; Supplementary Figure S6H).

45In summary, we developed an assay to systematically profile the diversity and relative frequencies of mutational events resulting from CRISPR/Cas9-mediated cleavage and NHEJ-mediated DSB repair of thousands of synthetic sequences. In applying this assay and analyzing the editing outcomes associated with 6,872 target sequences, we confirm that CRISPR/NHEJ-mediated repair outcomes for any given target sequence are reproducible, predictable, and largely shaped by the sequence context around the cleavage site.

46Our results also provide further insights into NHEJ-mediated repair of CRISPR/Cas9-mediated DSBs in human cell lines. First, we observe that insertion events are dominated by 1–2 bp insertions templated by the sequence immediately upstream of the cleavage site. Together with in vitro data from the literature, the data supports a model in which the sequence context around the DSB biases the extent to which cleavages are blunt-ended vs. include a 1–2 bp 5′ overhang. Such 5′ overhangs are repaired by gap-filling and ligation, resulting in asymmetrically templated 1–2 bp insertions. Second, we observe extensive usage of 1–4 bp microhomology in mediating deletion events, and furthermore show that repair outcomes can be strongly biased towards predictable outcomes by intentionally introducing MH tracts at specific distances from the DSB. Notably, however, the introduction of MH tracts did not substantially increase the proportion of MH-mediated events. Third, both MH and non-MH-mediated deletions were overwhelmingly unidirectional (i.e. extending either upstream or downstream from the DSB, rather than spanning it).

47Our assay has two main limitations. First, because of the locations of the PCR primer sites, we are only able to recover small deletions, and may be missing the rare, large deletion events that we and others have described. Greater knowledge of the frequency and determinants of large events is necessary to enable their prediction. Second, the lentiviral-based assay that we used fails to capture the influence of chromatin state on editing efficiency and repair outcomes. Lentivirus integrates to diverse locations across the genome, such that we are effectively observing an average, but integrations are biased towards open chromatin. As such, the patterns that we observe and model may be biased to this compartment. Furthermore, we are varying the sgRNA spacer sequence within the context of constant neighboring sequence, i.e. the lentiviral backbone. To the extent that this sequence biases nucleosome positioning, and that nucleosome positioning in turn influences Cas9 binding and cleavage, the patterns that we observe and model may be additionally biased. Additional systematic profiling of repair outcomes, in different compartments (e.g. open versus closed chromatin) and in different sequence contexts, will be necessary to understand the magnitude and nature of each of these potential sources of chromatin-mediated bias

48In addition to insights into into NHEJ-mediated repair of CRISPR/Cas9-mediated DSBs, our study also provides a new tool for sgRNA design for diverse goals. First, an important application of CRISPR/Cas9 is to achieve gene knockouts, a goal that depends on the efficient introduction of frameshifting indels. Dual cleavage with a variety of different nuclease systems has previously been shown to be an effective strategy for introducing frameshifting mutations. Our model's accuracy for predicting which sgRNAs/targets are likely to result in a high proportion of frameshifting indels will improve the viability of single cleavage with CRISPR/Cas9 for this same goal. Of note, our approach will not obviate the need for downstream validation to identify clones bearing the intended mutation, although it may reduce the number of clones that need to be screened. Second, for applications focused on mutation correction (e.g. using CRISPR/NHEJ to correct pathogenic mutations), the model may be useful for identifying sgRNAs/targets for which the desired outcome is predicted to occur at a high or sufficient frequency. Third, we and others have recently repurposed CRISPR/Cas9 as a tool for lineage tracing and/or molecular recording. For some goals (e.g. lineage tracing), the identification of ‘high entropy’ targets may critically enable the diversity necessary to uniquely label millions or billions of cells. For other goals (e.g. molecular recording), the design of ‘low entropy’ targets may facilitate predictable sequential editing. More generally, a deeper understanding of CRISPR/NHEJ-mediated mutations will strengthen our ability to precisely orchestrate not only the locations but also the outcomes of genome editing.

49Scripts and other software are available from our GitHub repository: https://github.com/shendurelab/Lindel.