A comprehensive approach for genome-wide efficiency profiling of DNA modifying enzymes

Abstract

A precise understanding of DNA methylation dynamics is of great importance for a variety of biological processes including cellular reprogramming and differentiation. To date, complex integration of multiple and distinct genome-wide datasets is required to realize this task. We present GwEEP (genome-wide epigenetic efficiency profiling) a versatile approach to infer dynamic efficiencies of DNA modifying enzymes. GwEEP relies on genome-wide hairpin datasets, which are translated by a hidden Markov model into quantitative enzyme efficiencies with reported confidence around the estimates. GwEEP predicts de novo and maintenance methylation efficiencies of Dnmts and furthermore the hydroxylation efficiency of Tets. Its design also allows capturing further oxidation processes given available data. We show that GwEEP predicts accurately the epigenetic changes of ESCs following a Serum-to-2i shift and applied to Tet TKO cells confirms the hypothesized mutual interference between Dnmts and Tets.

Keywords: 5hmC; 5mC; DNA demethylation; DNA methylation; Dnmts; Tets; epigenetics; hairpin sequencing; hidden Markov model.

Graphical abstract

Figure 1

GwEEP - Pipeline overview (A) Laboratory pipeline: (1) Genomic DNA is digested by endo nucleases followed by (2) Klenow exo-catalyzed A-tailing. (3) A-Tailed DNA molecules are subjected to sequencing adapter, hairpin linker ligation and (4) subsequent enrichment of hairpin-ligated molecules. (5) Half of the library is used for BS, the other half for oxBS treatment. (6) After amplification and indexing using PCR the libraries are sequenced on an Illumina platform with minimum 100 bp in a paired-end mode (created with BioRender.com). (B) Computational processing: Illumina raw data are processed into base calls (FASTQ) and trimmed for adapter and hairpin linker sequences. Bisulfite reads from the same molecules are paired to restore the genomic sequence for efficient mapping. Subsequently, the double-strand information is annotated and stored in DSI (double strand information) files. The HMM then derives 5mC and 5hmC distributions, as well as the efficiencies of Dnmt and Tets which are stored in the IGV file format. Both DSI and IGV files can be visualized using the IGV genome browser. (created with BioRender.com).

Figure 2

RRoxBS and HMM results (A) Average CpG methylation level based on uncorrected hairpin sequencing counts for WT ESCs. (B) Average nonCpG methylation level based on uncorrected hairpin sequencing counts for WT ESCs. (C) Demethylation rate in WT and Tet TKO cells. (D) Average CpG methylation level based on uncorrected hairpin sequencing counts for Tet TKO ESCs. (E) Average nonCpG methylation level based on uncorrected hairpin sequencing counts for TKO ESCs. (F) Relative difference in demethylation rate between WT and Tet TKO cells. (G) Estimated, and conversion error corrected, 5mC and 5hmC distribution after HMM's application. (H) HMM-derived maintenance methylation, de novo methylation and hydroxylation efficiencies. (I) Integrative Genomics Viewer (IGV) snapshot across a gemoic region located at the Tet3 gene showing the distribution of CpGs across CpG-rich and -poor regions of RRHPoxBS.

Figure 3

Model and computational methods (A) Transitions between methylation states of a single CpG dyad; u indicates an unmethylated, m a methylated, and h a hydroxylated state of a CpG. ${\mu }_d$ describes the efficiency of de novo methylation, ${\mu }_m$ the efficiency of maintenance methylation and η the hydroxylation efficiency. λ represents the overall methylation efficiency (maintenance $+$de novo) that is defined as $\lambda ={\mu }_m+{\mu }_d-{\mu }_m\cdot {\mu }_d$. The parameter p describes the probability that maintenance methylation does not consider hemihydroxylated sites. (B) Possible conversion errors during bisulfite and oxidative bisulfite sequencing. (C) Metropolis-Hastings update step: Assuming each efficiency is a linear time function, each next value is sampled using two truncated normal distributions in two consecutive steps. Step 1: ample the intercept ${\mathbf{y}}_{i-1}$ from the truncated normal with mean ${\mathbf{x}}_{i-1}$ and bounds [0, 1]. Step 2: sample the gradient ${\mathbf{y}}_i$ from the truncated normal distribution with mean ${\mathbf{x}}_i$ and bounds $[a_i,b_i]$, which depend on the sampled intercept ${\mathbf{y}}_{i-1}$ of Step 1. (D) Clustering of estimated enzymatic efficiency with intercept ${\beta }_0$ and gradient ${\beta }_1$ for CpGs A, B, C, D using k-means versus k-error algorithm. (E) Spatial auto- and cross-correlations of maintenance methylation, de novo methylation and hydroxylation efficiencies over the whole genome.

Figure 4

DNA methylation segmentation (A) Percentages of HMRs, PMDs, LMRs, and UMRs across the genome based on WBGS data derived from WT ESC by Ficz et al. (2013). (B) HMM estimated 5hmC distribution across the distinct segments based on RRHPoxBS of WT ESCs. (C) Size distribution of the individual segment types based on WBGS data derived from WT ESCs by Ficz et al. (2013). (D) HMM estimated methylation distribution of HMRs, PMDs, LMRs, and UMRs (methylation levels of the individual segment types for WT and Tet TKO ESCs). (E) Methylation level of HMRs, PMDs, LMRs, and UMRs based on WBGS data derived from WT ESC by Ficz et al. (2013). (F) Estimated HMM enzyme efficiencies in HMRs, PMDs, LMRs, and UMRs for WT and Tet TKO ESCs based on RRHPoxBS.

Figure 5

Enzymatic profiles across genes and TFBSs (A) Comparison between ChIP profiles of epigenetic writers and estimated enzymatic efficiencies in WT ESCs at day 0 (serum/LIF) across expressed and low/non-expressed genes. (B) Estimated enzymatic efficiencies in WT and Tet TKO ESCs across genes. (C) Estimated enzymatic efficiencies of WT and Tet TKO ESCs at Sox2-, Tet1-, and H3K4me3-enriched regions obtained from ENCODE.