The 15-mer sequence dataset used to train DeepSHM was generated as described in our previous publication (18) and is available at https://gitlab.com/maccarthyslab/deepshm. Germline IGHV reference sequences from the international ImMunoGeneTics information system (IMGT) (56) were downloaded and split into k-mers using a sliding window approach. Mutation frequencies for the central nucleotide in each k-mer were calculated by comparing against a B cell receptor (BCR) sequencing (BCR-seq) dataset from marginal zone, memory, and plasma B cells from healthy volunteers, as described in our previous study (57). To study intrinsic SHM patterns and avoid confounding issues arising from clonal selection in germinal centers, we used only non-productive sequences (containing internal frameshifts or stop codons) and clonally independent sequences (one sequence per clone, as assigned by Change-O (58), which uses CDR3s to segment clones) (57).

For this study, only 15-mers containing an AGCT motif, with either the G or the C as the central nucleotide, were used with DeepSHM to predict mutation frequencies. The total 15-mer dataset was processed to extract those containing AGCTs using custom Python scripts. Our statistical analysis of synonymous mutations (those that do not change the protein sequence) ablating the CAGCTG motif was done with productive, clonally independent sequences. All statistical tests were performed using SciPy (59).

DeepSHM (https://gitlab.com/maccarthyslab/deepshm) is a deep learning model that uses a convolutional neural network architecture to predict mutation frequency or substitution rate of the central nucleotide in a k-mer of size 5, 9, 15 or 21. We used the 15-mer mutation frequency model, which takes a DNA sequence of 15 nucleotides as input and outputs a predicted mutation frequency value between 0 and 1. The 15-mer sequences were encoded into a 4 x 15 binary matrix, with rows corresponding to the 4 nucleotides and columns to the 15 positions along the k-mer. In each column, a 1 was placed in the appropriate row to denote the base identity for that position, while the remaining rows were 0s. This procedure, called one-hot encoding, is a common method for converting categorical data (A, G, C, T) into a machine-readable format (0s and 1s).

We downloaded the h5 file containing the model (model_15_mf.h5) and used it with Python to predict mutation frequencies for AGCT 15-mers in our dataset.

Integrated gradients is an attribution method that measures the impact of individual inputs towards the output prediction of a deep learning model (60). It relies on a baseline value to compute a path integral of the model’s gradients with respect to its inputs, from the baseline to the input value. Since our input data is binary, we used a zero-matrix as our baseline with 50 steps taken from baseline to input. All 15 nucleotides in the 15-mer are considered as input features in the prediction of the mutation frequency of the central nucleotide, hence integrated gradients calculates a score for each base according to its impact on the output prediction.

The following GitHub repository was used to compute the integrated gradients scores for each of the DeepSHM predictions: https://github.com/hiranumn/IntegratedGradients. The repository was cloned and imported into the python script where the DeepSHM predictions were being run and used to compute integrated gradients scores for each input in each prediction. We generated sequence logo plots to visualize the frequency of nucleotides occurring at each position across 15-mers in each subregion using Logomaker (61), which is available at the following GitHub repository: https://github.com/jbkinney/logomaker.

MOODS (https://github.com/jhkorhonen/MOODS) is a position-weight matrix (PWM) matching algorithm that takes sequences and a counts matrix as inputs and outputs match scores for a segment of the sequences (62). The counts matrix is a 4×n matrix where the rows correspond to nucleotides (A, G, C, or T) and the columns correspond to positions along the TF binding motif, with the number of counts for each nucleotide in each position empirically obtained using SELEX and available on the JASPAR database (63). MOODS uses log-likelihood scoring to convert the counts matrix to a PWM, which it then compares against the sequence to generate a match score, only reporting scores at positions that exceed a P value cutoff of 0.001. We used MOODS to gauge the fidelity of our 15-mer sequences to binding motifs for the E-box TFs, E2A (https://jaspar2020.genereg.net/matrix/MA0522.2/) and TFAP4 (https://testjaspar.uio.no/matrix/MA1570.1/).

The E2A ChIP-seq data in Ramos cells was taken from a previous study (45) and is available at https://www.ncbi.nlm.nih.gov/bioproject/587064. The E2A ChIP-seq data in GM12878 cells was taken from ENCODE (64) and is available at https://www.encodeproject.org/experiments/ENCSR000BQT/. Both datasets were subject to the same analysis pipeline - a local alignment to hg38 of both case and control datasets using bowtie2 (65), sorting and indexing the resulting bam file using samtools (66). The callpeaks function of MACS2 (67) was run on the aligned bam files using the default q-value cutoff of 0.05 to call peaks. These peaks were then subject to motif enrichment analysis using the findMotifsGenome.pl function of the HOMER suite (68), using the hg38 genome and default size setting of 200 bp. The RPKM calculation was conducted using deeptools bamCoverage (69) with the bin size parameter set to 500 bp.

V regions can be structurally divided into antigen-binding complementarity-determining regions (CDR1–3) and intervening structural framework regions (FW1–3). AGCT is one of the most highly mutated WRCH motifs in V regions (70). The fact that AGCT is palindromic increases the probability that AID will deaminate cytosines on both the forward and the reverse strands (70). All human IGHV genes (at least the *01 IMGT alleles), except for three from the IGHV2 family, have one or more AGCT motifs near the 5’ end located in FW1 (57).

Our previous analysis of IGHV3–23*01 non-productive sequences showed higher mutability of AGCTs in the CDRs and a particularly low mutability of the 5’ AGCT in FW1 (14). To examine this differential mutability of AGCT motifs across all IGHV genes, we used DeepSHM to predict mutation frequencies of the central nucleotides in the AGCT 15-mers in our dataset and compared the results with the observed data using a correlation analysis. We achieved a Pearson correlation of r=0.92 for those with a central G site ( Figure 1A ) and r=0.86 for those with a central C ( Figure 1B ). These high correlations suggested that the neural network had identified sequence features that distinguish low from high mutation frequencies for AGCT sites. To confirm that the positional differences in AGCT mutability previously observed for IGHV3–23*01 applied to other IGHV alleles, we examined our dataset of 16,870 15-mers across 65 IGHV alleles and their associated mutation frequencies (57). We separated the AGCT 15-mers in our dataset by IMGT subregion and plotted the observed and predicted mutation frequencies for those with central Cs ( Figure 1C ) and central Gs ( Figure 1D ). We observed a statistically significant difference in mutation frequency between the AGCT motifs in FW1 and all other subregions for both observed and predicted datasets centered on the G (t-test, P<10^–30) ( Figure 1C ) and C (t-test, P<10^–20) ( Figure 1D ).

We conclude that that AGCTs in FW1 are significantly less mutated than those in other V subregions and that DeepSHM can recapitulate these observed positional differences in AGCT mutability.

To interrogate the specific sequence features associated with high or low mutation frequency predictions, we used an interpretability method, integrated gradients (22). Integrated gradients analyses involve computing the derivative of the output (mutation frequency prediction) with respect to the input (15-mer sequence) to ascribe importance to input features based on their impact on the output prediction. Specifically, a higher integrated gradients score would imply that a small change in input had a more positive contribution towards the output prediction. Conversely, lower integrated gradients scores indicate that changes in input features contributed negatively to the predicted output.

We generated integrated gradients scores for each nucleotide within the AGCT 15-mers for its prediction of (i.e. contribution towards) the mutation frequency of the central G or C within the hotspot. We plotted the range of scores for each position as boxplots which were further categorized based on the location of the 15-mers within FW and CDR subregions ( Figures 2A–J ).

As a positive control, we observed that nucleotides within the AGCT hotspot across all central G 15-mers ( Figures 2A–E ) and 85% of central C 15-mers ( Figures 2F–J ) had a positive integrated gradients score, meaning that the presence of these nucleotides increased the mutability prediction of the central G or C. We then examined the integrated gradients scores for FW1 AGCTs ( Figures 2A, F ) as they are significantly less mutated than AGCTs in other V sub-regions ( Figure 1C ). We found that 78% of the lowly mutating FW1 AGCTs were flanked by a 5’-C and 3’-G nucleotide, both of which have large negative integrated gradient scores ( Figure 2F ). This suggests that an extended CAGCTG motif context decreases the mutability prediction for the central G and C nucleotides within these lowly mutating FW1 AGCTs, implying that CAGCTG motifs may be less frequently targeted by AID.

To explore this idea further, we directly compared integrated gradients scores of the 5’ and 3’ nucleotides flanking AGCTs across all 15-mers. We found that for CAGCTG-containing 15-mers, the integrated gradients scores were almost always negative for the 5’-C (98%) and consistently negative for the 3’-G (100%), supporting the notion that the CAGCTG context has a predominantly negative influence on mutagenesis of AGCT ( Figure 3 ). Additionally, of all AGCT flanking nucleotide combinations, the 5’-C and 3’-G combinations were overwhelmingly within the FW1 region and were significantly less mutated than AGCTs with other flanking nucleotide combinations (Mann-Whitney test, P<10^–30) ( Figure 4 ). We note that although most of the FW1 AGCT motifs were flanked by 5’-C and 3’-G nucleotides, even those flanked by other nucleotide combinations tended to have lower mutation frequencies (blue dots in Figure 4 ). This suggests that the position of the AGCT within the V region may also have some influence on its mutability. In addition, a fraction of CAGCTG motifs in FW2 undergo higher mutation frequency than those in FW1 ( Figure 4 ), suggesting the presence of additional mechanisms, possibly involving differences in the larger sequence context of these motifs, that influence differential mutability, which we address in the following section.

We conclude that the weakly mutated AGCTs in FW1 are predominantly flanked by 5’C and 3’G nucleotides, implying that the CAGCTG sequence context correlates with reduced SHM of AGCT motifs.

CAGCTG corresponds to the CANNTG E-box binding motif of the basic helix-loop-helix TF family, which includes E2A (50). To predict binding probabilities of E2A to AGCT 15-mers, we used MOODS, a TF-binding prediction package which utilizes counts matrices obtained from empirical SELEX data, wherein higher MOODS scores reflect a stronger sequence match to a particular TF binding motif (62).

We found significant negative Pearson correlations of r=-0.64 and r=-0.65 between the MOODS binding scores for E2A motifs and the mutation frequencies of the central Gs ( Figure 5A ) and central Cs ( Figure 5B ) in the AGCT sites, respectively. Similar analysis for TFAP4, another E-box TF commonly expressed in B cells, showed weaker correlations ( Supplementary Figure 1 ). The predicted MOODS scores for the AGCT 15-mers fell roughly into three discrete tiers. Tier 1, having the highest MOODS scores but generally lower mutation frequencies, and consisting almost entirely of CAGCTG-containing 15-mers ( Figures 5A, B ). Tier 2, having intermediate MOODS scores with a wide range of mutation frequencies. Importantly, although this tier consists of a mixture of CAGCTG and non-CAGCTG 15-mers, the former showed a tendency to be less mutated than the latter ( Figures 5A, B ). Interestingly, the FW2 CAGCTG 15-mers observed to be highly mutating in Figure 4 fall into this tier and contain central Gs (green crosses in Figure 5A ). This indicates lower fidelity to the E2A motif than the lowly mutating FW1 CAGCTGs and may explain, in part, the higher mutation frequency due to diminished E2A binding. Tier 3, which harbored the lowest MOODS scores and generally higher mutation frequencies, consisted mostly of non-CAGCTG 15-mers ( Figures 5A, B ). These results suggest that the binding probability of E2A to an AGCT-centered 15-mer negatively correlates with mutation frequency of that AGCT.

Next, we determined the distribution of all potential E2A sites across all human IGHV genes. We segmented germline sequences for the 220 alleles obtained from the IMGT database into six subregions and counted the number of occurrences of CAGCTG ( Figure 6A ) and the more general CANNTG ( Figure 6B ) E-box motifs. CAGCTG motifs were mostly distributed in the FW1 region, with the IGHV2 family notably lacking them ( Figure 6A ). The IGHV4 family has the highest density of CANNTG motifs while most of IGHV2 family have E-box motifs in the FW3 region ( Figure 6B ). Overall, each of the 220 alleles had at least one E-box motif, with most of them in FW1 and/or the leader-intron-leader (L-intron-L) sequence which immediately precedes FW1 ( Figure 6B ).

To determine whether E2A associates with IGHV regions, we analyzed E2A ChIP-seq datasets derived from Ramos (45) and GM12878 (64) B cell lines. After aligning these data to the hg38 reference genome using bowtie2 (65), we used MACS2 (67) to call peaks and then conducted a motif enrichment analysis using HOMER (68). We saw that for both Ramos ( Supplementary Table 1 ) and GM12878 ( Supplementary Table 2 ) cell lines, the CAGCTG motif corresponding to the E2A TF binding motif was highly enriched and among the top two most significant results. We calculated the reads per kilobase million (RPKM) values of 500bp bins in both the E2A ChIP-seq and the IgG control ChIP-seq alignments and plotted their correlations. In both Ramos ( Figure 7A ) and GM12878 ( Figure 7B ) cells, the bins containing the rearranged V region (IGHV4–34 in Ramos (71) and IGHV3–21 in GM12878) (72)) were enriched for E2A binding, as was the bin containing the IGH Eμ enhancer, which serves as a positive control for E2A binding (43). However, bins containing a negative control region, TRBV20–1, a commonly used T-cell receptor V gene (73), showed no enrichment for E2A binding in either cell line ( Figures 7A, B ). Thus, E2A can directly associate with V regions.

Altogether, these results suggest that E2A association with CAGCTG motifs in FW1 suppress AID targeting to these AGCTs, thereby providing a plausible explanation for the strong negative correlation between this motif context and mutation frequency.

We additionally sought to compare the role of E2A binding with another sequence-level determinant of AGCT mutability proposed in a recent study (19), namely, the presence of PyPy dimers. In this study, a higher frequency of PyPy dimers in the 6 nucleotides 5’ to the AGCT was associated with increased mutation frequency of the central C residue (19). Therefore, we counted the number of PyPy dimers in the 6 nt region immediately upstream of AGCT motifs in our k-mer dataset and fit a linear regression model predicting mutation frequency, including an indicator variable for the presence of the E-box motif. This binary E-box indicator variable had a Pearson correlation of r=-0.61 with the mutation frequency, while the integer PyPy count had a correlation of 0.51 ( Table 1 ). Importantly, therefore, both variables individually correlate with mutation frequency in directions consistent with our expectations, that is, positive for PyPy counts, which increases mutation frequency, and negative for the presence of an E-box, which reduces mutation frequency.

Given this trend, we expected a linear model with both variables to have a higher performance than a model with either individual variable, as measured by R², which directly reflects the proportion of variance in the output variable (mutation frequency) explained by the input variables (E-box motif, PyPy richness, or both). The combined regression model achieved an R² of 0.45 meaning that 45% of the variance in mutation frequency is explained by the presence of E-box and PyPy motifs ( Table 1 ). The regression model with only the PyPy counts variable achieved an R² of 0.25 while the model with only the E-box indicator variable achieved a higher, and closer to the combined, R² value of 0.39 ( Table 1 ). Of note, the coefficients generated by the model had signs appropriate to the direction of correlation with mutation frequency, that is, positive for PyPy counts and negative for the E-box indicator ( Table 1 ).

These results lead to the conclusion that both mechanisms, decreasing mutation frequency of FW1 AGCTs, plausibly through E2A binding, and increasing mutation frequency of AGCTs through increased AID binding to flexible PyPy-rich DNA can contribute independently to the observed mutability. Importantly, however, the R² values observed from these analyses also imply that additional mechanisms are necessary to fully explain AGCT mutability.

To better understand the relationship between the E-box motif and the mutability of the central nucleotides, we examined mutations of the CAGCTG hotspot in FW1. We hypothesized that if this motif context negatively contributes to SHM, then naturally occurring mutations that ablate this context would be expected to increase mutation frequency of the AGCT within it.

Due to the paucity of non-productive sequences in our dataset, we examined synonymous mutations (i.e. those that do not cause changes in protein sequence) in productive BCR sequences to preclude any effect of affinity selection. Specifically, we focused on the G residues at positions 3 (G₃) and 6 (G₆) of CAGCTG. Importantly, the FW1 CAGCTG motif occurs at position 7–12 of the V segment and is always in frame, such that mutations at G₃ and G₆ are in the third position of their respective codons. Thus, G₃>A₃ mutations are synonymous since CAG and CAA are degenerate codons for glutamine. Similarly, G₆>H₆ mutations (where H = A/C/T) are also synonymous since CTG, CTA, CTC and CTT are degenerate codons for leucine. Importantly, G₆>H₆ mutations (CAGCTH) would ablate the E-box motif. Thus, we compared mutation frequency at the central AGCT in clonal groups having an unmutated CAGCTG in FW1 or a CAGCTH in the same position.

To prevent double counting of mutations occurring during clonal expansion, we selected a sequence at random from each clonal group (57). Our sequence data consisted of 642,367 clonal groups, of which 504,333 had a sequence identity in position 7–12 of the V segment, corresponding to one of the four motifs of interest: unmutated CAGCTG (G₃/G₆), single mutant CAACTG (A₃/G₆), single mutant CAGCTH (G₃/H₆) and double mutant CAACTH (A₃/H₆). We counted the number of occurrences of each ( Table 2 ) and used these numbers to calculate mutation frequencies for sites 3 and 6 ( Supplementary Table 3 ). We observed that the mutation frequency of G₃ increases from 9.6% when G₆ is unmutated to 21.5% when G₆ is mutated, a highly significant difference (Fisher test, P<10^–16) ( Figure 8 ). Thus, the presence of an intact CAGCTG motif is strongly associated with lower mutation of the AGCT within it.

Altogether, our results support the notion that FW1 AGCT motifs occurring in the context of the E-box binding motif, CAGCTG, lead to dampened SHM in these locales.