Comparison of sequencing-based methods to profile DNA methylation and identification of monoallelic epigenetic modifications

Abstract

Analysis of DNA methylation patterns relies increasingly on sequencing-based profiling methods. The four most frequently used sequencing-based technologies are the bisulfite-based methods MethylC-seq and reduced representation bisulfite sequencing (RRBS), and the enrichment-based techniques methylated DNA immunoprecipitation sequencing (MeDIP-seq) and methylated DNA binding domain sequencing (MBD-seq). We applied all four methods to biological replicates of human embryonic stem cells to assess their genome-wide CpG coverage, resolution, cost, concordance and the influence of CpG density and genomic context. The methylation levels assessed by the two bisulfite methods were concordant (their difference did not exceed a given threshold) for 82% for CpGs and 99% of the non-CpG cytosines. Using binary methylation calls, the two enrichment methods were 99% concordant and regions assessed by all four methods were 97% concordant. We combined MeDIP-seq with methylation-sensitive restriction enzyme (MRE-seq) sequencing for comprehensive methylome coverage at lower cost. This, along with RNA-seq and ChIP-seq of the ES cells enabled us to detect regions with allele-specific epigenetic states, identifying most known imprinted regions and new loci with monoallelic epigenetic marks and monoallelic expression.

Figure 1. Critical parameters in sequencing-based DNA methylation profiling

Sequencing statistics and CpG coverage are shown for MethylC-seq, (207 lanes, data analyzed here were from Lister, et al. .), RRBS (2 lanes), MeDIP-seq (4 lanes each), MBD-seq (3 lanes), and MRE-seq (3 lanes each). Since the amount of sequence produced per lane is increasing, we also provide “Gbp of sequence” as a measure of the relative cost of each method. The methods differ significantly in total bases generated by the Illumina sequencer, total high quality bases passing Illumina chastity filtering and mapping uniquely, and total bases used for generating methylome maps (high quality bases passing redundancy filters) (a) The H1 replicates assayed and the Gbp of sequence at successive processing stages by each method are shown. The bisulfite-based methods and MRE-seq resolve the methylation status of individual cytosines, whereas the MeDIP-seq and MBD-seq read mappings are extended to 150bp, resulting in a maximum resolution of 150bp. This extension is applied to calculations of CpG coverage but is not applied to the Gbp of sequence at the processing stages. Coverage information is shown for repeats (primarily transposon sequences) genome-wide. While maximum resolution of each method is reported, resolution can be assessed at various levels. As the level of resolution decreases, as a consequence of averaging of methylation scores over a window of larger size for example, imperfect coverage and limited accuracy become less limiting, provided that the average score is not affected by systematic biases in coverage and accuracy. Thus, methylome coverage and accuracy in methylation calls are a function of resolution. (b) The percentage of genome-wide CpGs (28,163,863) covered by multiple, single or no methods are shown. The percentage of CpGs covered genome-wide (c) or in CGI (d) are plotted as a function of read coverage threshold.

Figure 2. Comparison of bisulfite-based methods

