GPTomics/bio-methylation-dmr-detection

Version Compatibility

Reference examples tested with: dmrseq 1.22+, DSS 2.50+, methylKit 1.28+, bsseq 1.38+, DMRcate 2.16+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

The GENOME BUILD is a version that matters. methylKit/bsseq/DMRcate assembly= is metadata, but annotation packages (annotatr::build_annotations(genome='hg38'), TxDb.Hsapiens.UCSC.hg38.knownGene) are build-specific and must match the alignment genome. DMRcate arraytype='EPIC'/'450K' and the IlluminaHumanMethylation annotation package set the array CpG universe. DMRcate defaults shift across Bioconductor releases (the C kernel scaling has no single default - it is platform-dependent); confirm with ?dmrcate on the installed build.

DMR Detection

"Find differentially methylated regions" -> Detect candidate regions, then score them with a null that re-ran the selection - because a DMR is defined by a chain of thresholds, not found, and only a selection-aware q-value is honest.

• R: dmrseq(bs, testCovariate='condition') (selection-aware) ; DSS::callDMR(), methylKit::tileMethylCounts(), bsseq BSmooth, DMRcate::dmrcate() (combine-and-correct)

Scope: REGION-level differential methylation from any per-CpG methylation+coverage table (WGBS/RRBS counts, array beta/M, long-read modkit bedMethyl), plus single-sample domain segmentation. Per-site DMC/DMP testing -> differential-cpg-testing. The methylKit import/filter/unite object model -> methylkit-analysis. Long-read MM/ML calling that produces the counts -> long-read-sequencing/nanopore-methylation. Functional enrichment of DMR genes -> pathway-analysis/go-enrichment (with the CpG-bias correction noted below).

The Single Most Important Modern Insight -- A DMR Is DEFINED, Not Found, and the Region p-Value Is Only Honest If Its Null Re-Ran the Selection

The naive recipe - compute a per-CpG statistic, select runs of CpGs that look extreme to DEFINE candidate regions, then test those same regions and report a p/q on the same data - reuses the data twice. The regions were CHOSEN because they were extreme; scoring them with the data that selected them inflates significance and produces uncalibrated region FDR. This is post-selection inference, the field's original sin. No error is thrown; the q-values just lie. Three corollaries:

1. dmrseq neutralizes the selection; the others combine-and-correct but do not re-select. dmrseq (Korthauer 2019 Biostatistics 20:367) builds the null by PERMUTING the condition labels and RE-RUNNING the entire candidate-detection procedure on each permutation, pooling permuted region statistics into one genome-wide null - so its q-value is on REGIONS and accounts for selection. comb-p and DMRcate combine existing per-CpG p-values (Stouffer/Fisher/SLK) and apply a region multiplicity correction, but never re-run selection. methylKit tiles use FIXED windows (boundaries not chosen from the data, so the selection part is sidestepped) but ignore inter-tile correlation. DSS callDMR merges significant CpGs (threshold-then-define).

2. Region q-values are NOT comparable across tools. "I found 4,000 DMRs at q<0.05" is meaningless without naming the tool and what its q controls. Cross-tool OVERLAP is the real evidence - run two callers and intersect.

3. Thresholds are conventions, not biology. delta-beta 25%, min-CpGs 3, max-gap 1000bp, q<0.01 are tutorial folklore. The same data yields radically different DMR sets under different settings. Report and justify every knob; never present one config as correct.

Organize the analysis around defending these, not around listing functions.

Tool Taxonomy

