bio-methylation-dmr-detection

name	bio-methylation-dmr-detection
description	Detects differentially methylated regions (DMRs) from short-read bisulfite (WGBS/RRBS), array, and long-read methylation count tables using dmrseq (permutation region-FDR over the region selection), DSS callDMR (beta-binomial), methylKit tiles, bsseq BSmooth, DMRcate Gaussian-kernel smoothing, metilene, and comb-p. Covers why a DMR is DEFINED by arbitrary thresholds (min-CpGs, max-gap, delta-beta, q) and a smoothing bandwidth, why selecting extreme runs of CpGs then testing them on the same data is post-selection inference, why region q-values are not comparable across tools, and a single-sample domain-segmentation section (PMD, UMR/LMR, MethylSeekR, solo-WCGW) that must run before focal calling on cancer/aging genomes. Use when calling region-level methylation differences, choosing a DMR caller, controlling region-level FDR, or segmenting megabase methylation domains. For per-site testing see differential-cpg-testing; for the methylKit object model see methylkit-analysis.
tool_type	r
primary_tool	dmrseq

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  dmrseqbs testCovariate   cutoff  
sig  dmrsdmrsqval     

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   minlen   minCG  
                dis.merge   pct.sig     

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      

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     

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