Rarefaction of EM-seq and WGBS libraries to avoid analytical bias

Analysis of sequencing coverage revealed that EM-seq libraries had higher CpG coverage than sample-matched WGBS libraries. The modes of WGBS libraries ranged from 8–12×, while EM-seq had much higher modes of 10–40× (Fig. 1A). To reduce degree of coverage as a confounding factor for inter-library comparisons, raw reads—pre-trimmed, pre-mapped—were rarefied to match the coverage of the shallowest library (166 million reads; see Materials and Methods). This resulted in coverage patterns that were roughly equal among all libraries (8–10×; Fig. 1A; Supplementary File S1). All downstream analyses were based on these rarefied libraries.

Open in a separate window

Fig. 1

Methylation readouts from EM-seq libraries have better cross-sample correlations than WGBS libraries. (A) Original and rarefied coverages of EM-seq libraries (solid line) compared to WGBS libraries (dashed line), coloured by sample. (B) Pearson correlation was performed on rarefied data, comparing individual EM-seq (ending “-ER”) and WGBS (ending “-WR”) samples for 4 million randomly sampled CpG methylation beta values. The bottom left triangle contains pairwise sample comparisons represented as a 2D density plot. Red lines are lines of best fit, while the background shading indicates relative CpG density (blue: low; yellow: high). These sample comparisons are represented as Pearson correlations in the top right triangle, with *** indicating p-values<0.001. The diagonal represents the distribution of CpG methylation beta values for individual samples. (C) Pearson correlation comparisons were separated by sample (x-axis) and grouped into EM-seq (left panel) and WGBS (right panel) library types. Comparison of samples from the same patient are coloured blue, while comparisons against different patients are coloured red

EM-seq libraries are more consistent and better covered in high GC regions than WGBS

All eight libraries had a similar bimodal distribution of methylation beta values, with heavy concentration of beta around 0 and 1 (0 indicates a fully unmethylated position, and 1 a fully methylated position; Fig. 1B), in line with previous observations [34]. To assess inter-library beta value correlations, we calculated Pearson correlation values from 4 million randomly sampled CpG sites. This analysis revealed significant positive correlations of methylation beta values between all samples, irrespective of the library preparation method (r=0.826–0.906, all p values<0.001; Fig. 1B). These results closely match those obtained using all 28.7 million CpG sites for Pearson correlation analysis, suggesting that randomly sampling 4 million CpGs sufficiently represents the data (Supplementary Fig. S1). In addition, we compared the delta beta values of all 28.7 million CpGs between EM-seq and WGBS, revealing that 95.26% of CpG sites exhibited similar methylation values (delta beta<0.15, Supplementary Fig. S2). We further separated the Pearson correlations values by individuals and by library type. EM-seq libraries had significantly higher intra-method correlations than WGBS (mean r±s.e.m., 0.885±0.007 versus 0.844±0.007, two-tailed t-test p<0.01; Fig. 1C). Also, as expected, correlations of samples from the same patient (blue dots) were significantly higher than correlations with samples from different patients (red dots; 0.887±0.011 versus 0.854±0.007, two-tailed t-test p<0.05; Fig. 1C).

Read GC content in EM-seq and WGBS libraries were analysed to identify whether libraries had differences in coverage in areas of high and low GC content. In line with previous studies [19, 20], EM-seq had elevated normalised coverage compared to WGBS in high GC regions (55–95 GC%) while WGBS had slightly higher normalised coverage in low GC regions (20–35 GC%; Supplementary Fig. S3). In general, both libraries had low coverage in extreme low and high GC content regions (0 and 100 GC%; Supplementary Fig. S3).

CpGs with significantly different betas have strand-specific and motif biases