| Tool | Citation | Mechanism / role | When | |------|----------|------------------|------| | dmrseq | Korthauer 2019 Biostatistics 20:367 | GLS area statistic + permutation null that re-runs selection; smooths the difference internally | headline WGBS inference; calibrated region FDR; >=2 reps/group | | DSS callDMR | Feng 2014 Nucleic Acids Res 42:e69; Park & Wu 2016 Bioinformatics 32:1446 | per-CpG Bayesian beta-binomial dispersion shrinkage, then merge significant CpGs | small n, complex/multi-factor designs; low-coverage with smoothing | | methylKit tiles | Akalin 2012 Genome Biol 13:R87 | fixed windows + logistic/F test per tile | fast RRBS/WGBS screening; reuses the methylKit object model | | bsseq (BSmooth) | Hansen 2012 Genome Biol 13:R83 | per-sample local-likelihood smoothing -> smoothed t-statistic | low/uneven-coverage WGBS; superseded by dmrseq for calibrated FDR | | DMRcate | Peters 2015 Epigenetics Chromatin 8:6; Peters 2021 Nucleic Acids Res 49:e109 | Gaussian-kernel smoothing of the per-CpG statistic | arrays (450K/EPIC) and WGBS (different kernel) | | metilene | Juhling 2016 Genome Res 26:256 | binary segmentation + 2D Kolmogorov-Smirnov; standalone C CLI | fast whole-genome second caller; data-adaptive boundaries; tolerates missingness | | comb-p | Pedersen 2012 Bioinformatics 28:2986 | ACF + Stouffer-Liptak-Kechris combination + Sidak; Python CLI on any p-values | array EWAS regions; tool-agnostic corroboration of any per-CpG p-values |

Decision Tree by Scenario

| Scenario | Recommended | Why | |----------|-------------|-----| | WGBS, >=2 reps/group, headline region inference | dmrseq | only caller whose region FDR accounts for selection | | WGBS, small n, multi-factor / covariates | DSS (DMLtest.multiFactor -> callDMR) | beta-binomial dispersion shrinkage + model formula | | Low / uneven coverage WGBS | dmrseq (smooths internally) or DSS smoothing=TRUE | smoothing borrows strength across CpGs | | RRBS, quick region screen | methylKit tiles (filtered, cov.bases>=3) | fixed windows; fast; CpG-island enriched; screening q, not selection-corrected | | EPIC/450K array | DMRcate array mode (cpg.annotate('array')) or comb-p on limma p | kernel tuned for array spacing; methylKit/bsseq/dmrseq are count-based | | Corroborate any DMR set | run a second caller, intersect | region q is not comparable across tools | | Cancer / aging / placenta / cultured-cell WGBS | segment PMDs FIRST (see domain section) | a focal caller manufactures fake hypo-DMRs from PMD background | | Per-CpG, not region | -> differential-cpg-testing | site-level test before aggregating to regions | | Long-read modkit bedMethyl input | -> long-read-sequencing/nanopore-methylation, then any caller here | modkit calls; feed Nmod + Nvalid_cov here for region statistics |

dmrseq: The Selection-Aware Headline Caller

Goal: Call WGBS DMRs with a region-level FDR that survives the region-selection step.

Approach: Build a BSseq object from raw counts (do NOT pre-smooth), filter loci with zero coverage in any sample, then run dmrseq, which detects candidate regions (runs of CpGs whose smoothed methylation-difference coefficient exceeds cutoff), computes one GLS area statistic per region, and generates the null by permuting labels and re-running detection.

library(dmrseq)
library(bsseq)

bs <- read.bismark(c('ctrl1.cov.gz', 'ctrl2.cov.gz', 'treat1.cov.gz', 'treat2.cov.gz'),
                   colData = DataFrame(condition = c('ctrl', 'ctrl', 'treat', 'treat')),
                   rmZeroCov = TRUE, strandCollapse = TRUE)

# dmrseq requires every locus to have non-zero coverage in EVERY sample.
bs <- bs[rowSums(getCoverage(bs) == 0) == 0, ]

# cutoff=0.1 only SEEDS candidate detection (10% smoothed difference); significance
# comes from the permutation statistic, so dmrseq does NOT hard-threshold delta-beta.
dmrs <- dmrseq(bs, testCovariate = 'condition', cutoff = 0.1)
sig <- dmrs[dmrs$qval < 0.05]   # qval is a REGION FDR that accounts for selection

