Benchmarking edgeR and methylKit for the Detection of Differential DNA Methylation: A Methodological Evaluation

TL;DR

This study compares two tools, edgeR and methylKit, for detecting DNA methylation changes in mouse embryonic stem cells after morphine exposure.

Contribution

The study evaluates and compares the performance of edgeR and methylKit in detecting differential DNA methylation in a morphine exposure model.

Findings

Both edgeR and methylKit detected global hypomethylation patterns despite differences in the number of DMCs identified.

DMCs were predominantly found in intergenic, intronic, and open sea genomic regions.

High concordance at the gene level (~90%) suggests overlapping sets of differentially methylated genes despite tool differences.

Abstract

Despite the improvements in tool development for DNA methylation analysis, there is a lack of a consensus on computational and statistical models used for differentially methylated cytosine (DMC) identification. This variability complicates the interpretation of findings and raises concerns about the reproducibility and biological significance of the detected results. In this regard, here we conducted a comparative evaluation of edgeR and methylKit tools to assess their performance, concordance, and biological relevance in detecting DMCs following a morphine exposure model in mouse embryonic stem cells (mESCs). Both pipelines were applied to the same WGBS dataset (GEO accession number: GSE292082), and concordance was calculated at both single-base and gene levels. Although the total number of DMCs identified differed between tools, both pipelines detected a global hypomethylation pattern. Genomic distribution analysis revealed that DMCs predominantly localized to intergenic and intronic regions, as well as to open sea regions. Despite differences in sensitivity, both pipelines demonstrated moderate concordance at the DMC level (~56%) and high concordance at the gene level (~90%), identifying largely overlapping sets of differentially methylated genes (DMGs). Comparative assessments further showed that the choice of statistical metric can influence the perceived magnitude of biological effects. Sensitivity analyses indicated that threshold selection and normalization methods influence DMC detection, whereas aggregation at gene level reduces discrepancies. Overall, our findings underscore the complementary strengths of methylKit and edgeR and highlight the importance of careful tool selection for epigenetic studies. As a conclusion, we recommend integrating both pipelines to ensure a balanced interpretation of effect sizes, particularly in studies with complex experimental designs.

Figures

Full text

1. Introduction

Epigenetics is the study of heritable changes in gene expression that do not involve alterations to the underlying DNA sequence. The primary mechanisms of epigenetic regulation include DNA methylation, histone modifications, chromatin remodelling, RNA-based regulation (involving both small and long non-coding RNA) and the emerging role of non-canonical DNA structures [1,2,3]. These modifications are essential for cellular differentiation and function, as well as for normal development, aging and a number of diseases [1,2].

Among previously mentioned mechanisms, DNA methylation is the most stable type of epigenetic modification that plays a significant role in regulating gene expression, maintaining genomic integrity and influencing developmental processes and disease in a cell-type specific manner [4,5,6]. Specifically, this modification primarily consists of the addition of a methyl group to the fifth carbon position of cytosine (5mC) [7,8]. While 5mC is the most prevalent form in vertebrates, it is increasingly recognized that methylation can also occur on other bases, such as N6-methyladenine (6mA), which plays important regulatory roles despite being less frequent and less studied in eukaryotes [9]. It is worth noting that DNA methylation exerts diverse regulatory effects depending on its genomic context since it is not uniformly distributed across the genome. While hypermethylation of CpG-rich promoters is linked to gene silencing in vertebrates (implicated in processes such as genomic imprinting and X-chromosome inactivation), global hypomethylation is associated with genomic instability and aberrant gene activation, contributing to disease development, such as tumorigenesis, neurodegenerative diseases and neurological disorders [10,11,12].

With the fast development of next-generation sequencing technologies, currently the most reliable and widely used method for measuring DNA methylation is bisulfite sequencing (BS-Seq), as it enables a detailed exploration of the epigenome at single-nucleotide resolution [13,14]. In particular, whole-genome bisulfite sequencing (WGBS) has become the gold standard for profiling cytosine methylation as it provides a high-resolution map of DNA methylation across the entire genome, allowing for the identification of different methylation patterns genome-wide [15]. In this technique, DNA is treated with sodium bisulfite, converting unmethylated cytosine residues to uracil (and ultimately to thymine when PCR amplification occurs), whilst maintaining methylated cytosines unchanged and finally providing single-base, quantitative cytosine methylation levels [16]. However, the analysis of methylation sequencing data presents several computational challenges, including the increasing diversity of methodologies and library preparation protocols for WGBS, the complexity of bisulfite-converted data and the need for precise identification of differential methylated cytosines (DMCs) and differential methylated regions (DMRs) [17,18,19].

