Comparison and imputation-aided integration of five commercial platforms for targeted DNA methylome analysis

Targeted bisulfite sequencing (TBS) has become the method of choice for the cost-effective, targeted analysis of the human methylome at base-pair resolution. In this study, we benchmarked five commercially available TBS platforms—three hybridization capture-based (Agilent, Roche and Illumina) and two reduced-representation-based (Diagenode and NuGen)—across 11 samples. Two samples were also compared with whole-genome DNA methylation sequencing with the Illumina and Oxford Nanopore platforms. We assessed workflow complexity, on/off-target performance, coverage, accuracy and reproducibility. Although all platforms produced robust and reproducible data, major differences in the number and identity of the CpG sites covered make it difficult to compare datasets generated on different platforms. To overcome this limitation, we applied imputation and show that it improves interoperability from an average of 10.35% (0.8 million) to 97% (7.6 million) common CpG sites. Our study provides guidance on which TBS platform to use for different methylome features and offers an imputation-based harmonization solution that allows comparative, integrative analysis. Differences in CpG coverage of targeted bisulfite sequencing methods can be overcome using imputation.

D NA methylation is an indispensable epigenetic mark for many biological processes, such as development, differentiation and maintenance of cell-type-specific states 1 . Interrogating the changes in DNA methylation patterns is essential to better understand the biology of normal and pathologic states and to identify clinically relevant biomarkers. In the past decade, methods for methylome analysis have moved away from semi-quantitative methods with coarse resolution (MeDIP and MRE-seq) toward methods with the single-base resolution based on bisulfite conversion, such as microarrays and next-generation sequencing (NGS) 2,3 . NGS offers several advantages, including single-molecule analysis and read phasing to study sample heterogeneity and epiallele composition. 4,5 Whole-genome bisulfite sequencing (WGBS) remains the gold standard method for studying DNA methylation at a single base pair (bp) resolution, although it is associated with high costs and requires considerable computational resources. Targeted bisulfite sequencing (TBS) directs the sequencing to more informative parts of the genome, reducing the costs. The analysis of DNA methylome by TBS can be achieved either through target-specific enrichment of regions containing CpG sites of interest using probe hybridization capture (HC) methods 6 or through non-specific enrichment of CpG-dense regions by reduced-representation bisulfite sequencing (RRBS) mediated by a restriction enzyme recognition site containing a CG motif 7 .
Currently, five commercial manufacturers are producing off-the-shelf kits, offering standardized reagents and conditions for genome-wide TBS using HC (Agilent SureSelect Methyl-Seq, Roche NimbleGen SeqCap EpiGiant and Illumina TruSeq Methyl Capture EPIC) or RRBS (Diagenode Premium RRBS and NuGen Ovation RRBS Methyl-Seq). The five platforms employ different experimental strategies to generate sequencing libraries and differ in the scope of regions they cover. The diversity of these platforms' characteristics, along with the absence of a thorough comparison 8 of their output and performance, warrants a comprehensive benchmarking to provide guidance for users to select the most appropriate platform based on each one's strengths and limitations.
Here we systematically compare performance of these five platforms in terms of sequencing output, target capture efficiency, genomic features coverage, CpG coverage similarity, intra-platform reproducibility, between-platform concordance and differential methylation, across a set of 11 samples. We benchmark each of the five TBS platforms to gold standard data generated using Illumina WGBS and Nanopore sequencing. Finally, we evaluate interoperability of the TBS platforms and provide guidance for their harmonization.

Results
DNA methylation data generation. To ensure reproducibility and data availability under open-access agreement, we selected for the study a set of 11 well-characterized and/or commercially sourced samples. These included common reference standards: Coriell-NA12878 and HeLa cell lines processed in duplicate; a human reference genomic DNA (gDNA) sample isolated from pooled healthy peripheral blood (Ref.gDNA), representing a heterogeneous mix of cell types, at two different DNA inputs (recommended by the manufacturer and 500 ng) in duplicate; four DNA methylation standards generated from commercial (ZYMO) fully methylated and unmethylated control samples; a pair of genetically