We sought to identify CpG dinucleotides with significant beta value differences between EM-seq and WGBS libraries, and better understand factors associated with large discrepancies in estimated DNA methylation rates. From 28,704,358 CpGs analysed, only 124 CpGs (0.00043%) had significant differences in methylation betas across the library preparation methods (Supplementary File S2), indicating that both WGBS and EM-seq libraries overall have high CpG methylation beta concordance. We further analysed the discordance by clustering the beta values from the 124 CpGs in a heatmap with four other computed metrics representing strand-specific methylation, overall coverage, strand-specific coverage, and the immediate sequence context (3 bp either side of a CpG site). The clustering produced two distinct groups: one with significantly higher betas in WGBS than in EM-seq, and another vice versa (Fig. 2A).

Open in a separate window

Fig. 2

Strand-specific and motif biases associated with discordant methylation readouts. (A) Beta values from CpGs that were significantly different between EM-seq and WGBS libraries clustered into two groups in a heatmap. Top annotation bars represent library type, EM-seq (green) and WGBS (purple). Each column represents a different library; each row represents an individual CpG site. Dinucleotides with higher betas in WGBS are above, while those higher in EM-seq are below. Overall and strand-specific methylation levels are represented by “Beta” and “Absolute delta beta” respectively; overall and strand-specific coverages are represented by “Log2 coverage” and “Evenness” respectively. Sequence logos of 8 bp CpG motifs for both groups are placed beside the left dendrogram, while GC% of the motifs are on the rightmost column (“Motif GC%”). White boxes in the heatmap represent missing values. (B) Discordant CpGs were initially split by GC% of 8 bp motifs and subsequently by library type, EM-seq (green) and WGBS (purple). Individual bins were plotted against coverage (top panel) and beta (bottom panel). Sequence logos for each GC% bin is shown below the x-axis

We sought to understand whether these discordant CpGs had strand-specific methylation readouts that departed from the assumed norm of being symmetrical (i.e., betas of Watson and Crick strands should be identical). We calculated two metrics for each CpG, one for beta and one for coverage, which we termed “absolute delta beta” and “evenness”, respectively. The former was calculated using the absolute differences in the strand-specific beta values from the Watson and Crick strands. We observed that absolute delta betas correlated with the methylation state at the CpG site, irrespective of library type (Fig. 2A). This previous finding is likely linked to the latter metric “evenness”, which was calculated using the ratio of the maximal single-strand coverage divided by combined coverage. If coverage was perfectly evenly split across Watson and Crick, then “evenness” is 0.5; if coverage was completely dominated by one strand (perfectly uneven), “evenness” is 1. The majority of differentially methylated CpG sites had coverage confined to a single strand, again irrespective of library type (Fig. 2A). Finally, coverages appeared to be lower for CpGs with significantly higher betas in WGBS (Fig. 2A). This suggests that the differential methylation readouts across both methods is an artefact resulting from coverage differences, where reads map asymmetrically to a single strand.

Next, we investigated whether the discordance in methylation readouts is linked to immediate sequence context. For each of the 124 CpGs, we derived an 8 bp motif comprising the CpG site with flanking 3 bp sequences (5’-NNNCGNNN-3’). Typically, CpGs with significantly higher betas in WGBS have motifs with high GC (Fig. 2A). In addition, these CpGs have a slight propensity to be flanked with “C” and “G” in the third and sixth positions of the motif. Conversely, CpGs with significantly higher betas in EM-seq have motifs with lower GC (Fig. 2A).

To further explore whether local GC content affects methylation betas and coverages, the 124 significant CpGs were split by the GC% value of their corresponding 8 bp CpG motifs (into seven bins) and by library preparation method (EM-seq or WGBS). EM-seq libraries had overall higher median coverage and lower methylation betas than WGBS, apart from coverage for motifs with GC% of 25% (Fig. 2B). Furthermore, these coverage differences were largest in motifs with high GC (62.5–100%), and very low GC (25%). This relationship is reflected in the methylation betas, where lower coverage (for GC% of 62.5–100%) resulted in more binary methylation betas of 0 and 1 for WGBS, compared to EM-seq which had a more even spread of methylation betas (Fig. 2B). If binary methylation betas were excluded, the remaining non-binary beta values were evenly distributed for EM-seq, compared to WGBS where values were mainly confined to betas of <0.5. Interestingly, motifs with low GC (25–50%) preferred to be flanked by either “A” or “T” 3 bp homopolymers. Finally, motifs with 50% GC percentage have the most even beta and coverage between library preparation methods, and strong preference for adenine in the −1 position (5’-NNACGNNN-3’). Overall, there appears to be a relationship between coverage and motif GC% which contributes to differences in methylation readouts between EM-seq and WGBS.