The first step of analyzing bisulfite-sequencing data is to align reads to a reference genome and to count the number of C-to-T conversions, using read mapping and methylation-calling aligner tools, such as Bismark [14], Methylcoder [20], BRAT [21], RMAP [22], BSSeeker [23,24] and BSMAP [25], among others. Additionally, within the past few years, multiple computational approaches have been developed for downstream DNA methylation studies, which often involve detecting DMCs and DMRs [26]. Some of them implement linear regression designs, including methylKit [27], RnBeads [28] and BSmooth [29]. Other commonly used software are based in binomial model, for example, edgeR [30], DSS [31], BiSeq [32], MOABS [19], RADMeth [33], Bisulfighter [34], methylSig [35], DMRfinder [36], and HMM-DM [37].

Despite the improvements in tool development, a major challenge in DNA methylation analysis arises from the lack of a consensus on statistical and computational tools used for DNA methylation analysis. Indeed, there is no standard protocol to carry out this process [38,39,40]. Specifically, given that these methods rely on distinct statistical models and normalization and dispersion estimation strategies, they lead to variable DMC or DMR outputs with no large percentage of agreement, even when applied to the same dataset. This variability complicates the interpretation of findings and raises concerns about the reproducibility and biological significance of the detected regions.

To address this gap and in order to have a deeper understanding of DMCs identification, in this article, we performed a methodological comparison between two widely used tools to test predefined genomic regions for differential methylation in BS-seq data: methylKit [27] and edgeR [30,41]. Both are available within the statistical packages provided by the Bioconductor software on the R platform and have shown high sensitivity and/or low false-positive rates. However, while edgeR was initially designed for RNA-seq and later adapted for NGS-based methylation data, methylKit was developed specifically for bisulfite sequencing applications. Both tools support differential methylation analysis, even if using different statistical and conceptual frameworks. methylKit [27] employs Fisher’s exact test or logistic regression with overdispersion correction to calculate p-values adjusted to q-values for multiple test correction using a SLIM approach [42]. In comparison, edgeR [30,41] is based on the linear model framework (GLM), quasi-likelihood F tests (QLF) or an exact test that is conceptually similar to Fisher’s exact test, but tailored for overdispersion [43], and models the variation between biological replicates through the negative binomial dispersion supporting both competitive and self-contained tests. This difference in modelling also extends to the reporting metrics. methylKit provides the percentage of methylation change—a highly intuitive metric for biological interpretation—whereas edgeR reports log-fold changes, requiring an extra step of contextualization to relate it to absolute methylation levels. It is worth mentioning that a key divergence was identified in their regional analysis capabilities. methylKit inherently supports de novo DMR calling through genome tiling windows, whereas edgeR’s capacity for regional analysis is generally dependent on predefined genomic intervals. To avoid systematic and base-calling errors, both tools implement different quality control and preprocessing strategies. Specifically, edgeR relies on coverage normalization and removes low-coverage regions to ensure accurate modelling of count data, whereas methylKit applies library-size normalization and uses similar criteria for filtering and removing potential duplication bias, enhancing the reliability of downstream statistical testing.

In this study, a publicly available dataset (GEO accession number: GSE292082) was used as a benchmark to quantitatively evaluate and compare the performance of edgeR and methylKit under controlled and reproducible conditions to share a new perspective on the use of an appropriate tool or tool combination for the identification of DMCs and differentially methylated genes (DMGs). While the dataset used was originally generated in a previous study from our group [44], where the genome-wide effect of chronic morphine exposure on DNA methylation in mESCs was described, the present work focuses specifically on methodological aspects. Here, we aim to rigorously assess and compare the analytical capabilities of these pipelines for detecting DMCs and DMGs, contributing to the reproducibility and transparency of DNA methylation research. Our findings highlight the importance of analytical choices in methylome interpretation and provide a reproducible benchmark framework that can guide researchers in selecting appropriate tools for differential methylation analysis. In fact, it is expected that the combination of different advanced analytical tools will expand the frontiers of methylome analysis, with increasingly relevant applications in environmental epigenetic, biomedical research and precision medicine [45].