Comparison and imputation-aided integration of five commercial platforms for targeted DNA methylome analysis
Sequencing output variability. Sequencing cost-efficiency depends on several factors, including read quality, mapping rates, duplication rates and, in case of hybridization capture methods, target capture efficiency. To simulate real-life biological experimental setup and to enable an unbiased comparison between platforms, equimolar amounts of each sample library were pooled together and sequenced to achieve an average of 30 million uniquely mapped de-duplicated reads (UMDRs). We compared the quality of raw and pre-processed data (Extended Data Fig. 1 and Supplementary Fig. 1) and the technical variability of sequencing output (Fig. 2a and  Supplementary Table 2). To achieve the targeted number of UMDRs, between ~55 million reads (for NuGen) and over ~80 million reads (for Illumina) had to be sequenced. Due to reduced sequence diversity of bisulfite-converted DNA, alignment rates are generally lower than for genome sequencing. Mapping rates for HC-based methods were similar between platforms (~79%) and higher than for RRBS-based platforms (~64%). Bisulfite conversion (BC) is a harsh chemical process that degrades DNA, reducing the library complexity and, if combined with a high number of polymerase chain reaction (PCR) cycles, may lead to a large number of duplicate reads skewing the methylation readout 9,10 . De-duplication can be achieved by relying on fragment genomic coordinates 11 (Agilent, Roche and Illumina) or using unique molecular identifiers (UMIs) 12 (Nugen) but not for RRBS without UMIs (Diagenode). Duplication rates were the highest for Illumina's platform (~52%) and the lowest for NuGenʼs (~7%). Overall, the highest ratio of UMDRs to passing filter (PF) reads was observed for NuGen's platform and the lowest for Illuminaʼs.
Target capture efficiency. Sequencing cost-effectiveness for HC methods also depends on the percentage of sequencing directed to regions of interest versus the amount of off-target sequencing. For this analysis, we omitted RRBS-based platforms because their enrichment is not target specific. The performance of the wet-lab protocol for the three HC platformsʼ on-target capture efficiency was evaluated based on coverage of all targeted bases according to manufacturers' design files (Supplementary Table 3). Overall, on-target capture efficiency measured as the percentage of uniquely aligned PF bases that mapped on or near bait was highest for Illumina's platform (~90.6%), followed by Agilentʼs (~78.2%) and Rocheʼs (~61.5%), with the two latter platforms marked by high inter-sample variations (Fig. 2b). Fold enrichment by which the baited region has been amplified above the genomic background was similar between platforms (21.3-24.3) (Fig. 2c).
Target coverage was assessed as a proxy for how well the data are likely to perform in downstream applications. Mean target coverage was similar for all platforms, from 35.4× (Roche) to 37.4× (Illumina) (Fig. 2d). We observed higher uniformity of coverage for Rocheʼs and Illumina's platforms (~2.3 and ~2.4, respectively) compared to Agilentʼs (~3.3) (Fig. 2e). The uniformity of coverage translates into the fraction of targeted regions being covered at a specific depth (Extended Data Fig. 2). Intra-platform variability was pronounced for all three HC platforms, including technical replicates, with the average percentage of targets covered ≥10× ranging from ~80.9% for Agilent to ~88.1% for Illumina (Fig. 2f).
CpG coverage. An important consideration when choosing a platform for methylome analysis is the number of CpG sites it interrogates. The number of CpGs covered by RRBS is a function of the depth of sequencing and sequencing conditions (read length and single-end or paired-end mode). Although each MspI-generated fragment has a CpG site at its ends, the number of CpGs within each fragment is not uniform due to the random selection of recognition sites by the enzyme. In contrast, for HC methods, the number of CpGs covered is determined by design. However, we observed a much higher number of experimentally covered CpGs at 1× than expected for all three HC platforms, consistent with high levels of off-target sequencing (Fig. 3a). The average number of CpG sites at 5× drops to ~4 million CpGs for all platforms but remains higher for RRBS platforms at 10× and 30× compared to HC methods. Both RRBS platforms showed similar CpG coverage, ranging from over 5 million CpGs at 1× to ~2 million CpGs at 30×. Mean CpG depth of coverage was higher for RRBS methods (Diagenode ~33× and NuGen ~31×) compared to HC methods (Agilent ~13×, Roche ~10× and Illumina ~13×), reflecting a more uniform coverage for the former (Supplementary Table 2).
Next, we determined the degree of CpG coverage similarity between and within platforms, in terms of the number of overlapping CpGs covered ≥10×. For each sample platform, we calculated a pairwise Jaccard similarity index defined as the number of commonly interrogated CpG sites (intersection) over the total number of CpGs (union). Unsupervised hierarchical clustering indicated that, by and large, RRBS-based and HC-capture-based platforms cover different CpG sites (Fig. 3b). All samples clustered according to the platform of origin, with Diagenode showing the highest degree of intra-platform similarity. As expected, NuGen and Diagenode RRBS platforms covered mostly the same CpG sites, whereas, among HC methods, Agilent and Roche platforms had a higher degree of similarity compared to Illumina.
Theoretical and empirical coverage of genomic features. The CpG distribution is not uniform throughout the genome; thus, an important consideration for choosing a platform is the coverage of genomic features of interest. HC platforms target CpG islands (CGIs), shores and shelves, regulatory regions and known differentially methylated regions (DMRs) to a different extent by design. We compared the size of targeted genomic features and the number of  CpGs targeted in each platform ( Supplementary Fig. 2). Illumina's platform targets ~20 megabases (Mb) more than other platforms, including the highest portion of FANTOM5 (ref. 13 )-defined enhancers (>40%) and insulator regions. Agilent's platform targets the largest region and the highest number of CpGs within CGIs, promoters and exons. Interestingly, although Illumina's platform targets the largest portion of inter-CGI, intron, and heterochromatin regions, RRBS encompasses more CpGs in those regions.
To determine how well each platform experimentally covers these regions, we counted the number of CpGs covered ≥10× per platform in each genomic/regulatory feature (Fig. 3c). In general, we observed a pronounced variation for all platforms except Diagenode. Striking differences in theoretical versus empirical coverage were observed, especially for Agilent's platform, reflecting lower uniformity of coverage, with many CpGs not reaching the 10× threshold. Overall, RRBS platforms covered the highest number of CpGs located in CGIs and shores, whereas shelves and open-sea regions were best represented by Roche. Illumina interrogated over twice as many CpGs located in FANTOM5 enhancer regions compared to other platforms. Roche showed the highest coverage of promoters, exons, 5′ untranslated regions (UTRs) and 3′ UTRs, whereas RRBS platforms targeted the most CpGs in introns. When looking at regulatory features defined by ChromHMM 14 where promoters were stratified by their activity, HC methods covered more CpGs in active promoter regions, whereas RRBS methods had higher coverage of weak promoters. There were no significant differences between platforms in the coverage of either strong or weak enhancers, whereas Illumina covered the most CpGs located in insulator regions.

