# Bio Applied Dmr Analysis

> 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.

- Skill: `pavel-kravchenko/bio-applied-dmr-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-dmr-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-dmr-analysis/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Integrations & APIs
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```r
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.

```python
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).

```python
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`