2. Results

2.1. Comparative Assessment of DMC Detection Tools

First, to evaluate their suitability for our dataset, a systematic comparative analysis of their methodological frameworks was performed on both R packages, edgeR and methylKit. Supplementary Table S1 summarizes the main features of each tool used for analysing DMCs, including platform, source and implementation language, tool origin, alignment ability, single- or paired-end capacity, major function, DMR calling capabilities, concept, model or test, statistical method, reporting metrics, quality control and preprocessing step, smoothing and further analysis options. Although both packages are compatible with major sequencing protocols and platforms, they diverge in their statistical modelling and reporting metrics (absolute percentages in methylKit vs. log-fold changes in edgeR). After confirming that both tools implement robust but distinct quality control and normalization strategies to mitigate technical biases, their performance in detecting DMCs was compared.

2.2. Similarity Analysis of DMC Detection Tools

To evaluate the concordance between edgeR and methylKit in the detection of DMCs, both tools were applied to DNA methylation data derived from control and morphine-treated mESC samples. edgeR identified a total of 203,337 DMCs, while methylKit detected 223,280 DMCs. To further characterize these findings in depth, changing trends in DMCs were identified. Of these, edgeR classified 60,704 DMCs as hypermethylated (29.85%) and 142,633 as hypomethylated (70.15%). Similarly, methylKit identified 67,054 hypermethylated DMCs (30.03%) and 156,226 hypomethylated DMCs (69.97%) (Table 1 and Figure 1A). Then the overlap between DMCs detected by both tools was identified, amounting to 153,394 DMCs, corresponding to 56.14% of the total unique DMCs identified by either method (n = 273,223), calculated as the union of both sets to avoid double-counting overlapping sites. Within this overlapping set, 44,641 common DMCs (29.10%) were classified as hypermethylated, while 108,753 common DMCs (70.90%) were identified as hypomethylated. This was equivalent to a 53.71% and a 57.20% fraction of the total unique hypermethylated (n = 83,117) and hypomethylated (n = 190,106) DMCs, respectively (Table 1 and Figure 1A). To further evaluate the technical concordance and ensure the detection was not skewed towards a specific direction, a cross-platform comparison and a Spearman correlation analysis was performed between the methylation metrics of both tools (Supplementary Figure S2). For the common DMCs, a moderate-to-strong positive correlation was observed (ρ = 0.62, p < 0.0001), confirming high directional consistency. Furthermore, non-common DMCs (sites detected by only one tool) exhibited a weaker correlation (ρ = 0.49, p < 0.0001), indicating that even for tool-specific sites, the underlying biological signal is preserved. Additionally, the mean and the standard deviation were calculated for both tools to further estimate the coefficient of variation. This value was less than 5% (0.05) in all cases, indicating low relative variability. In order to focus on the commonly identified sites, the distribution of overlapping DMCs was analysed for each chromosome. To explore the chromosomal distribution of commonly identified DMCs, their frequency across chromosomes was analysed. The majority of the DMCs were significantly located on chromosomes 8 and 11. In contrast, chromosomes 2, 3, 4, 7, 12, 13, 18, 19 and sex chromosomes (X and Y) exhibited a notably lower number of DMCs (Supplementary Figure S3). However, it is worth mentioning that a consistent pattern was observed across all chromosomes, and the number of hypomethylated cytosines was always higher than hypermethylated ones.