Intra-platform reproducibility and inter-platform concordance.
To estimate the level of confidence in obtained results for any given biological experiment and to set appropriate cutoffs for calling differential DNA methylation, it is important to know the approximate level of technical variation in the data. To assess platform reproducibility, we compared DNA methylation data from samples for which the same DNA was used to generate duplicate sequencing libraries (Ref.gDNA, Coriell and HeLa). DNA methylation levels were calculated as the fraction of methylated CpGs over the sum of methylated and unmethylated CpGs. We constricted the analysis to CpG sites covered ≥10× in both samples to calculate pairwise Pearson correlation coefficients.
All platforms showed a very high overall correlation in DNA methylation levels of replicate samples, with Diagenode slightly underperforming ( Fig. 4a and Supplementary Fig. 3). The effect of DNA input was estimated using Ref.gDNA libraries prepared with 1,000 ng of DNA for HC or 100 ng input for RRBS methods compared to 500 ng common to all platforms. All platforms, apart from Diagenode, showed very high intra-platform and inter-platform correlation in DNA methylation levels (>0.98) for different DNA inputs ( Supplementary Fig. 4).
We next looked at the between-platform concordance in methylation levels for each sample in a pairwise fashion over commonly covered CpGs (94,000-450,000) to enable direct comparison ( Supplementary Fig. 5). Methylation levels for each sample per platform were compared to all other platforms to calculate the average correlation over common CpGs (Fig. 4b). Roche HeLa-1 replicate was excluded from the analysis because the hybridization capture failed. We observed a high correlation in called meth- ylation levels for all samples (R: 0.85-0.99) except for DNAm-5pct, DNAm-10pct and ZYMO fully methylated control. To examine this further, we compared the concordance between platforms for overall DNA methylation levels on a set of DNA methylation standards generated by mixing DNA from fully unmethylated (ZYMO-UM) and fully methylated controls (ZYMO-FM) at known ratios (0% ZYMO-FM, 5% ZYMO-FM, 10% ZYMO-FM and 100% ZYMO-FM). ZYMO-UM control DNA is isolated from DNMT1/DNMT3b knockout cell line known to have less than 5% methylated DNA 15 . All technologies appeared to slightly overestimate methylation levels for ZYMO-UM, and all apart from NuGen tended to underestimate methylation levels of ZYMO-FM (Fig. 4c). This was most pronounced for Diagenode's platform with 83% estimated methylation for ZYMO-FM, most likely due to the lack of a de-duplication step that removes PCR duplicate reads, driving down the inter-platform concordance for these samples.