dmrseq smooths the DIFFERENCE internally (bpSpan/minInSpan/maxGapSmooth); running BSmooth() first double-smooths and invalidates the model. The permutation null is COARSE at 2-vs-2 (few distinct label permutations) - the package pools across all candidate regions to compensate, but more replicates give a finer null. Use adjustCovariate for nuisance variables and block = TRUE for large-scale differential blocks.

DSS: Beta-Binomial Dispersion Shrinkage

Goal: Call DMRs with per-CpG dispersion shrinkage for small n or a multi-factor design.

Approach: DMLtest does per-CpG Wald tests with a Bayesian beta-binomial dispersion estimate, then callDMR merges significant CpGs into regions.

library(DSS)

bs <- makeBSseqData(list(c1, c2, t1, t2), c('C1', 'C2', 'T1', 'T2'))   # each: chr/pos/N/X data.frame
dml <- DMLtest(bs, group1 = c('C1', 'C2'), group2 = c('T1', 'T2'), smoothing = TRUE)   # TRUE for low-cov WGBS

# callDMR defaults (verify on installed build): delta=0, p.threshold=1e-5,
# minlen=50, minCG=3, dis.merge=100, pct.sig=0.5.
# delta=0 means NO effect-size floor - SET it explicitly so tiny shifts are not called.
dmrs <- callDMR(dml, delta = 0.1, p.threshold = 1e-5, minlen = 50, minCG = 3,
                dis.merge = 100, pct.sig = 0.5)   # pct.sig=0.5: >=50% of region CpGs individually significant

callDMR merges significant CpGs (threshold-then-define), so its region p is NOT selection-corrected; its strengths are dispersion shrinkage and multi-factor support (DMLtest.multiFactor with a model formula).

bsseq BSmooth

Goal: Call DMRs on low/uneven-coverage WGBS by smoothing each sample before testing.

Approach: Smooth per sample, compute the smoothed t-statistic, then threshold it into regions.

library(bsseq)

bs_smooth <- BSmooth(bs, BPPARAM = MulticoreParam(4), verbose = TRUE)
keep <- rowSums(getCoverage(bs_smooth) >= 2) == ncol(bs_smooth)   # >=2x in every sample
bs_filt <- bs_smooth[keep, ]

tstat <- BSmooth.tstat(bs_filt, group1 = c('C1', 'C2'), group2 = c('T1', 'T2'),
                       estimate.var = 'same', mc.cores = 4)
dmrs <- dmrFinder(tstat, cutoff = c(-4.6, 4.6))   # t-stat cutoff (Hansen 2012 uses quantile-based cutoffs)

dmrFinder needs BSmooth.tstat output, not the smoothed BSseq object directly. BSmooth gives a RANKED DMR list, not a calibrated region FDR - dmrseq (same lab lineage) supersedes it for region inference. Over-smoothing (too-wide bandwidth) washes out focal promoter DMRs.

DMRcate: The Array-vs-WGBS Fork

Goal: Call DMRs by Gaussian-kernel smoothing of the per-CpG statistic, with the correct kernel for the platform.

Approach: Annotate per-CpG statistics through the array OR sequencing entry point, then smooth and extract.

library(DMRcate)

# ARRAY (450K/EPIC): beta/M matrix; arraytype sets the CpG universe.
design <- model.matrix(~ condition)
ann_array <- cpg.annotate('array', m_values, what = 'M', arraytype = 'EPIC',
                          analysis.type = 'differential', design = design, coef = 2)
dmrs_array <- extractRanges(dmrcate(ann_array, lambda = 1000, C = 2))   # array kernel

# WGBS: different entry point AND a much smaller kernel - the array C=2 over-smooths
# dense sequencing CpGs ~25x. Annotate from a count/edgeR-DSS path, then C=50.
ann_seq <- sequencing.annotate(bs, design = design, coef = 2)
dmrs_seq <- extractRanges(dmrcate(ann_seq, C = 50))   # WGBS kernel (Peters 2021)

The default lambda=1000, C=2 is ARRAY-only; applying it to WGBS produces massively over-smoothed, merged, inflated DMRs. pcutoff='fdr' returns no DMRs if the upstream limma/DSS yields no significant CpGs.