Given the functional effect of DNA methylation on gene expression, the next step was to examine each tool’s capacity to detect DMCs across specific gene-associated features. Particular attention was given to DMCs located in promoter regions due to their central role in transcriptional regulation. The following gene features were included in the analysis: promoters (≤1 kb, 1–2 kb and 2–3 kb upstream of transcription start sites), 5′UTR, 3′UTR, exons, introns, downstream regions and distal intergenic regions (Figure 1B). In the promoter region edgeR identified a total of 34,817 DMCs (14,015 DMCs within ≤1 kb, 11,216 DMCs in 1–2 kb and 9586 DMCs in 2–3 kb), while methylKit detected a total of 37,699 DMCs (14,974 DMCs within ≤1 kb, 12,140 DMCs in 1–2 kb and 10,585 DMCs in 2–3 kb). The identified number of 5′UTR and 3′UTR DMCs were very similar in both tools, specifically 129 and 4460 DMCs in edgeR and 149 and 4920 in methylKit. Furthermore, for exonic and intronic regions, edgeR reported 9146 and 83,043 DMCs, whereas methylKit identified 10,024 and 92,460 DMCs, respectively. Similarly, for downstream and distal intergenic regions, edgeR identified 2284 and 69,458 DMCs, while methylKit reported 2540 and 75,488 DMCs respectively. As shown in Figure 1B, most DMCs were located in introns and distal intergenic regions, corresponding to 40.84% and 34.16% of gene features in edgeR, and conforming 41.41% and 33.81% in methylKit. Promoter regions and adjacent areas collectively accounted for approximately 17% of DMCs in both tools. In any case, the distribution of DMCs across genomic features was consistent between tools, with only minor variations (Table 2 and Figure 1B). Then, to assess overlap between both tools, common DMCs within specific gene features were identified. For the promoter region a total of 26,031 DMCs were identified (10,324 DMCs within ≤1 kb, 8412 DMCs in 1–2 kb and 7295 DMCs in 2–3 kb); 100 DMCs in 5′UTR; 3403 DMCs in 3′UTR; 6909 DMCs in exons; 63,657DMCs in introns; 1699 DMCs in downstream regions; and 51,595 DMCs in distal intergenic regions. In all cases, the fraction of overlapped DMCs remained close to 56%, in the same way as observed in the general dataset. Additionally, the coefficient of variation remained below 5% (0.05) in the majority of cases, reinforcing the reproducibility of the findings (Table 2).

DNA methylation is not uniformly distributed across the genome. Instead, CpG sites tend to cluster in distinct genomic contexts, including CpG islands (CGIs), shores, shelves and open sea regions (Figure 1C), each with potentially different genetic regulatory roles. Consequently, the way forward was to evaluate the ability of each tool to identify DMCs within these CpG-related features. Both tools predominantly identified DMCs in the area known as open sea (over 90% of all DMCs); that is, edgeR detected 184,129 DMCs, and methylKit identified 203,048 DMCs outside the CpG islands area. Likewise, edgeR identified 985 DMCs in CGI, 9638 DMCs in shores and 8585 DMCs in shelves (summing up to 9.44%), while methylKit detected 849 DMCs in CGI, 9888 DMCs in shores and 9495 DMCs in shelves (amounting to 9.06%). In fact, the 1% of DMCs located in CGIs that simultaneously coincided with promoter regions ranged from 25.79% in edgeR to 22.73% in methylKit (Table 3 and Figure 1C). Next, overlap analysis revealed 566 common DMCs in CGIs, 6839 DMCs in shores, 6604 DMCs in shelves and 139,385 in open sea regions. As with previous comparisons, the proportion of overlapping DMCs remained close to 56% across all categories, except for CGIs, which exhibited a lower overlap of 44.64%. A coefficient of variation under 5% (0.05) was observed in the vast majority of cases, suggesting again low variability between replicates (Table 3).

2.3. Similarity Analysis of DMG Detection Tools