TBS comparison to WGBS and Oxford Nanopore methylome.
To examine how well each of the more commonly used TBS platforms compares to WGBS as the current gold standard method for near-complete characterization of DNA methylome at single-bp resolution, and Oxford Nanopore as an emerging third-generation sequencing technology that allows for a direct detection of modified bases 16   with very high degree of overlap (~24 million CpGs) and methylation level concordance (ρ ~0.92) (Extended Data Fig. 3), all TBS platforms displayed, on average, a similar correlation to WGBS datasets (0.96 ± 0.02) and to Nanopore data (0.94 ± 0.02), albeit with the higher CpG overlap with the latter (~2.9 million versus ~ 3.4 million).
Calling differentially methylated sites and regions. Ultimately, many studies aim to identify differentially methylated CpG (DMC) sites or DMRs between samples. To assess the between-platform concordance for calling differential methylation, we compared two pairs of bladder cancer cell lines-one pair with a different genotype and diverse phenotypes (253J and T24) and a second pair of isogenic cell lines with divergent sensitivity to cisplatin (RT112 and RT112-CP). Because the observed technical variability was as high as 5% for some platforms, we set the threshold to >10% methylation difference and q ≤ 0.05 to call DMCs and DMRs. The number of DMCs depends on the number of CpGs assayed, their location with respect to genomic features and inherent technical biases for each of the platforms. To calculate DMCs, a set of overlapping CpGs within the pair is analyzed per platform but different sets of CpGs between platforms as a consequence of differences in platform CpG coverage. As expected, the bulk of identified DMCs for all platforms showed larger methylation differences for the genetically distinct pair of cell lines compared to the isogenic pair, with only a small subset of DMCs commonly identified by all five platforms for both pairs of cell lines (Fig. 5a). Surprisingly, for the RT112 versus RT112-CP pair using Agilent's platform data, we failed to call any significant DMCs, although the other four platforms shared most of the identified significant DMCs. Aggregating methylation levels using tiling regions is a strategy commonly used to deal with sparse data where there is no coverage on exactly the same CpGs. When comparing aggregated methylation levels over 1,000-kilobase (kb) tiling windows to call DMRs (Fig. 5b), high concordance was observed among all five platforms for T24 versus 253J cells, with NuGen identifying an additional subset of highly significant DMRs. Again, Agilent's platform did not identify any of the DMRs called by the other four platforms for the RT112 versus RT112-CP pair. Notably, methylation differences for significant DMRs were lower than methylation differences observed for DMCs. Imputation to improve platform interoperability. The platform similarity as a function of the number of overlapping CpG sites is of particular concern as downstream analyses regularly require subsetting to a set of CpG sites that are present in all datasets. Furthermore, the effect of low overlap in the CpG sites selected is further exacerbated as the number of samples/datasets increases. Imputation can be used to estimate the methylation value for each missing CpG site 18,19 and, thereby, can be used to help increase platforms' intra-operability and inter-operability. Samples with duplicate sequencing libraries, one clonal (Coriell-NA12878) and one heterogeneous admixture of blood cell types (Ref.gDNA), were further analyzed to measure the effect of imputation on CpG methylation concordance and overlap. We used BoostMe 19 to impute the vast majority of missing CpG sites, leveraging information learned from two neighboring CpGs of the same dataset, with or without a distance threshold (Fig. 5c-f). Imputation without distance threshold increased the overlap between platforms almost ten-fold, from 0.79 million (10.35%) to 7.6 million (97%) CpG sites (Fig. 5c, Extended Data Fig. 4 and Supplementary Fig. 6). However, it came at the cost of reduced DNA methylation concordance from an average of 0.97 to 0.80 after imputation ( Fig. 5d and Supplementary Fig. 7). Mean absolute error (MAE), which estimates by how much each methylation value is off, increased from an average of 0.05 to 0.14 ( Fig. 5e and Supplementary Fig. 8), and root mean square error (RMSE) increased from 0.09 to 0.20 ( Fig. 5f and Supplementary Fig. 9), whereas MAE corrected for the number of CpGs present in the   dataset decreased from an average of 0.12 to 0.002 ( Supplementary  Fig. 10). When a distance threshold of <25 bp was applied to impute neighboring CpGs, the total number of overlapping CpGs was 2.5 million (32%), whereas the intra-platform concordance after imputation remained high (ρ = 0.89, MAE = 0.10).