We also conducted a differential analysis prefiltered to only consider CpG sites with delta beta>0.1 (EM-seq - WGBS), reducing the total from 28,031,345 to 3,576,895 sites. The goal here was to only examine sites with some difference and to reduce the penalty for multiplicity correction. This filtering increased the number of significant sites from 124 to 292. Despite this increase, replotting the filtered data (see Supplementary Fig. S4) confirmed that the overall patterns and conclusions remain consistent: coverage asymmetries and motif GC% primarily drive the observed beta differences between EM-seq and WGBS.

EM-seq readouts do not support previously reported TET2 biases

Echoing our motif-related observations, two recent studies have experimentally characterised the preference bias of TET enzymes for specific motifs surrounding methylated CpG dinucleotides [35, 36]. Specifically, TET2 (a component of EM-seq) has a binding preference for the 4-mer motif 5’-MCGW-3’, where M=A/C and W=A/T. Biases in TET2 activity are a concern, as its failure in catalysing 5-methylcytosine to 5-hydroxymethylcytosine will lead to APOBEC converting 5-methylcytosine to uracil (i.e., base will be misclassified as unmethylated).

We analysed the dataset hoping to observe MCGW-driven differences in methylation properties that is exclusive to EM-seq and absent in WGBS. We were not able to find any meaningful differences in beta values, absolute delta beta, and evenness (Supplementary Fig. S5). Overall, our dataset indicates that the beta readouts from EM-seq do not support the previously reported TET2 preferences for MCGW.

Methylation readouts from EPIC arrays equally comparable to EM-seq and WGBS

Using a principal components analysis, we observed that the variance in betas were predominantly due to library preparation method, not by sample ID (Supplementary Fig. S6). To simplify downstream analyses, we computed mean betas across all four samples from each method (EM-seq, WGBS and EPIC; circles in Supplementary Fig. S6). Readouts were largely similar: strong correlations were observed from all three between-method pairwise comparisons (r>0.96, Fig. 3A, C and E). Whilst differences are slight, the short-read methods do agree with each other best (r=0.97), followed by EPIC and EM-seq, and lastly EPIC and WGBS (both r=0.96). The comparisons involving EPIC data had the greatest departures from the identity line (y=x), likely due to EPIC not being able to report perfectly unmethylated or methylated betas as 0 and 1 respectively due to autofluorescence during the scanning procedure (observed betas ranged from 0.0059 to 0.9933), and the addition of a default “stabilising constant” of 100 to the denominator in beta calculations (which prevents division by zero, but beta will always be <1).

Open in a separate window

Fig. 3

Per-position methylation levels assayed using EM-seq, WGBS and EPIC are well correlated and mostly independent of GC% context. (A) and (B) are for EM-seq vs. WGBS; (B) and (D) are for EPIC vs. EM-seq; while (E) and (F) are for EPIC vs. WGBS. Each point in the plots on the left (A, C and E) represent a single cytosine, and its position on the plot corresponds to the mean methylation level assayed using the method labelled on the axes. It is further coloured by the GC% of its local sequence context (±50 bp) to demonstrate context-dependent biases in methylation level measurements. For each point, the residual (differences in assayed methylation levels) was computed and plotted against the same GC% value on the right (B, D and F). All three trend lines are flat and close to origin. This implies that for well-covered cytosines, methylation readouts were almost independent of GC% context that are typical to these three methods (20–80%). At high GC values (>75%), methylation levels appear to be higher in EPIC than either of EM-seq or WGBS (more points above the trendline in D and F; red points in the top right-hand corner of C and E)