metilene and comb-p (Second Callers)

metilene is a standalone C CLI taking one tab table (chrom, pos, per-sample methylation rate); binary segmentation + a 2D-KS test find data-adaptive boundaries. comb-p is a Python CLI that takes a BED of per-CpG p-values from ANY upstream test, estimates the p-value autocorrelation, does a Stouffer-Liptak-Kechris combination of neighbors, groups regions, and applies a one-step Sidak correction.

metilene -M 1000 -m 10 -d 0.1 -a g1 -b g2 input.tsv | metilene_output.pl   # -m min CpGs, -d min mean diff
comb-p pipeline -c 4 --seed 0.01 --dist 500 --step 50 -p out methyl_pvals.bed   # --seed = p to start a region

Both combine-and-correct rather than re-select; use them as fast corroboration and intersect with dmrseq.

Thresholds Are Conventions; Region FDR Is Tool-Specific

Every caller exposes the same coupled knobs under different names: min-CpGs (minNumRegion/minCG/-m/min.cpgs/cov.bases), max-gap (maxGap/dis.merge/-M), delta-beta (cutoff/delta/-d/betacutoff/difference), and a significance cutoff. Shrinking max-gap, raising min-CpGs, and raising delta all reduce the DMR count, and the same data yields wildly different DMR sets. The phrase "region-level FDR" means three different objects: a selection-aware permutation FDR (dmrseq), a BH/Sidak correction on combined per-CpG p-values (DMRcate/comb-p), or per-unit q on independent tiles/CpGs (methylKit/DSS). Report all knobs, name the tool, and use cross-tool overlap as the evidence statement.

DMR-to-Gene Mapping and the CpG-Density Enrichment Bias

Goal: Interpret DMRs without inflating enrichment from CpG-rich genes.

Approach: Annotate DMRs to features (annotatr returns one row per DMR-feature overlap; genomation collapses by precedence), then run enrichment with a method that corrects for CpG/probe count.

library(annotatr)
annots <- build_annotations(genome = 'hg38', annotations = c('hg38_basicgenes', 'hg38_cpg_islands'))
dmr_ann <- annotate_regions(regions = sig, annotations = annots, ignore.strand = TRUE)   # one row per overlap

# Enrichment: methylation has a CpG-density bias (CpG-rich genes harbor DMRs by chance),
# so a plain hypergeometric GO test is biased. Use missMethyl goregion (probe/CpG-bias-aware).
# missMethyl::goregion(sig_ranges, all.cpg=..., collection='GO', array.type='EPIC')

A single DMR commonly overlaps or sits between several genes; mapping DMR -> gene (nearest TSS vs overlap vs within-X-kb) is a modeling choice that changes the gene list. Hand the corrected enrichment to pathway-analysis/go-enrichment, flagging that methylation input needs a CpG-bias-aware method (missMethyl gometh/goregion), not a generic hypergeometric test.

Single-Sample Domain Structure (NOT Differential)

This is a DIFFERENT problem from the focal between-group callers above. The mammalian methylome partitions at MEGABASE scale into Highly Methylated Domains (HMDs, ~80-90%, ordered) and Partially Methylated Domains (PMDs, ~40-70%, disordered, high-variance), and PMDs coincide with late replication, Lamina-Associated Domains, and the Hi-C B-compartment (Lister 2009 Nature 462:315; Berman 2012 Nat Genet 44:40). Cancer "global hypomethylation" is a DOMAIN phenomenon - focal CpG-island hypermethylation sitting ON a background of megabase PMD hypomethylation - not a focal one. Domain structure is a SINGLE-SAMPLE, structural question answered by SEGMENTERS, not by any between-group DMR caller (methylKit/DSS/dmrseq/DMRcate/metilene/comb-p have no single-sample segmentation mode).