Discussion
To answer important biological questions in the epigenetics of health and disease, we must understand the possibilities and limitations of the tools that we use. Over the past decade, many studies harnessed the power of bisulfite sequencing methods to advance knowledge on epigenetics of non-communicable diseases, to identify clinical biomarkers and to understand normal development [20][21][22][23][24]   using scarce clinical samples. NuGen RRBS offers the shortest and least complex protocol that requires only 100 ng of DNA, suitable for high-throughput applications using scarce clinical samples. On the other hand, smaller-scale studies focused on in vitro experiments using cell lines, animal models or abundant clinical specimens have fewer constraints, allowing the use of HC platforms. Among HC platforms, Illumina's protocol can be completed in 2 days and requires half the amount of DNA per sample (500 ng) by multiplexing four samples, albeit reducing the flexibility in experimental design. Differences in design, size of the targeted region and experimental strategies for library preparation and target capture translate into stark differences in the cost-effectiveness of sequencing. To achieve a similar number of uniquely mapped de-duplicated reads, HC platforms required ~20% more sequencing than RRBS platforms. In addition to RRBS platforms, only the Roche HC platform targets both strands at the same ratio, whereas Illumina covers both strands only for common single-nucleotide polymorphisms (SNPs), permitting the distinction between a C-to-T SNP and a methylated CpG site.
Estimates of 5-methylcytosine (5mC) levels may be biased by experimental and bioinformatic factors. Our results reveal high within-platform and between-platform concordance in measured 5mC levels of the overlapping CpG sites. Although all platforms slightly underestimated global 5mC levels for fully methylated DNA standards, this was particularly pronounced for Diagenode's platform, reflected by the lower between-platform concordance of methylated ZYMO controls. Because the sample source and sequencing conditions were the same for all, the bias must have been introduced during the library preparation and/or data de-duplication. Olova et al. 25 found that the BC is the main driver of sequencing biases that are further amplified by PCR amplification. Although all platforms perform adaptor tagging before BC, Agilent and Illumina perform hybridization capture before BC, whereas Roche performs hybridization capture after BC. 26 Although the data processing of HC platforms included a de-duplication step based on read coordinates, only NuGen's platform incorporates UMIs, enabling precise read de-duplication to account for bias introduced by PCR amplification. In contrast, Diagenode's platform does not allow for PCR de-duplication, rendering it susceptible to biased estimation of 5mC levels. Conversion efficiency was above 99% for all platforms except for Roche, where the rate was over 98%. Over-conversion may lead to underestimation of 5mC levels; however, fully methylated spike-in was not included in all platforms, so we could not compare the over-conversion rates.
The choice of method defines the space for the discovery of CpG loci or regions associated with a specific phenotype. Each of the platforms targeted different regions of the genome, and the overlap between CpG loci covered by all was modest. Overall, RRBS methods covered more CpGs ≥10×, which were not targeted by HC platforms. We provide a breakdown of genomic feature coverage by platform that can be used as a guide for platform selection depending on the biological question.
Calling of DMCs is limited on the CpGs covered in common between the samples compared. Given a pronounced variation in the number of covered CpGs within each platform, the larger the number of samples, the lower the overlap. To fully capitalize on a dataset, two strategies can be employed. First, methylation values can be aggregated in defined regions (for example, promoters) or tiling regions (for example, 1,000-bp windows) to call DMRs. We found that, although many of the DMCs called by the five platforms did not overlap, most DMRs were commonly identified by all platforms. Alternatively, to take the most advantage of the base-level resolution offered by TBS, missing CpGs can be imputed. We showed that imputation of missing CpG sites is a viable option to increase inter-platform operability while maintaining high concordance for downstream analysis. Although intra-sample heterogeneity and discordant methylation may pose a challenge for precise imputation, by selecting the stringency of the distance threshold and associated MAE, users can determine the appropriate methylation difference cutoff for calling DMCs.
In summary, our study provides an in-depth analysis of the TBS platforms' comparative performance and characteristics, allowing users to determine the optimal technology for methylome analysis depending on their needs and restrictions. It also provides guidance for cross-platform data integration.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41587-022-01336-9.

Methods
Sample origin and processing. The reference gDNA sample (Ref.gDNA) was obtained by pooling DNA extracted from peripheral whole blood cells of 20 adult healthy volunteers who provided informed consent under 15/YH/0311 of the UCL BioBank for Health and Disease ethics committee. Blood samples were collected in EDTA-treated 10-ml BD Vacutainer tubes and centrifuged at 1,900g for 10 minutes at 4 °C using a refrigerated centrifuge to separate blood cells from plasma. Buffy coat (0.5 ml) was transferred using a Pasteur blob to a clean Falcon tube containing 4.5 ml of HEMAgene•BUFFY COAT DNA stabilizing reagent (Oragene), and the DNA was extracted using QIAamp DNA Blood Mini Kit (Qiagen). The Coriell-NA12878 DNA sample was purchased from Coriell Cell Repositories. HeLa cell gDNA was purchased from New England Biolabs. Bladder cancer cell lines T24, 235J and RT112 were acquired from the American Type Culture Collection (TCP-1020). RT112-CP was a gift from Jim Catto at Sheffield University. gDNA was extracted using QIAamp DNA Mini Kit (Qiagen). The Human Methylated & Non-methylated DNA Set (ZYMO) was used to generate DNA methylation standards: fully methylated (ZYMO-FM), 10% and 5% methylated (DNAm-10pct and DNAm-5pct) and unmethylated (ZYMO-UM). The negative control (ZYMO-UM) DNA in the set comes from a human HCT116 DKO non-methylated DNA that contains genetic knockouts of both DNA methyltransferases DNMT1 (−/−) and DNMT3b (−/−) and has low levels of DNA methylation 15 , whereas the fully methylated DNA control was generated by the manufacturer by treating the HCT116 DKO with M.SssI methyltransferase. Sample quality was determined using an Agilent 2100 Bioanalyzer 1000 DNA chip to assess the DNA integrity, and sample purity was estimated using a NanoDrop spectrophotometer. Quantification was performed using Qubit dsDNA BR Assay (Invitrogen).