Overall, GC% context had little-to-no influence on the eventual readouts of most positions common to these three methods. We defined local GC context as ±50 bp around the methylated cytosine. Between-method methylation readouts were mostly independent of typical GC contexts (20–80%; Fig. 3B, D and F). At high GC (>75%), EPIC overestimates methylation levels relative to both EM-seq and WGBS (Fig. 3C, D, E and F); or both short-read methods were underestimating actual methylation levels relative to EPIC. However, as there were few positions (n=235) with GC% contexts of >75%, broader conclusions could only be drawn from loci-specific experiments on a GC-rich loci.

EPIC probes were designed to target a subset of methylated CpGs in the human genome with cost-effectiveness in mind. As such, targeted positions are primarily located in biologically relevant and non-repetitive regions: mostly within gene bodies, or in promoter regions. Methylation readouts from all three methods were thus contrasted in a context-dependent manner (Fig. 4). Whilst correlation between all three methods were high in all six studied genomic contexts (r>0.93), the short read methods (EM-seq and WGBS) were again in closer agreement than either was with EPIC (Fig. 4), like the earlier context-insensitive per-position analysis.

Open in a separate window

Fig. 4

Methylation levels are well correlated in all biologically relevant genomic contexts across EM-seq WGBS and EPIC. Positions were annotated and binned by their genomic contexts. The first three categories are exclusive—all positions are either in a CpG island, CpG shore, or neither. The next three categories are not exclusive as the human genome contains overlapping genes: a position could be in the gene body of the upstream gene, and in the promoter region of the downstream gene. Across all biologically relevant contexts with varying beta value and GC% ranges, the short read methods (EM-seq and WGBS) produced beta values that were in better agreement with each other than with EPIC.

EPIC vs. EM-seq vs. WGBS

Microarrays are a cost-effective tool for obtaining methylation information from larger sample sizes. Methylation readouts from EPIC have been shown to strongly correlate with both WGBS [15, 37] and EM-seq [37], with correlation coefficients all exceeding 0.85 from the cited studies. Our work builds upon these studies by confirming the strong correlations, and further framing methylation readouts across these three methods with local GC context.

As the EPIC protocol uses bisulphite-converted input DNA, like WGBS, we anticipated that beta values from EPIC would be more similar to WGBS instead of EM-seq. We observed that EPIC–EM-seq consistently had fractionally higher correlation values in all analysis than EPIC–WGBS. These observations were recapitulated in a similar analysis within the Supplementary section of Foox et al. [37], where EPIC–EM-seq had r=0.853 and EPIC–WGBS had r=0.852 (Swift WGBS was labelled as “MethylSeq”). While differences are minor, coverage is likely the reason driving this difference. As EM-seq has a more uniform coverage distribution across the genome than WGBS [19, 20, 37], beta values from the former are less “digital” (weaker clustering around common fractional values e.g., 0, 0.25, 0.33, 0.5, 1). This produces higher correlation values against “analogue” EPIC readouts (fluorescence readouts rarely cluster around common fractions).

For the subset of positions (n=103,670) with sufficient data across all 12 datasets from three methods, betas were independent of window GC% (Fig. 3). However, for a small number of probes (n=235) at high GC (>75%), EPIC overestimates betas relative to both EM-seq and WGBS. This could be due to probes being cross-reactive, i.e., they hybridise to, and produce readouts from, off-targets elsewhere in the genome [15, 41]. Amongst the high GC probes, there were significantly more cross-reactive probes (n=20 out of 235, 8.5%) than expected (n=4,805 of 103,670, 4.6%; p<0.01, Fisher’s exact test). To determine the broader biological implications for these observations, we extended the correlation analysis by grouping positions by their annotated genomic contexts. Due to the technical limitations leading to EPIC being less capable to report extreme beta values (of 0 and 1), the GC-rich, lowly methylated CpG islands had higher beta readouts in EPIC than both short-read methods. Overall, readouts remain strongly correlated across all genomic contexts in all three methods (r>0.93), with EM-seq–WGBS having the strongest correlations.