When identifying DMCs, it is crucial that a tool also facilitates the identification of DMGs, as the latter provide a more biologically meaningful and functionally relevant perspective. While DMCs offer site-specific information on methylation changes, DMGs enable the contextualization of these modifications at the gene level, linking methylation patterns to biological processes and transcriptional effect, thereby facilitating interpretation within a physiological and pathological context. Therefore, a comparative evaluation of both tools in DMG detection was conducted. edgeR identified a total of 17,657 DMGs, whereas methylKit reported 17,772 DMGs. Within these, edgeR classified 13,128 DMGs as hypermethylated, 16,313 DMGs as hypomethylated and 11,783 DMGs with both alterations (66.7%). Similarly, methylKit identified 13,426 hypermethylated DMGs, 16,429 hypomethylated DMGs and 12,081 with both marks (68%) (Table 4 and Figure 2A). Analysis of the overlap revealed 16,357 DMGs detected by both tools, corresponding to 87.8% of the total unique DMGs identified by either method (n = 19,072), calculated as the union of both sets. Of these, 11,594 were commonly classified as hypermethylated DMGs and 14,954 as hypomethylated DMGs, while 10,191 common DMGs displayed both methylation alterations (62.30%). This was equivalent to 80.00%, 86.20% and 74.53% of the total unique DMGs in each respective category (Table 4 and Figure 2A). In all comparisons, the coefficient of variation remained below 5% (0.05), reinforcing the reproducibility of the observations (Table 4).

Significantly, while the general overlap in DMCs between edgeR and methylKit was 56%, this value increased substantially to 87.8% when focusing on the gene level (DMGs). To enhance the biological relevance of the results and minimize false positives, equivalent stringency thresholds were applied, that is, a minimum fold-change difference of 2 in edgeR and minimum methylation difference of 25% in methylKit (with significance set at p ≤ 0.05 in both cases) (Figure 2B,C). These parameters reduced the total number of DMCs to 121,122 in edgeR and 156,050 for methylKit. The resulting intersection of both tools identified 64,138 common DMCs, corresponding to 30.11% of the total unique DMCs (n = 213,034). Despite this reduction, the distribution of methylation trends remained consistent with previously described results since common hypermethylated DMCs were 28.31% and common hypomethylated DMCs were 71.69% (Figure 2B). Focusing on the genetic impact, these thresholds decreased the number of DMGs to 16,362 for edgeR and 16,538 for methylKit, with an overlap of 14,896 DMGs. It is worth mentioning that this overlap was 82.74% of the total unique DMGs (n = 18,004). To measure the biological relevance of the identified DMGs, gene ontology enrichment analyses were performed separately on the DMGs identified by each tool, as well as on the set of DMGs commonly detected by both tools, revealing substantial biological overlap. Both edgeR and methylKit identified processes related to cellular development, metabolic activity, chemical stimulus response, cellular localization and regulation of cell adhesion, reflecting core biological functions across datasets. The common list reinforced these findings by highlighting biological processes such as cellular localization, regulation of metabolic and cellular processes, regulation of response to stimulus, and multicellular organismal processes, thereby emphasizing the robustness and the biological relevance of the shared results. While some differences in the ranking of top GO terms were observed between the individual and common sets (Figure 2C), these results offered complementary insights suitable for integrative interpretation.

3. Discussion

This study provides a comprehensive evaluation and comparison analysis between methylKit and edgeR, two widely used computational tools for detecting DMCs in WGBS data. In brief, both methods proved effective in identifying DMCs across various genomic contexts, showing that both tools can be complementary, although notable differences were observed in the total number detected and the way the results were presented. These findings highlight the importance of considering methodological integration in methylation analyses and underscore the need for standardized protocols to ensure reproducibility and consistency in findings with biological relevance.

Although methylKit and edgeR share similar objectives, they differ considerably in their statistical frameworks and normalization strategies. edgeR employs a generalized linear model based on a negative binomial distribution with empirical Bayes dispersion estimation and TMM (trimmed mean of M-values) normalization (which may reduce false positives) [30], while methylKit implements Fisher’s exact test or logistic regression with overdispersion correction, applying multiple testing correction via the SLIM method [27]. These methodological differences were reflected in how each tool reported DMCs, as methylKit provided percentage change values, whereas edgeR presented fold-change values. However, this did not present any impediment to the inter-tool comparison performed in this study, nor in the establishment of an equivalent threshold between both tools to identify the most relevant methylation changes [27,30]. While the metrics differ in nature, they allowed for a synchronized filtering of the data for downstream analysis. Interestingly, methylKit identified more total DMCs than edgeR (223,280 vs. 203,337 DMCs). On the one hand, this may be potentially related to the correction for overdispersion in normalization, which can lead to a more conservative DMC detection in edgeR. On the other hand, that identification difference may reflect the methylKit design as a methylation-specific tool [27], unlike edgeR, which was initially designed to analyse RNA-seq [41] and later modified to identify methylation changes [30]. These differences support previous findings indicating that the choice of analytical pipeline can significantly affect both the sensitivity and specificity of methylation calls [39,46]. Importantly, the cross-platform comparison revealed a moderate level of concordance (56.14%) in DMC detection, with a coefficient of variation below 5%, indicating a strong technical robustness and reproducibility of results, even with a limited number of replicates. This result aligns with earlier benchmarking studies comparing methylation tools, which have noted similarly limited overlap among DMC detection methods using distinct statistical models [39,46].