• MethylSeekR (Burger 2013 Nucleic Acids Res 41:e155) segments one WGBS methylome into UMRs (CpG-rich unmethylated = promoters/CGIs), LMRs (CpG-poor low-methylated ~30% = distal enhancers), and PMDs. Pipeline: readMethylome() -> plotAlphaDistributionOneChr() (diagnostic: does the sample have PMDs?) -> segmentPMDs() (2-state Gaussian HMM, 101-CpG windows) -> calculateFDRs() -> segmentUMRsLMRs(m=0.5, n=..., pmdGRanges=...). PMDs MUST be masked before UMR/LMR calling or PMD disorder spawns spurious LMRs.
• solo-WCGW (Zhou 2018 Nat Genet 50:591) - an isolated CpG in [A/T]CG[A/T] context - loses methylation fastest and most monotonically with cell division and is the most sensitive PMD/mitotic-clock readout, detecting PMD hypomethylation even in near-normal tissue. Quantify as the mean over the published common-PMD solo-WCGW CpG set, not as a DMR. See epigenetic-clocks for the broader clock taxonomy.

The warning: running a focal DMR caller on a PMD-bearing genome manufactures thousands of fake hypo-DMRs that are really one phenomenon - the PMD background shifting - chopped into pieces by the max-gap/min-CpG knobs. Segment domains FIRST, then EXCLUDE PMD intervals from focal calling or STRATIFY every DMR by in-PMD vs out-of-PMD and report the fraction that is PMD background.

Per-Method Failure Modes

PMD background reported as DMRs

Trigger: focal caller on tumor/aged/placenta/cultured WGBS without domain screening. Mechanism: megabase PMD hypomethylation chopped into pieces by max-gap/min-CpG. Symptom: thousands of large hypo-DMRs in gene-desert, late-replicating, low-CpG-density coordinates. Fix: segment PMDs (MethylSeekR) first; exclude or stratify; report the PMD fraction.

Pre-smoothing before dmrseq

Trigger: BSmooth() then feeding the smoothed object to dmrseq. Mechanism: dmrseq smooths the difference internally; pre-smoothing double-smooths. Symptom: distorted candidate regions and invalid statistics. Fix: feed dmrseq the raw BSseq counts.

DMRcate array defaults on WGBS

Trigger: copying lambda=1000, C=2 onto sequencing data. Mechanism: the array kernel is ~25x too wide for dense WGBS CpGs. Symptom: massively over-smoothed, merged, inflated DMRs. Fix: sequencing.annotate() + C=50 for WGBS (Peters 2021).

Threshold-then-test reported as region FDR

Trigger: greping runs of significant per-CpG calls and reporting the per-CpG q. Mechanism: the regions were selected for extremeness, then tested on the same data. Symptom: anti-conservative, uncalibrated region q. Fix: use dmrseq (selection-aware) for the headline; at minimum state that a tile/merge q is a screening q.

Single-CpG tiles

Trigger: tileMethylCounts at the default cov.bases=0. Mechanism: a window with one covered CpG becomes a "DMR." Symptom: thousands of single-CpG noisy regions. Fix: set cov.bases >= 3.

Cross-tool count comparison

Trigger: comparing "N DMRs at q<0.05" between callers. Mechanism: each tool's q controls a different object. Symptom: apparent disagreement that is really an FDR-definition mismatch. Fix: compare OVERLAP, not counts.

Quantitative Thresholds