ONT vs. EM-seq vs. WGBS

The direct detection of methylated cytosines from ONT reads is still undergoing rapid improvements. This is not limited to only improvements in hardware including pore chemistry and library consumables, but also software such as modified base callers [28, 29]. Furthermore, recent variations in ONT sequencing e.g., Cas9-driven enrichment [42] and adaptive sequencing [43] provides alternative avenues to enrich for loci of interest, enabling higher coverages (and thus more accurate methylation readouts) for little-to-no additional effort.

By utilising an amplification-free, Cas9-driven enrichment of a GC-rich rDNA loci, we were able to compare methylation readouts from ONT direct detection against EM-seq and WGBS in this challenging context. In addition, the selection of a multicopy gene resulted in coverages that exceeded 500× for all three techniques, which increases accuracy of the methylation values. Our analysis revealed that EM-seq and WGBS had coverage values that were inversely correlated to local GC%, with WGBS having much stronger dips in coverage within extremely GC-rich regions (Fig. 5; Supplementary Fig. S9D, S9F). Of importance, we observed that ONT coverages remain unaffected by the same GC context bias, due to the long-read nature of this technology (Fig. 5). Previous literature also supports our findings, in the context of ONT sequencing of bacterial genomes with varying GC content [32]. Therefore, ONT is the optimal method for conducting studies involving genomes with large variations in GC contexts across the genome. In addition, shallower sequencing (<30× coverage) due to cost considerations would still provide adequately informative methylation readouts in GC-rich contexts.

One challenge not commonly articulated about ONT modified base calling is that calls are based on confidence, rather than binary calls expected from sequencing converted DNA. For now, most methylation callers set confidence thresholds (e.g., 0–20% as unmethylated, 20–80% as undetermined, 80–100% as methylated) to generate beta values for comparative purposes. For this work, as the comparisons involved methods that generate binary calls, we used the default confidence thresholds of megalodon to generate ONT beta values. Correlations for ONT–EM-seq (r=0.58) and ONT–WGBS (r=0.54) were lower than EM-seq–WGBS (r=0.77), most likely driven by greater technical commonalities between EM-seq and WGBS. Another reason would be that ONT readouts are the least biased by GC-context, followed by EM-seq and WGBS, with the latter having the most severe bias. This can be seen by uncertainties in measuring betas, quantified as the discrepancies in beta values (Fig. 6B, D and F), rising quickly in extremely GC-rich (>75%) contexts. The overestimation of betas by WGBS in GC-rich contexts is in line with previous work that observed greater recovery of fully methylated fragments than fully unmethylated ones after bisulphite treatment [21].

It is important to note that the calling of modified bases from raw ONT signals is a rapidly developing field. Beta values differ based on methylation caller [29] and modified base models. Our current analysis used “remora”, a model trained purely on synthetic datasets (M.SssI-converted methylated DNA and PCR-amplified unmethylated DNA). Previously, when we used “rerio”, trained on a mix of WGBS data and synthetic datasets, ONT readouts behaved more like WGBS at extremely GC-rich contexts (ONT, personal communication; data available on GitHub). While future algorithmic refinement and improvement in ONT pore/chemistry could affect direct detection of modified bases, we remain confident that ONT methylation readouts would be less affected by GC-context biases than short-read methods.

Practical considerations across all four methods

Ultimately, the choice of detection method is dependent on the biological question or translational needs. We focused on the relative performances of each method in GC-rich regions due to their relevance in cancer biology [22–24] and monogenic disorders [38]. We recommend selecting the method which is most cost-efficient and produces the highest quality data, especially when the question or need does not involve GC-rich regions. For brevity, this comparison is presented as a table (Table 1).

Table 1

Picking the right tool for the job. Practical considerations involved in all four methods, as well as their relative strengths and weaknesses

