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

Decision Tree by Scenario

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.