| Threshold | Source | Rationale | |-----------|--------|-----------| | dmrseq cutoff 0.1 (candidate seed only) | Korthauer 2019 | seeds detection; significance is the permutation statistic, NOT a delta floor | | DSS callDMR(delta=0) default -> SET it | Park & Wu 2016; DSS docs | delta=0 calls regions with no effect-size floor; set ~0.1 to require a real shift | | min-CpGs per region 3-5 | convention | single-CpG "regions" are DMPs in disguise; trades sensitivity vs specificity | | delta-beta 25% ("moderate") | methylKit tutorial folklore | NOT derived; biologically meaningful delta is feature- and purity-dependent | | methylKit tileMethylCounts(cov.bases>=3) | nuance (default is 0) | the default 0 lets single-CpG tiles through | | coverage floor ~10x per CpG | field standard | a single-CpG beta below ~10x is a coin flip | | DMRcate WGBS C=50 (array C=2) | Peters 2021 Nucleic Acids Res 49:e109 | dense WGBS CpGs need a far narrower kernel than array probes | | MethylSeekR m=0.5, FDR<5% | Burger 2013 | methylation cutoff for hypomethylated regions; FDR target picks the CpG-count threshold n |

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | dmrseq error about zero-coverage loci | a locus has 0 coverage in some sample | filter rowSums(getCoverage(bs)==0)==0 first | | Thousands of huge hypo-DMRs in cancer WGBS | PMD background not segmented | MethylSeekR segmentPMDs first; exclude/stratify | | Over-merged WGBS DMRs with DMRcate | array kernel on sequencing | sequencing.annotate() + C=50 | | dmrFinder errors on a smoothed object | needs BSmooth.tstat output | run BSmooth.tstat before dmrFinder | | DMRcate returns no DMRs | pcutoff='fdr' and no significant upstream CpGs | check the upstream limma/DSS result first | | Biased GO enrichment of DMR genes | plain hypergeometric ignores CpG density | use missMethyl goregion/gometh |

References

• Korthauer K, Chakraborty S, Benjamini Y, Irizarry RA. 2019. Detection and accurate false discovery rate control of differentially methylated regions from whole genome bisulfite sequencing. Biostatistics 20:367-383.
• Feng H, Conneely KN, Wu H. 2014. A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data. Nucleic Acids Res 42:e69.
• Park Y, Wu H. 2016. Differential methylation analysis for BS-seq data under general experimental design. Bioinformatics 32:1446-1453.
• Hansen KD, Langmead B, Irizarry RA. 2012. BSmooth: from whole genome bisulfite sequencing reads to differentially methylated regions. Genome Biol 13:R83.
• Akalin A, Kormaksson M, Li S, et al. 2012. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol 13:R87.
• Peters TJ, Buckley MJ, Statham AL, et al. 2015. De novo identification of differentially methylated regions in the human genome. Epigenetics Chromatin 8:6.
• Peters TJ, Buckley MJ, Chen Y, et al. 2021. Calling differentially methylated regions from whole genome bisulphite sequencing with DMRcate. Nucleic Acids Res 49:e109.
• Juhling F, Kretzmer H, Bernhart SH, Otto C, Stadler PF, Hoffmann S. 2016. metilene: fast and sensitive calling of differentially methylated regions from bisulfite sequencing data. Genome Res 26:256-262.
• Pedersen BS, Schwartz DA, Yang IV, Kechris KJ. 2012. Comb-p: software for combining, analyzing, grouping and correcting spatially correlated P-values. Bioinformatics 28:2986-2988.
• Lister R, Pelizzola M, Dowen RH, et al. 2009. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature 462:315-322.
• Berman BP, Weisenberger DJ, Aman JF, et al. 2012. Regions of focal DNA hypermethylation and long-range hypomethylation in colorectal cancer coincide with nuclear lamina-associated domains. Nat Genet 44:40-46.
• Zhou W, Dinh HQ, Ramjan Z, et al. 2018. DNA methylation loss in late-replicating domains is linked to mitotic cell division. Nat Genet 50:591-602.
• Burger L, Gaidatzis D, Schubeler D, Stadler MB. 2013. Identification of active regulatory regions from DNA methylation data. Nucleic Acids Res 41:e155.

Related Skills

• differential-cpg-testing - Per-site testing before region aggregation
• methylkit-analysis - methylKit object model and tile construction
• methylation-calling - Produces the input count tables
• array-preprocessing - Array beta/M-value input for DMRcate array mode
• epigenetic-clocks - Mitotic-clock / solo-WCGW overlap (domain section)
• pathway-analysis/go-enrichment - CpG-bias-aware enrichment of DMR genes (missMethyl gometh)
• long-read-sequencing/nanopore-methylation - Pipe modkit bedMethyl counts here for region statistics
• workflows/methylation-pipeline - End-to-end bisulfite pipeline