Library preparation and target enrichment with Agilent SureSelect Methyl-Seq kit.
To generate SureSelect Methyl-Seq sequencing libraries (Agilent), we followed the manufacturer's protocol. In brief, we sheared 1 µg of the gDNA on Covaris S2 ultrasonicator system to obtain ~150-bp fragments. Fragmented DNA was cleaned up using AMPure XP beads (Beckman Coulter), followed by end-repair, A-tailing, second cleanup and adapter ligation before the final cleanup. Quality and quantity of the libraries were checked using a Bioanalyzer dsDNA High Sensitivity chip (Agilent) and Qubit HS dsDNA Assay (Life Technologies, Thermo Fisher Scientific) before overnight hybridization (minimum 16 hours) with biotinylated oligo RNA baits, followed by streptavidin-conjugated magnetic bead pulldown and wash steps. The captured library was then bisulfite-converted using EZ DNA Methylation-Gold Kit (ZYMO), followed by eight cycles of PCR amplification and AMPure XP bead cleanup. Second PCR amplification with six cycles of the bisulfite-converted libraries was performed to introduce sample indices, followed by final AMPure XP bead cleanup. Library concentration was measured using Qubit dsDNA HS kit, and the library size distribution was checked on a Bioanalyzer High Sensitivity DNA chip.

Library preparation and target enrichment with Roche NimbleGen SeqCap
EpiGiant kit. NimbleGen SeqCap EpiGiant (Roche) libraries were prepared according to manufacturer instructions, as follows: 1 µg of the gDNA sample was sheared using Covaris S2 ultrasonicator system to an average DNA fragment size of 200 bp; sheared DNA was end-repaired and cleaned up with magnetic beads, followed by A-tailing of the 3′ end, second cleanup, adapter ligation and double SPRI cleanup. Libraries were bisulfite-converted using EZ DNA Methylation Lightning Kit (ZYMO), followed by 12 cycles of Pre-Capture LM-PCR amplification using Kapa HiFi Uracil+ Polymerase, cleanup and fragment distribution analysis on a Bioanalyzer High Sensitivity DNA chip. The bisulfite-converted library was hybridized with biotinylated DNA capture probes over 72 hours and washed with streptavidin-conjugated magnetic beads. The enriched libraries were amplified by a second round of LM-PCR with 16 cycles using Kapa HiFi Uracil+ Polymerase (Roche) and purified with AMPure XP magnetic beads. Library concentration was measured using the Qubit dsDNA HS kit; purity was measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific); and the library size distribution was checked on a Bioanalyzer High Sensitivity DNA chip.

Library preparation and target enrichment with Illumina TruSeq Methyl
Capture Epic kit. The TruSeq Methyl Capture Epic kit was a gift from Illumina. Sequencing libraries were prepared according to the manufacturer's instructions. In brief, 500 ng of each gDNA sample was sheared on a Covaris S2 ultrasonicator to a median size of 160 bp, followed by a cleanup step using AMPure XP magnetic beads, end-repair reaction, second magnetic bead cleanup, A-tailing, adaptor ligation and final cleanup. Four samples containing adaptors were pooled together to a total mass of 2 µg and hybridized with a biotinylated DNA capture panel for 35 minutes at 58 °C. Target DNA fragments were captured by streptavidin-conjugated magnetic beads and washed, followed by a second overnight hybridization (touchdown from 95 °C to 58 °C), streptavidin-conjugated magnetic bead pulldown and wash step. Enriched libraries were bisulfite-converted using the EZ DNA Methylation Lightning kit with magnetic bead de-sulphonation and cleanup. Eleven cycles of PCR amplification were performed using Kapa HiFi Uracil+ Polymerase, available separately from the kit. Final libraries were bead-purified to remove adapters and quantified using the Qubit dsDNA HS kit, and the library size distribution was checked on a Bioanalyzer High Sensitivity DNA chip.

RRBS library preparation using Diagenode's RRBS Premium kit.
The protocol for Premium RRBS kit (Diagenode) library preparation included enzymatic digestion of 100 ng of gDNA with MspI, followed by end-repair, adaptor ligation and size selection by bead purification. Each library was quantified by quantitative PCR (qPCR) to determine library concentration, and eight samples were pooled in equimolar amounts determined by an Excel pooling tool provided by the manufacturer. Pooled libraries were bisulfite-converted using manufacturer-provided BS Conversion Reagent and de-sulphonated on columns. The second qPCR was performed to determine the optimal number of cycles for the final amplification step for each pool (11 and 12 cycles) before final cleanup. Library concentration was measured using the Qubit dsDNA HS kit, and the library size distribution was checked on a Bioanalyzer High Sensitivity DNA chip.