This interpretation is supported by the correlation analysis, which revealed that both common and non-common DMCs maintain a consistent directional trend. The fact that non-common DMCs still display a correlation of ρ = 0.49 suggests that many tool-specific calls are not biological contradictions but rather marginal cases where one algorithm is more conservative than the other. Specifically, the observed instances where edgeR reports high effect sizes (logFC) relative to moderate methylation differences in methylKit can be attributed to their mathematical frameworks. While methylKit measures absolute shifts in methylation, edgeR is highly sensitive to relative changes; thus, small absolute increments can manifest as disproportionately large fold changes. This divergence is particularly evident in loci where edgeR reports logFC values exceeding 4 despite methylKit showing only moderate percentage shifts (0–40%). Such a pattern highlights a decoupling between statistical sensitivity and biological magnitude: while edgeR is highly effective at identifying statistically significant relative changes in read distribution, these variations do not necessarily translate into large-scale absolute changes in the methylome. That is why tool integration remains essential to contextualize these mathematical artefacts and distinguish high-sensitivity relative signals from substantial physical changes in DNA methylation [47]. Consequently, while the overlap is moderate, the high directional agreement reinforces the reliability of the detected information. Despite these differences, both tools exhibited a shared methylation directionality trend, with a predominant hypomethylation pattern (approximately 70% of DMCs) not only globally but also across individual chromosomes. The relative absence of differential methylation in these regions is in line with the absence of repressive methylation pattern in undifferentiated mESCs, reflecting the maintenance of the pluripotent state [48]. Altogether, these results further support the notion that both tools provide complementary perspectives and that their combined use is critical for robust biological interpretation in methylome studies.

The distribution of DMCs was largely consistent between both platforms, particularly within intronic and distal intergenic regions, which together represented over 70% of total DMCs. This distribution agrees with previous epigenomic landscape analyses showing that DNA methylation variability is widely distributed throughout the genome (not confined to promoter regions) and may regulate enhancer activity and long-range chromatin interactions [6,49,50,51]. When looking to the rest of the gene features, both tools demonstrated similar distribution, further underscoring their robustness despite differences in statistical methodology. Interestingly, although promoter regions accounted for a smaller proportion of total DMCs in both tools (~17%), this result has mechanistic relevance since promoter methylation is linked to gene silencing [52,53]. This trend was also observed by DMCs distributed across different CpG contexts. Specifically, a minor fraction mapped to CGIs, shores and shelves in both tools, and the vast majority of identified DMCs resided in open sea regions (over 90%). These results reinforce the evidence that open sea methylation may act as a dynamic regulatory element sensitive to epigenetic modulation [54]. Although DMCs within CGIs were less frequent, they displayed a high co-occurrence with promoter regions (approximately 25% of CGI DMCs overlapped with promoters), highlighting their potential regulatory relevance. This is in line with the view that methylation outside of canonical CGIs, particularly in enhancer regions, may provide more sensitive markers of cell state and environmental perturbation [55]. These results suggest that while both tools can detect DMCs in high-density CpG contexts, their primary sensitivity lies in identifying methylation variation in lower-density regions. These results align with studies suggesting that integrative approaches, combining multiple algorithms, can provide the most reliable view of methylation landscapes [38,39,40].