Criteria	EM-seq	WGBS	EPIC	ONT
Flexibility in DNA conversion	NEB-only	Any bisulphite conversion kit		N/A
Flexibility in library construction	NEB-only	More options	Illumina-only	ONT-only
Flexibility in sequencing	Illumina-only	Depends on library type	Illumina-only	ONT-only
Experimental complexity	Well-established protocols which can be performed by trained scientists.			Protocols actively being developed and slightly more complex.
Data analysis complexity	Robust and mature packages/pipelines available.			Pipelines are still in flux.
Turnaround time (from DNA extracts)	2–4 days		3 days	1–2 days (data is streamed)
Relative costs (per sample)				(for ONT Cas9) (for whole genome)
Strengths	Cheaper than WGBS. Coverage more evenly distributed across genome. Data quality better from GC-rich loci than WGBS.	Easier to compare against publicly available data (most are WGBS/RRBS). Bisulphite conversion (without library building) cheaper than enzymatic conversion, better suited for translation into amplicon-based assays.	Very cost effective for getting a subset of methylated and biologically relevant positions across more samples. Ideal for model organisms.	Almost unbiased coverage regardless of context. Quickest turnaround time. Least affected by GC-context biases.
Weaknesses	Increased laboratory time than WGBS. Comparisons against existing data should consider readout divergences at GC-rich loci.	Coverage and methylation readouts biases very pronounced at GC-rich loci.	Not very practical for non-model organisms. Custom panels possible but less cost effective and less reliable.	Higher inputs required. Methylation data from whole genome possible, but more costly. Methylation calls are not binary, unlike bulk of existing data. Higher complexity in sequencing and in analysis.

Open in a separate window

To conclude, both EM-seq and WGBS produces concordant methylation readouts, making either method a reliable choice for obtaining accurate data across the majority of cytosines in most genomes. While construction of EM-seq libraries is more time-consuming, they offer enhanced accuracies for readouts in GC-rich loci or genomes. On the other hand, EPIC presents a more cost-effective option, with the limitation of interrogating a smaller, predetermined set of positions in specific organisms. For researchers seeking unbiased coverage and rapid results, ONT emerges as a promising technology, particularly suitable for extremely GC-rich loci. However, it should be noted that ONT requires a slightly higher level of laboratory expertise and its analytical pipelines are still in the process of maturation.

Methods

Sample collection and DNA extraction

Four whole blood samples were obtained from two participants, WR025 and WR069 (both females, aged 32 and 28 respectively), at two collection timepoints: V1, at the start of the study and V9, 16 weeks into the study. Samples were labelled with participant ID and visit ID, i.e., “WR025V1”, “WR025V9”, “WR069V1” and “WR069V9”.

