| name | bio-applied-dmr-analysis |
| description | Call DMRs from WGBS/RRBS beta values via BSmooth smoothing/t-stats or DSS/methylKit (R); annotate to promoters/CpG islands, correlate with RNA-seq log2FC. Use for DMR calling, DSS callDMR, or methylation-expression integration. |
| tool_type | r |
| primary_tool | DSS |
Differentially Methylated Regions (DMRs)
When to Use
- You have per-CpG methylation calls (Bismark cytosine report / methylKit object) for case vs control and need region-level calls, not single-CpG DMPs.
- You need to run or interpret
DSS::callDMR, dmrseq, or methylKit::calculateDiffMeth output.
- You need to annotate DMRs to promoters, CpG islands, or regulatory elements to judge functional relevance.
- You want to validate DMRs by correlating promoter methylation change with matched RNA-seq log2FC.
- You need a heatmap of beta values across samples for the top DMRs.
Version Compatibility
- R / Bioconductor 3.18–3.19: DSS ≥2.50, methylKit ≥1.28, dmrseq ≥1.26, bsseq ≥1.38, genomation ≥1.34, annotatr ≥1.28 (R ≥4.3)
- Python: numpy ≥1.24, pandas ≥2.0, scipy ≥1.11, pyranges ≥0.0.129
Prerequisites
- Per-CpG methylation counts (methylated/coverage) from Bismark or methylKit — see
bio-methylation-analysis-bismark-alignment, bio-methylation-analysis-methylation-calling.
- Matched RNA-seq differential expression (log2FC) if doing expression integration.
- Basic familiarity with beta-binomial models and multiple-testing correction.
DMR Calling with DSS (R)
Goal: call genome-wide DMRs from WGBS counts while modeling biological overdispersion.
Approach: build BSseq objects per group, run DMLtest with smoothing, then callDMR on the DML statistics.
library(DSS)
library(bsseq)
bs_ctrl <- makeBSseqData(list(ctrl1, ctrl2, ctrl3), sampleNames = c("C1", "C2", "C3"))
bs_treat <- makeBSseqData(list(trt1, trt2, trt3), sampleNames = c("T1", "T2", "T3"))
dml_test <- DMLtest(bs_ctrl, bs_treat, smoothing = TRUE, smoothing.span = 500)
dmrs <- callDMRdml_test p.threshold delta minlen minCG
headdmrsorderdmrsareaStat
Key parameters: smoothing.span (bp bandwidth, 200–500 typical), p.threshold (per-CpG Wald p for DMR seeding), delta (minimum mean methylation difference — avoid calling DMRs with negligible effect size), minlen/minCG (minimum DMR length/CpG count). DSS fits a beta-binomial model per CpG (Y_i ~ Binomial(n_i, p_i), logit(p_i) = mu + epsilon_i) so overdispersion across replicates is captured, unlike a naive binomial test. Because genome-wide DMR calling involves ~25M CpGs, per-CpG p-values aren't directly usable for FDR on regions — DSS/dmrseq instead rank DMRs by the area statistic (sum of per-CpG test statistics across the region) and estimate FDR via permutation of condition labels.
DMR Calling from Beta Values (Python)
Goal: reproduce the BSmooth-style approach when you only have a beta-value matrix (e.g., array/methylKit export) instead of raw counts.
Approach: smooth each replicate across neighboring CpGs, compute a per-CpG t-statistic on the smoothed means, then call DMRs as contiguous runs exceeding a threshold.
import numpy as np
import pandas as pd
from scipy.ndimage import uniform_filter1d
def simulate_region_betas(n_cpgs=400, n_ctrl=3, n_treat=3, seed=1):
"""Simulate per-CpG beta values for a region containing a hyper- and a hypo-DMR."""
rng = np.random.default_rng(seed)
positions = np.sort(rng.choice(np.arange(0, n_cpgs * 30, 30), n_cpgs, replace=False))
baseline = rng.beta(2, 2, n_cpgs)
ctrl_betas = np.clip(baseline + rng.normal(0, 0.05, (n_ctrl, n_cpgs)), 0, 1)
delta = np.zeros(n_cpgs)
delta[80:120] = 0.40
delta[250:290] = -0.35
treat_betas = np.clip(baseline + delta + rng.normal(0, 0.05, (n_treat, n_cpgs)), 0, 1)
return positions, ctrl_betas, treat_betas
def smooth_betas(beta_matrix, window=11):
"""BSmooth-style local smoothing across the CpG axis (borrows strength from neighbours)."""
return np.array([uniform_filter1d(row, size=window, mode="nearest") for row in beta_matrix])
():
dmrs, in_dmr, start = [], ,
i, t (t_stats):
(t) >= threshold:
in_dmr:
in_dmr, start = , i
in_dmr:
length = i - start
length >= min_cpgs:
dmrs.append({: start, : i - , : length,
: t_stats[start:i].()})
in_dmr =
in_dmr ((t_stats) - start) >= min_cpgs:
dmrs.append({: start, : (t_stats) - ,
: (t_stats) - start, : t_stats[start:].()})
pd.DataFrame(dmrs)
positions, ctrl_betas, treat_betas = simulate_region_betas()
ctrl_smooth, treat_smooth = smooth_betas(ctrl_betas), smooth_betas(treat_betas)
mean_ctrl, mean_treat = ctrl_smooth.mean(), treat_smooth.mean()
delta_beta = mean_treat - mean_ctrl
pooled_se = np.sqrt(
ctrl_smooth.var(, ddof=) / ctrl_smooth.shape[]
+ treat_smooth.var(, ddof=) / treat_smooth.shape[]
+
)
t_stat = delta_beta / pooled_se
dmrs = call_dmrs(t_stat, threshold=, min_cpgs=)
dmrs[], dmrs[] = positions[dmrs[]], positions[dmrs[]]
(dmrs[[, , , ]])
DMR Annotation and Expression Integration (Python)
Goal: annotate DMRs to genomic features (promoters, CpG islands) and check whether hypermethylated promoter DMRs correlate with reduced expression.
Approach: interval-join DMRs against a feature table with pyranges, then Pearson-correlate promoter Δβ against matched RNA-seq log2FC (expect a negative slope — hypermethylation silences).
import pyranges as pr
from scipy import stats
def annotate_dmrs_to_features(dmr_df, feature_df):
"""
Left-join DMRs onto genomic features (promoters, CGI/shore/shelf, enhancers).
dmr_df / feature_df need columns: chrom, start, end (+ any extra metadata columns).
"""
dmr_gr = pr.PyRanges(dmr_df.rename(columns={"chrom": "Chromosome", "start": "Start", "end": "End"}))
feat_gr = pr.PyRanges(feature_df.rename(columns={"chrom": "Chromosome", "start": "Start", "end": "End"}))
return dmr_gr.join(feat_gr, how="left").df
def correlate_methylation_expression(promo_delta_beta, rna_log2fc):
"""Pearson r between promoter Delta-beta (methylation change) and RNA-seq log2FC.
Expect r < 0: hypermethylated promoters -> downregulated genes (epigenetic silencing)."""
r, p_val = stats.pearsonr(promo_delta_beta, rna_log2fc)
return r, p_val
Promoter methylation silences transcription because CpG-island methylation blocks TF binding and recruits MBD proteins (MeCP2, MBD1), which in turn recruit HDACs — methylation → deacetylation → compact chromatin → no transcription. Note gene-body methylation is positively correlated with expression (opposite direction from promoters), so only annotate/correlate at the promoter (TSS ± 2 kb).
Pitfalls
- Coordinate systems: BED is 0-based half-open; VCF/GFF/most R packages are 1-based inclusive — mixing them causes off-by-one DMR boundaries.
- Batch effects: check for batch confounding (bisulfite conversion batch, array chip) before interpreting biological signal.
- Multiple testing: never apply per-CpG FDR to ~25M genome-wide tests directly interpreted as region significance — use the DMR area statistic + permutation FDR (DSS/dmrseq), not naive BH on single CpGs.
- Delta-beta thresholding: a statistically significant DMR with
delta near 0 is rarely biologically meaningful — always filter on both p-value/area-stat and minimum mean methylation difference.
- Gene body vs promoter methylation: they correlate with expression in opposite directions; don't average methylation across a whole gene when the question is promoter silencing.
- Low coverage: below ~10x WGBS coverage, per-CpG betas are noisy — rely on smoothing (BSmooth/DSS) rather than single-CpG binomial tests.
See Also
bio-methylation-analysis-dmr-detectionbio-methylation-analysis-methylkit-analysisbio-methylation-analysis-bismark-alignmentbio-differential-expression-de-results