In terms of biological relevance, the most striking evidence for complementarity between the tools was identified from the gene -level analysis. Despite the moderate overlap in DMCs, both edgeR and methylKit identified a high concordance (87.8%) in DMG identification. This indicates that the biological impact of methylation changes is preserved at the gene level, even if the precise cytosine positions differ between tools. Furthermore, gene ontology enrichment analyses of the DMGs underlined processes related to stimulus response, cell adhesion, metabolism and development pathways in each tool separately but also after overlapping results, functions that are known to be affected by opioid exposure in other studies [56,57,58]. Notably, while the core biological signal remained consistent, minor shifts in GO rankings between individual and common sets were observed, reflecting the sensitivity of enrichment statistics to list size [59,60]. However, the consensus set acts as a stringency filter that highlights the most robust biological signals by reducing background noise. Importantly, the observed increase in overlap when analysing DMGs (compared to DMCs) suggests that gene-level aggregation may reduce false positives and increase biological interpretability, which supports its use as a complementary tool in methylation studies [61].

While this study provides valuable insights into the comparative performance of two key tools, several limitations should be acknowledged. First, as noted in our comparative assessment, the choice of tools carries inherent trade-offs. While edgeR is robust for handling biological variance, its origin as an RNA-seq tool means it lacks certain native methylation features found in methylKit, such as the ability to perform de novo DMR calling via tiling windows. Second, although both tools allow for effective data filtering, the mathematical nature of logFC in edgeR can lead to a potential overestimation of the biological effect size, particularly in regions with low baseline methylation. In these cases, a high log-fold change may represent a significant relative shift in read ratios without corresponding to a large-scale absolute change in DNA methylation. Therefore, we emphasize that while edgeR excels in detecting statistical significance, its effect size should be contextualized alongside methylKit’s percentage-based metrics to ensure a realistic biological interpretation. Third, the analysis was restricted to one experimental model (mESCs, single cell type and one treatment condition, morphine treatment), which is why results may differ in other biological contexts or with different exposure durations and doses. Fourth, we focused on DMCs and DMGs rather than DMRs; thereby, future work should consider DMRs, which may provide additional site-level information worth taking into account. Fifth, the inclusion of only two replicates per condition may limit statistical power and increase sensitivity to outliers. Future studies should incorporate more replicates and explore integrative frameworks that combine WGBS with transcriptomic or proteomic data for functional validation.

Taken together, these findings underscore the benefits of integrative and complementary strategies to provide a more comprehensive picture of methylation dynamics. While no single tool can provide a complete and accurate representation of the methylome, integrating results from multiple statistical approaches can strengthen analytical confidence and enable more reliable identification of functionally relevant regions, thus reinforcing biologically relevant interpretation. Future directions should involve benchmarking these and other tools across a broader range of biological systems and experimental designs to establish more standardized protocols and promote reproducibility in methylomic research.

4. Materials and Methods

4.1. DNA Processing and Sequencing Workflow

The used DNA samples and sequencing workflow was previously published by our research group and was generated under the same experimental procedure (GEO accession number: GSE292082) [44]. Briefly, genomic DNA from control and morphine-treated (10 μM, Alcaliber; 24-h treatment period) mouse embryonic stem cells (mESCs; Oct4-GFP cell line/PCEMM08, PrimCells, San Diego, CA, USA) was isolated and extracted, using a classic phenol–chloroform/isoamyl methodology with phenol (P4557, Sigma, Saint Louise, MO, USA), chloroform (CL01981000, Scharlau, Sentmenat, Spain) and isoamyl alcohol (BP1150, Fisher BioReagents, Waltham, MA, USA) and following the manufacturer’s instructions. DNA concentration and purity were measured using the Nanodrop Spectrophotometer ND-1000 (Thermo Fisher Scientific, Wilmingyon, DE, USA) by quantifying the 260/280 absorbance ratio. The amount of extracted DNA from two samples was preserved at −80° prior to library preparation.

