Probing DNA shape and methylation state on a genomic scale with DNase I

Abstract

DNA binding proteins find their cognate sequences within genomic DNA through recognition of specific chemical and structural features. Here we demonstrate that high-resolution DNase I cleavage profiles can provide detailed information about the shape and chemical modification status of genomic DNA. Analyzing millions of DNA backbone hydrolysis events on naked genomic DNA, we show that the intrinsic rate of cleavage by DNase I closely tracks the width of the minor groove. Integration of these DNase I cleavage data with bisulfite sequencing data for the same cell type's genome reveals that cleavage directly adjacent to cytosine-phosphate-guanine (CpG) dinucleotides is enhanced at least eightfold by cytosine methylation. This phenomenon we show to be attributable to methylation-induced narrowing of the minor groove. Furthermore, we demonstrate that it enables simultaneous mapping of DNase I hypersensitivity and regional DNA methylation levels using dense in vivo cleavage data. Taken together, our results suggest a general mechanism by which CpG methylation can modulate protein-DNA interaction strength via the remodeling of DNA shape.

Fig. 1.

Deep sequencing reveals striking positional dependencies between nucleotide positions within the DNase I recognition site. (A) Position-specific relative cleavage rate parameters as derived from DNase I digestion of human genomic DNA (normokaryotypic IMR90 fibroblasts) under the assumption of independence between nucleotides. Dependence on local sequence context is largely limited to a hexamer centered at the cleaved backbone bond. (B) Comparison between cleavage rates for pairs of hexamers that are related by a single-nucleotide substitution. The slope of the dashed line corresponds to the position-specific cleavage rate in panel A, and is directly related to the “unconditional” ΔΔG, the change in binding free energy associated with the point mutation. The fold change in cleavage rate due to a mutation from G to T at position −1 is largely independent of the base identity of the five neighboring nucleotides. (C) Breakdown of the independence assumption (dashed line). The effect on cleavage rate of a point mutation from A to C at position +2 is highly dependent on the base identity at the “modulating” position +1. Using a “conditional” ΔΔG for each possible base at position +1 (colored lines) provides a far more accurate description. (D) The strength of the positional dependencies can be quantified in terms of a new quantity “ΔΔΔG,” defined as the difference between the conditional and unconditional ΔΔG. The values in the highlighted row and columns correspond to the ratio in slope between each of the colored solid lines and the dashed line in C. Far away from the diagonal, ΔΔΔG becomes numerically small (white in heat map), indicating an increasing degree of independence.

Fig. 2.

MGW is predictive of DNase I cleavage rate. (A) ΔΔG derived from the negative logarithm of cleavage rate as a function of MGW at the six positions of all 4,096 unique hexamers. MGW of this region was predicted for naked binding sites based on a pentamer-based HT shape prediction approach (SI Methods). HT predictions for all possible 16 dinucleotide flanks were averaged and values of MGW that fall within intervals of 0.3 Å assigned to groups of sequences for which cleavage rates are shown as box plots. (B) DNase I–DNA complex based on crystal structure (Protein Data Bank ID code 2DNJ). Base pairs at positions −3 and −2, where DNase I cleavage anticorrelates with MGW, are highlighted in blue. Base pairs at positions −1 and +1, where DNase I cleavage correlates positively with MGW, are highlighted in green. Regions where no correlation could be detected are shown in gray. The color code of the base pairs in the crystal structure is equivalent to the one used for the box plots. (C) DNase I–minor groove contacts within a distance of 5 Å from any base atom are shown for the same crystal structure. Arg41 and Arg9 bind upstream of the cleavage site, where MGW anticorrelates with DNase I cleavage (blue base pairs). This anticorrelation likely arises from the attraction between the positively charged arginine residues and the locally enhanced negative electrostatic potential. The cleavage site (indicated by the orange arrow), by contrast, is located in a region where MGW correlates positively with DNase I cleavage (green base pairs).

Fig. 3.

Observation and analysis of the effect of methylation on DNase I cleavage rate. (A) The rate of cleavage depends strongly on the DNA methylation status. We used a positional map of DNA methylation in IMR90 (36) to delineate subsets of genomic positions with low/high degrees of CpG methylation, respectively. Comparison between the hexamer cleavage rates derived from these respective subsets shows an ∼eightfold increase in cleavage rate for hexamers with a methylated CpG immediately downstream of the cleaved phosphate (red points). (B) Interplay between DNA sequence and methylation status, DNA geometry, and DNase I cleavage suggested by our analysis. (C) Roll and MGW of methylated and unmethylated versions of the same hexamer based on the average of MC predictions for three different flanking sequences (see SI Methods for details). Methylation leads to an increase in the positive roll angle at the CpG dinucleotide and a narrowing of the MGW at position −2 by roughly 0.5 Å. (D) The effect of methylation on DNase I cleavage can be predicted in silico by training a model to predict the cleavage rates of unmethylated DNA sequences of type NNNCGN using information on DNA MGW and roll angle along these same unmethylated sequences. An increase in cleavage rate (i.e., data points shifting downward) is predicted when MGWs and roll angles for the methylated versions of the sequences are supplied as input to the model.

Fig. 4.

Genomic methylation status can be predicted from dense in vivo DNase I footprints. Starting from in vivo DNase I footprinting data for the IMR90 cell line, a set of nonoverlapping windows (2,500 bp long) containing at least 400 cleavage events with hexamer context NNNpCGN was identified. Next, for each window, the observed number of cleavages upstream of CpG dinucleotides was compared with the expected number. This allowed us to infer the methylation status of the corresponding DNA. To validate our predictions, we ranked all windows by their actual degree of methylation as measured by Lister et al. (36). Shown is the distribution of ranks for the subset of windows predicted to be hypomethylated (red) and hypermethylated (blue), respectively.