Bio applied dmr analysis

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)

# ctrl1..3 / trt1..3 are data.frames with columns: chr, pos, N (coverage), X (methylated)
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 (differentially methylated loci) test with local smoothing
dml_test <- DMLtest(bs_ctrl, bs_treat, smoothing = TRUE, smoothing.span = 500)

# Call DMRs from DML results
dmrs <- callDMR(dml_test, p.threshold = 0.001, delta = 0.1, minlen = 50, minCG = 3)
head(dmrs[order(-dmrs$areaStat), ])

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    # hypermethylated DMR
    delta[250:290] = -0.35  # hypomethylated DMR
    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])


def call_dmrs(t_stats, threshold=3.5, min_cpgs=5):
    """Call DMRs as contiguous runs of CpGs with |t| >= threshold spanning >= min_cpgs CpGs."""
    dmrs, in_dmr, start = [], False, None
    for i, t in enumerate(t_stats):
        if abs(t) >= threshold:
            if not in_dmr:
                in_dmr, start = True, i
        elif in_dmr:
            length = i - start
            if length >= min_cpgs:
                dmrs.append({"start_idx": start, "end_idx": i - 1, "n_cpgs": length,
                             "sum_t": t_stats[start:i].sum()})
            in_dmr = False
    if in_dmr and (len(t_stats) - start) >= min_cpgs:  # DMR running off the end of the array
        dmrs.append({"start_idx": start, "end_idx": len(t_stats) - 1,
                     "n_cpgs": len(t_stats) - start, "sum_t": t_stats[start:].sum()})
    return 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(0), treat_smooth.mean(0)
delta_beta = mean_treat - mean_ctrl

pooled_se = np.sqrt(
    ctrl_smooth.var(0, ddof=1) / ctrl_smooth.shape[0]
    + treat_smooth.var(0, ddof=1) / treat_smooth.shape[0]
    + 1e-6
)
t_stat = delta_beta / pooled_se

dmrs = call_dmrs(t_stat, threshold=3.5, min_cpgs=5)
dmrs["start_bp"], dmrs["end_bp"] = positions[dmrs["start_idx"]], positions[dmrs["end_idx"]]
print(dmrs[["start_bp", "end_bp", "n_cpgs", "sum_t"]])

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-detection
• bio-methylation-analysis-methylkit-analysis
• bio-methylation-analysis-bismark-alignment
• bio-differential-expression-de-results