(a) Calls of highly/partially/weakly methylated (0.80-0.20 or 0.75-0.25 cutoff) or highly/weakly methylated (0.20 cutoff) were made for CpGs covered at several minimum read depths by MethylC-seq and by RRBS (both on replicate #3). The number and percent of genome-wide CpGs covered and the percent of concordant calls are shown for each minimum read depth and methylation call cutoff. (b) Differences (MethylC-seq - RRBS) in methylated proportions (methylated reads/(methylated reads + unmethylated reads)) for CpGs with a minimum coverage of 5 reads by both methods. Percentages of concordant and discordant methylation were determined at cutoffs of + and − 0.1 (green dashed lines) and 0.25 (red dashed lines). (c) CpG density in a 400bp window and (d) genomic context of concordant and discordant CpGs at the 0.25 cutoff.

Figure 3. Comparison of methylated DNA enrichment methods

(a) Calls of highly/weakly methylated were made by averaging methylation scores for CpGs covered at varying minimum read depths by MeDIP-seq or MBD-seq in 1000bp and 200bp windows. The number of windows, percent of genome-wide CpGs covered, and the percent of concordant calls are shown for each minimum read depth and window size. For the 1000bp windows with a minimum read depth of 5, the (b) CpG density and (c) genomic context of the concordant and discordant windows are shown. The inset in (b) shows a close-up of the concordance/discordance of CpG densities consistent with CGIs. (d) For the 1000bp windows with a minimum read depth of 5, MethylC-seq methylation proportions for CpGs and non-CpG cytosines covered at a minimum read depth of 5, 444,590 windows, were summed and the windows were binned by the sum. For each of these bins, the number of windows called highly methylated by MeDIP-seq or MBD-seq is shown on the left y-axis and the percent of total windows with calls of highly methylated is shown on the right y-axis. Windows with a MethylC-seq methylation proportion sum greater than 15, representing 83% of all windows, were called highly methylated by MeDIP-seq and MBD-seq in 99.9% of cases. The windows with a methylation proportion sum from 1–15, representing 17% of all windows, were called highly methylated by MeDIP-seq and MBD-seq in at least 99.1% of cases.

Figure 4. Comparison of all methods

(a) The table shows the percentage of 1000bp windows with concordant and discordant MethylC-seq (replicate #3), RRBS (replicate #3), MeDIP-seq (replicate #2) and MBD-seq (replicate #2) calls at minimum read depths of 5 and 10. Methods making the same call are grouped together in parentheses. Calls were made for MethylC-seq and RRBS by averaging the methylation proportion of CpGs within the window that were covered at the minimum read depth and applying a highly/weakly methylated cutoff of 0.2. Calls were made for MeDIP-seq and MBD-seq by averaging the methylation score of CpGs within the window that were covered at the minimum read depth. (b) Genome browser view of the 100kb CpG rich Protocadherin alpha cluster (PCDHA) exemplifying the significant concordance in methylation status seen on a genome-wide level. For Methyl-seq and RRBS, the Y- axis displays methylation scores of individual CpGs. Scores range between -500 (unmethylated) and 500 (methylated), and the zero line is equivalent to 50% methylated. Negative scores are displayed as green bars and positive scores are displayed as orange bars. For MeDIP-seq 1, MeDIP-seq 2 and MBD-seq, the Y-axis indicates extended read density. Browsable genome-wide views of these datasets are available at http://www.genboree.org and http://genome.ucsc.edu/.

Figure 5. Integrative method increases methylome coverage and enables identification of a DMR

(a) MRE-seq involves parallel digests with methylation sensitive restriction enzymes (HpaII, AciI, and Hin6I), selection of cut fragments of approximately 50bp–300bp, pooling the digests, library construction, and sequencing. For every 600bp window along chromosome 21, MeDIP-seq scores were plotted against MRE-seq scores. The plot depicts the inverse relationship between MRE-seq and Me-DIP-seq signals. (b) Coverage of CpGs in the human genome by MeDIP-seq alone (red), MRE-seq alone (green), both (yellow), or neither method (no fill). Sequence from replicate #1 and #2 were used in these calculations (c) UCSC Genome Browser view of ZNF331in H1 ESC, showing overlap of Me-DIP-seq, MRE-seq and H3K4me3 (from ChIP-seq) signals at bisulfite region 1 and only MeDIP-seq signal at bisulfite region 2. (d) Clonal bisulfite sequencing results for specified regions in ESC from replicate #1. A filled circle represents a methylated CpG and an open circle indicates an unmethylated CpG.

Figure 6. Allelic DNA methylation, histone methylation, and gene expression in embryonic stem cells

(a) Venn diagram summarizing the number of loci exhibiting monoallelic DNA methylation, histone methylation or monoallelic expression, and their overlap. Allelic methylation or expression was determined by comparison of genetic variation between pairs of assays. To further evaluate the top 1000 loci (average size of 2.9kb and encompassing a CGI) with potential allelic DNA methylation, the following pairs of assays were used: MRE-Seq and MeDIP-Seq for allelic DNA methylation within the loci, MethylC-seq and expression data for monoallelic expression of genes associated (+/−50kb) with the loci, MethylC-seq and histone modifications H3K4me3 and H3K9me3 for monoallelic histone methylation within 1kb from the loci. Evidence of genomic variation was found for 119 of 1000 loci; 4 loci in 3 genomic locations provided confirmation of monoallelic status from all three pairs of assays; an additional 26 loci were confirmed by two pairs of assays. (b–c) Validation of known and novel DMRs identified from MeDIP-seq and MRE-seq. DMRs are presented in a UCSC Genome Browser window with MeDIP-seq and MRE-seq signals in human H1 ESC, along with bisulfite sequencing results. The results from the biological replicates (#1 and #2) were very similar. (b) Imprinted gene GRB10 including a known DMR (Bisulfite region 1) and an upstream unmethylated CGI (Bisulfite region 2). (c) Novel DMR upstream of POTEB which exhibits allele specific DNA methylation. Open circle indicates an unmethy-lated CpG site. Filled circle represents a methylated CpG site. “x” indicates absence of a CpG site due to a heterozygous SNP which destroyed the 28^th CpG. All clones without the CpG were unmethylated, while all the clones containing the CpG were methylated. Furthermore, the alleles could be distinguished in the sequence reads from MeDIP-seq (G, 9 of 9 reads) and MRE-seq (A, 30 of 30 reads).