RRBS library preparation using NuGen's Ovation RRBS Methyl-Seq System 1-16.
The Ovation RRBS Methyl-Seq System 1-16 (NuGen, now Tecan) user manual was followed to generate sequencing libraries by enzymatically digesting 100 ng of gDNA using MspI, followed by end-repair, adaptor ligation and a final repair step. Generated libraries were bisulfite-converted using EpiTect Fast DNA Bisulfite Kit (Qiagen) purchased separately from the kit. Converted libraries were amplified by PCR using 12 cycles and purified using Agencourt RNAClean XP magnetic beads, and fragment distribution was checked on a Bioanalyzer High Sensitivity DNA chip. An RRBS design BED file was created by in silico cutting of hg19 sequence at CG within CCGG motif and filtering to >100-bp length ranges and 100-bp window from each end.
Sequencing on an Illumina HiSeq 2500 instrument. All generated, targeted bisulfite-converted libraries were sequenced on the HiSeq 2500 using HiSeq v4 SBS reagents (Illumina) at 2 × 100 bp. Barcoded samples per each library were multiplexed, and libraries were denatured with NaOH and diluted to a final concentration of 2 nM. Denatured libraries were loaded onto cBot for cluster generation on separate flowcell lanes in HighOutput or RapidRun mode according to the manufacturer's protocol. Loading concentrations and percent of PhiX spike-in libraries for each platform are as follows: Agilent at 15 pM with 5% of equimolar PhiX, Roche at 16.5 pM with 5% equimolar PhiX, Illumina at 11 pM with 10% equimolar PhiX, Diagenode at 16 pM with 5% of 12-pM PhiX and NuGen at 11.5 pM with 5% of equimolar PhiX. Sample pools were sequenced over several lanes/runs, and raw FASTQ reads were merged before the analysis.