DNA was extracted from approximately 3 ml of whole blood using the Gentra Puregene Blood Kits (#158,467; Qiagen, Hilden, Germany), following the “Whole Blood” subsection in the manufacturer-provided handbook. The optional RNase digestion step was carried out with the kit-provided RNase A solution. DNA quality and yields were assessed on a NanoDrop 1000 and Qubit 4 Fluorometer with Qubit dsDNA BR Assay kit (Thermo Fisher Scientific, Waltham, MA).

Genome-wide/loci-specific methylation profiling with EM-seq and WGBS

For each sample, 1 µg of genomic DNA was pooled with 1:20 dilutions of unmethylated lambda (1 µl of 0.1 ng/µl) and methylated pUC19 control DNA (1 ul of 0.005 ng/µl) from the EM-seq kit (E7120L; NEB, Ipswich, MA). Volumes were made up to 50 ul with 0.1x TE buffer (Sigma-Aldrich, Burlington, MA). The pooled DNA samples were sheared to an average insert size of approximately 350 bp to 400 bp, using a Bioruptor (UCD400; Diagenode, Denville, NJ) on the high setting for 20–35 cycles (30 s on, 30 s off). Sheared DNA were run on a gel, and samples were subjected to additional shearing cycles if they were not sufficiently sheared. DNA concentrations were again quantified with a Qubit 4 Fluorometer (Thermo Fisher Scientific, Waltham, MA) to ensure DNA inputs for library creation were consistent across samples.

For EM-seq, 200 ng of sheared DNA was processed using the NEBNext Enzymatic Methyl-seq Kit (E7120; NEB, Ipswich, MA) following the manufacturer’s instructions for large insert libraries. For WGBS, 100 ng of sheared DNA was bisulphite-converted using the EZ DNA Methylation-Gold Kit (D5005; Zymo Research, Irvine, CA). Subsequently, the converted total DNA was processed using the Accel-NGS Methyl-Seq DNA Library Kit (#30,024; Swift Biosciences, Ann Arbor, MI). EM-seq and WGBS libraries were quantified using the KAPA library quantification kit (KK4854; Roche Molecular Systems, Pleasanton, CA) and pooled in equimolar amounts. Pooled libraries were sequenced using the NovaSeq 6000 S4 2×150 bp flowcell (Illumina, San Diego, CA) at the Ramaciotti Centre for Genomics (UNSW, Sydney, Australia), aiming for 30× coverage in perfectly pooled samples.

Analyses comparing WGBS against EM-seq in e.g., mapping rates, coverages, dinucleotide compositions, and per-position methylation levels are confounded by sequencing depth—for example, beta values are more accurate with greater sequencing depths. To remove this confounding effect, we equalised the sequencing depths of all short-read datasets used in this work to the shallowest one (166,282,895 reads). This rarefaction was carried out with a Python script (https://github.com/lyijin/common/blob/master/subsample_fastq.py), prior to read trimming.

The rarefied FASTQ files were processed using a self-written pipeline in Snakemake (https://github.com/lyijin/bismsmark) on the CSIRO High Performance Computing clusters. The pipeline depends on bismark v0.23.1 [44] and trim-galore v.0.6.7 (https://github.com/FelixKrueger/TrimGalore), while automatically applying tool-author-recommended command-line flags to deal with quirks associated with each method (https://github.com/FelixKrueger/Bismark/tree/master/Docs). For the WGBS datasets, the low-complexity bases added by the adaptase in the Accel-NGS kit necessitates the trimming of 10 bp from both ends of R1 and R2, and another 5 bp from off the 5’ end of R2 (i.e., 15 bp in total on the 5’ end of R2). For the EM-seq datasets, as the method-specific flags were added during manuscript preparation, they were treated as normal bisulphite-seq data: no special flags during trimming, but during extraction of methylation levels, methylation information in the first two bases on the 5’ end of R2 were discarded.

Genome-wide methylation levels were obtained by mapping the data (with default parameters) against the human GRCh38 patch 13 genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39/) with ALT and unplaced contigs removed. This was chosen as the updated Infinium MethylationEPIC array annotations (next section) is based on this version, allowing for comparisons that are free of annotation differences.

Loci-specific methylation levels were mapped against KY962518, a more modern 45S reference sequence produced with single-molecule sequencing [45], instead of U13369 that was pieced together from Sanger sequencing data from multiple labs [46]. The full KY962518 sequence contained a~13 kb transcribed region and a~32 kb intergenic spacer. As we were interested in methylation in the former region, we modified the sequence by placing the last 1 kb of the intergenic spacer (putative promoter region) in front of the ~13 kb transcribed region, and discarded the remaining~31 kb of intergenic spacer. To deal with the slightly heterogenous rDNA reads arising from 45S genes that are not fully identical [47], we relaxed mapping parameters with --score-min L,0,-0.6 (bismark default is L,0,-0.2), which allowed reads with more mismatches to map to the 45S reference sequence. This resulted in higher coverages across the locus, and in most CpG dinucleotides, the methylation levels of cytosines on the Watson strand was closer to that on the Crick strand (more concordant methylation levels on both strands; data not shown).

Code written to parse bismark intermediate files into tabular form is at https://github.com/lyijin/cpgberus/tree/master/04_parse_bismark_covs, and code for EM-seq to WGBS statistical comparisons analysed using bsseq (version 1.22.0) [48] and DSS R packages (version 2.34.0) [49] is at https://github.com/lyijin/cpgberus/tree/master/05_CpG_sequence_context with no smoothing or coverage cut-offs. The DSS package models count methylation data using a Beta-Binomial distribution, where biological variations are represented by the dispersion parameter. These dispersion parameters are estimated using a shrinkage estimator based on a Bayesian hierarchical model. Subsequently, a Wald test is conducted at each CpG site [49]. Code that investigated potential MCGW-driven biases in EM-seq readouts is at https://github.com/lyijin/cpgberus/tree/master/13_check_mcgw_emseq_wgbs. Figures 1 and ​and22 were plotted using ggplot2 (version 3.3.3) [50], GGally (version 2.1.1) https://github.com/ggobi/ggally, ComplexHeatmap (version 2.2.0) [51] and motifStack (version 1.30.0) [52].

Genome-wide methylation profiling with Illumina MethylationEPIC arrays

High molecular weight DNA (500–1,000 ng) was sent to the Australian Genome Research Facility (AGRF), Melbourne, Australia. DNA was subjected to bisulphite conversion, and methylation levels of over 850,000 sites were assayed with the Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA).

Following receipt of data, the methylation levels associated with the four samples “WR025V1”, “WR025V9”, “WR069V1” and “WR069V9” were extracted using code documented in https://github.com/lyijin/cpgberus/tree/master/02_process_methepic_data. Raw data was subjected to noob correction [53] and then the beta values extracted and annotated with the Illumina manifest v1.0 B5 GRCh38 genome positions to maintain compatibility with WGBS and EM-seq read mapping. Some array features were discarded. These included 990 probes noted as high variability after a manufacturing change, 38 probes that did not have a GRCh38 genome location, 1,407 probes missing from the v1.0 B5 manifest and 36 pairs of probes mapping to the same GRCh38 coordinates. Illumina also supplied the genome coordinates as 0-based, so these were adjusted to 1-based coordinates.

Three-way analysis of EPIC, EM-seq and WGBS data

The nature of the readouts was a key consideration in this three-way analysis: methylation readouts from EPIC were more “analogue”; while EM-seq and WGBS were more “digital” (e.g., methylation beta values of 0.50 are more common for the short-read methods than EPIC, as this results from having equal numbers of methylated and unmethylated reads). To reduce this “digital” effect, we picked positions that were covered≥5 times in 3 of 4 EM-seq samples and similarly≥5 times in 3 of 4 WGBS samples. This cut-off was stringent enough to minimise the “digital” effect, and fits with observations where gains in sensitivity is greatest going from 1× to 5× [54] yet lenient enough in not forcing all 4 samples to require the minimum 5× coverage that allowed for more positions for downstream analysis. Coverage values in WGBS samples were consistently lower than those for EM-seq: 27.6 million positions were covered≥5 times in 3 of 4 EM-seq samples, while 12.4 million positions were covered≥5 times in 3 of 4 WGBS samples. When further intersected with the >850,000 positions from EPIC, we ended up with a common set of n=103,670 positions with sufficient coverage for all downstream analysis. For these positions, the overall mean coverages from both short-read methods were comparable: 7.59 for EM-seq, 6.70 for WGBS. Our analysis proceeded with the assumption that this gap in coverage did not overly influence the per-position mean betas.

Code written to perform the three-way analysis of EPIC, EM-seq and WGBS data are available at https://github.com/lyijin/cpgberus/tree/master/14_methepic_vs_emseq_wgbs. Code used to annotate genomic context of probed positions in MethylationEPIC arrays utilised parsed databases generated with code at https://github.com/lyijin/cpgberus/tree/master/01_txdb.

We thank Jane Bowen and colleagues for sample collection; Thu Ho for assistance in DNA extraction; the Australian Genome Research Facility (Melbourne, Australia) for conducting the microarray experiments; Ramaciotti Centre for Genomics (UNSW Sydney, Australia) for short-read library sequencing; and Tonia Russell for carrying out the ONT Cas9 work.