Genomic DNA was sonicated using a Soniprep 150 (MSE Ldt, London, UK) to produce fragments of approximately 300 base pairs (bps). These fragments were subsequently denatured and prepared for bisulfite conversion, using the EZ DNA Methylation-lightning Kit (Zymo Research, Irvine, CA, USA) to facilitate downstream methylation analysis. DNA fragments were then subjected to library preparation, employing the KAPA Library Preparation Kit (Roche, Pleasanton, CA, USA) along with xGenTM Methyl UDI-UMI Adapters (IDT, Coralville, IA, USA). The quality and integrity of the libraries were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) with the DNA 7500 assay. Following library preparation, sequencing was performed on an Illumina NovaSeq 6000 S4 platform (paired-end) (Illumina, San Diego, CA, USA), with 4× multiplexing to ensure robust coverage, achieving a minimum of 50,000 reads per replicate for each sample. All cytosine residues across both duplicates for each sample were included in subsequent genome-wide methylation analyses.

4.2. Pre-Processing

The library fragment adapters were trimmed using Trim Galore! (v.0.6.2) [62], and subsequently, a FastQC High-Throughput Sequence Quality Control report (v.0.11.6) [63] was generated to evaluate the quality of the WGBS FASTQ files, indicating a quality score exceeding 30 and foreseeing a suitable mapping efficiency. The resultant reads were concatenated using Cat (v.8.22) [64]. Spearman correlation confirmed the reproducibility of the two biological replicates in methylation samples, demonstrating high concordance (Supplementary Figure S1A).

4.3. Cross-Platform Analysis of m5C Sequencing Data

The WGBS data was mapped to a bisulfite-converted reference genome file for mouse (UCSC GRm38/mm10) with Bowtie2 [65,66] on Bismarck (v.0.22.1) [67] with default parameters, and the resulting binary alignment files were sorted and indexed with SAMtools (v.14.0) [68] to enable quicker access. To be unable to use the post-processing programs that follow the WGBS bioinformatics analysis, a MethylDackel (v.0.5.1) [69] tool was used to estimate the percentage coverage and percentage methylation for all CpG sites on a genome-wide scale. Then, to evaluate the effect of morphine on differentially methylated cytosines (DMC), two different statistical tools were chosen (under R v4.0.3): (1) edgeR (v.3.32.1) [30,41], which uses the negative binomial distribution to assess statistical significance, with TMM normalization accounting for library size and (2) methylKit (v.1.16.1) [27], which normalizes read coverage distribution between samples. In both cases regions with a p-value <  0.05 were considered differentially methylated. The PCA (principal component analysis) in each of the tools confirmed the similarity between replicates and distinguishable differences between samples (control and morphine-treated duplicate samples). After applying TMM normalization (correcting for library size), a read count distribution plot showed a correct data normalization and Volcano plot analysis confirmed the differences in these DMCs analysed with two tools (Supplementary Figure S1B). It is worth mentioning that both tools present the DMCs in different ways: edgeR uses foldchange values, and methylKit applies percentage values. The methylation landscape is viewable via the UCSC genome browser [70]. Integrative analyses between both tools were performed using Venny tools (v.2.1.0) [71] and the Gene Ontology Resource from the GO Consortium (https://geneontology.org/, v.16.1.0) [72,73] was used to identify the biological functions. The computational workflow summarizing all the steps is included in Figure 3.

5. Conclusions

In conclusion, this comparative study demonstrates that while edgeR and methylKit each have unique advantages, their complementary use enhances the robustness and biological interpretability of differential DNA methylation analysis. The convergence in general trends, specifically the predominance of hypomethylation in response to morphine, the regional consistency of DMCs across gene features and the high concordance of overlapped DMGs identification, supports the validity of findings derived from both statistical frameworks. However, we found that tool-specific metrics used for results quantification significantly impact biological interpretation. Recognizing these mathematical differences is essential to avoid overestimating effect sizes. Therefore, the proposed tool integration is necessary to ensure that methodological consistency reinforces the biological relevance of the shared results.

Given the lack of a standardized pipeline in the field, cross-validation across multiple platforms remains essential to improve accuracy in DMC detection and robust and reproducible epigenetic profiling. As DNA methylation continues to gain visibility as a key epigenetic biomarker, computational strategies that combine sensitivity, specificity, and flexibility will be essential. Future studies should continue exploring the impact of pipeline design, normalization strategy, and statistical modelling in large-scale WGBS datasets to enable more comprehensive and biologically meaningful methylome interpretations.