TBS data analysis.
Raw sequencing data processing. For all targeted bisulfite sequencing platforms, Bismark 27 data analysis pipeline was applied with few platform-specific modifications. Samples were demultiplexed either onboard Illumina HiSeq2500 or offline using bcl2fastq Conversion Software version 1.8.4 (Illumina). The quality of raw reads in FASTQ format was checked using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) version 0.11.9 ( Supplementary Fig. 1).
For each sample platform, FASTQ files from different sequencing runs were merged before alignment. NuGen FASTQ reads contain an additional 1-3 inserted nucleotides and were trimmed using a diversity-trimming Python script provided by the manufacturer (https://github.com/NuGentechnologies/NuMetRRBS). Any contaminating adaptor sequence was removed using TrimGalore (https://github. com/FelixKrueger/TrimGalore TrimGalore) version 0.6.5 in paired-end mode with default conditions for adaptors and -trim1 option for Agilent, Roche and Illumina, whereas, for Diagenode's platform, -rrbs option was used, and, for NuGen, read 2 adaptor was specified. Bismark version 0.10.1 (ref. 27 ) was used for alignment to the human genome build hg19 in paired-end mode using Bowtie2 version 2.3.4.2 and default settings.
After alignment, PCR duplicates in Agilent, Roche and Illumina were marked and discarded using deduplicate_bismark perl script. For RRBS data, this form of de-duplication is not recommended as all reads start with the same sequence. However, the NuGen platform contains UMIs that enable for precise removal of duplicates using Duplicate Marking Python script provided by the manufacturer (https://github.com/NuGentechnologies/NuMetRRBS). Bismark's module bismark_methylation extractor with options --ignore_r2 2, --no_overlap, --comprehensive, --merge_nonCpG, --bedGraph and --gzip was used to generate output files with methylation status for each individual CpG, coverage and genomic coordinates. Bismark's coverage2citosine tool with and without --merge_ CpG option was used to obtain methylation status for each CpG in the genome, both strands, along with trinucleotide context, strand, coverage and methylation percent. The efficiency of bisulfite conversion was verified by the ratio of C-to-T conversion of CHG and CHH (non-CG) dinucleotides. MultiQC 28 version 1.10.1 was used for visualizing aggregated quality metrics.
DNA methylation data analysis. Downstream analysis was performed in Unix (RHE) and R environment (version 3.4.1). We did not use Repeat masking for any of the analyses. For hybridization capture methods, the Picard (https://broadinstitute.github.io/picard/index.html) CollectHsMetrics tool (version 2.3.0) was used to obtain statistics on HC efficiency. The uniformity of coverage of the target region was represented as the theoretical increase of non-zero read coverage necessary to raise 80% of bases to the mean coverage levels in those targets (fold-80 penalty), with lower values indicating less variability in target coverage. BEDTools 29 and custom sed and awk scripts were used to calculate and visualize the target region coverage as a function of read depth. To compare target breadth of coverage, UpSet plot from ComplexHeatmap 30 and Venn diagram from ChipPeakAnno 31 R packages were used to visualize the intersection of genomic ranges and CpG sites.
Bismark-produced strand-specific methylation calls were imported into the MethylKit 32 R package for downstream analysis, including CpGs with more than ten reads. To get per-CpG methylation levels independent on a strand, counts on both strands of a CpG dinucleotide were merged. To perform pairwise comparisons in methylation levels, samples of interest were merged into one object containing common CpGs. Pearson correlation was used to assess concordance in DNA methylation levels in a pairwise comparison matrix. Differential methylation analysis was performed using Fisher's exact test or logistic regression to calculate P values per CpG and over 1,000-kb tiling windows, corrected by the SLIM method 33 to obtain q values. DMRs/CpGs were selected based on q < 0.05 and percent methylation difference larger than 10%.
Annotation. Feature annotation was done using anotatr 34 and analyzed with the GenomicRanges 35 R package. The CpG islands for hg19 were retrieved through the AnnotationHub 36 package, and CpG shores were defined as 2 kb upstream/ downstream from the ends of the CpG islands, less the CpG islands; CpG shelves as another 2 kb upstream/downstream of the farthest upstream/downstream limits of the CpG shores, less the CpG islands; and CpG shores. The remaining genomic regions make up the inter-CGI annotation. The genic annotations are determined by functions from GenomicFeatures 35 and data from the TxDb and org.hg19.eg.db packages. Genic annotations include 1-5 kb upstream of the transcription start site (TSS), the promoter (<1 kb upstream of the TSS), 5′ UTR, first exons, exons, introns, coding sequences, 3′ UTR and intergenic regions (the intergenic regions exclude the previous list of annotations). FANTOM5 permissive enhancers determined from bi-directional CAGE transcription 13 were added to genic annotations. Chromatin states were determined by hidden Markov models (chromHMM) on ENCODE histone mark data for nine cell lines: Gm12878, H1hesc, Hepg2, Hmec, Hsmm, Huvec, K562, Nhek and Nhlf. These were used to check the theoretical coverage of regulatory features by each platform 14 . GRangesList object containing all cell lines data was stratified by chromatin state. Because cells have different regulatory states, overlapping ranges within each chromatin state were reduced ignoring strand so that each state in the plot represents the union of ranges for all cell types. One-way analysis of variance (ANOVA) in R was used to test for statistically significant differences between group means.
To calculate the breadth of coverage for the design regions determined based on manufacturer-provided BED files with annotation tracks, overlapping region pairs were intersected ignoring strand. To count the number of CpGs in regions of interest, the corresponding hg19 sequence was retrieved from BSgenome. Hsapiens.UCSC.hg19 using the DNA strings R package, reduced to contain only non-overlapping regions independent of DNA strand, and number of occurrences of CG was counted. The number of CpGs covered at a specific depth per platform was calculated using a custom R script. Jaccard similarity index was calculated using functions from the HelloRanges 37 R package. The Intervene shiny app was used to visualize the Jaccard similarity index matrix 38 .
WGBS data. WGBS library for Ref.gDNA sample (500 ng) was generated using the TrueMethyl Whole Genome Workflow Kit (CEGX) with an average insert size of 600 bp. Library was sequenced on Illumina X Ten at 2 × 150 bp with 1,000 million PE reads per lane. Read 2 was of low quality, reducing the mapping down to 35% in PE mode. Therefore, only read 1 was used for mapping with a 78% alignment rate. FASTQ files for the Coriell-NA12878 samples sequenced on Illumina HiSeq 2500 in 2 × 125-nucleotide mode were downloaded from the ENCODE Project 39 database (WGBS_EC) and from Illumina BaseSpace Hub Methyl-Seq. Eight replicates of Coriell-NA12878 were prepared with the TruSeq Methylation kit, run on one flowcell of the HiSeq 4000 at 2 × 76 (WGBS_IL) and reprocessed using the nf-core methyl-seq nextflow pipeline (https://nf-co.re/methylseq). The dataset was filtered to exclude CpG sites with low (<10×) or extremely high (>3 s.d.) coverage in comparison to other platforms.
Nanopore sequencing and data processing. One microgram of gDNA from the lymphoblastoid Coriell-NA12878 cell line was sheared using g-Tube (Covaris) to obtain high-molecular-weight fragments ranging from 50 kb to 160 kb. Fragment ends were repaired using NEBNext End Prep kit (New England Biolabs), and AMX adapters were introduced using LSK109 ligation sequencing kit (Oxford Nanopore). Approximately 150 GB of data were generated using two flowcells on PromethION. Nanopore data were processed using the nanopolish 18 pipeline. The dataset was filtered to exclude CpG sites with low (<10×) or extremely high (>3 s.d.) coverage in comparison to other platforms. Fig. 1 | Sequencing data processing quality metrics produced by MultiQC. a Bismark alignment rates for uniquely, ambiguously, or unaligned reads for each sample by platform; b Percent of reads aligning to top or bottom DNA strand for each sample by the platform; c Global methylation levels of CpG dinucleotides for each sample by the platform; d The global cytosine methylation level in CHG context for each sample by the platform used an estimate of sodium bisulfite under-conversion rates; e The global cytosine methylation level in CHH context for each sample by the platform used an estimate of sodium bisulfite under-conversion rates.; f -g M-bias plot shows the average percentage methylation and coverage across read length for each sample. Each line represents a sample. Methylation bias for the forward sequencing read by platform (f); Methylation bias for the reverse sequencing read by platform (g).