# Dna Methylation

> Align WGBS/RRBS bisulfite reads with Bismark; call DMRs/DMPs with methylKit or BSmooth. Use when analyzing bisulfite sequencing, CpG beta values, .cov/cytosine_report files, or DNAm epigenetic age (Horvath/GrimAge).

- Skill: `pavel-kravchenko/dna-methylation` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/dna-methylation`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/dna-methylation/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/dna-methylation

---


# dna-methylation

## When to Use
- Processing raw WGBS or RRBS FASTQ reads into per-CpG methylation calls (Bismark pipeline)
- Calling differentially methylated regions/positions (DMRs/DMPs) between two or more conditions
- Loading and QC-filtering Bismark coverage (`.cov`) or cytosine report files in Python/R
- Estimating biological age from methylation beta values (Horvath, Hannum, PhenoAge, GrimAge clocks)
- Deciding between array (EPIC/450K) vs WGBS/RRBS coverage, or interpreting non-conversion QC metrics

## Version Compatibility
- Bismark ≥0.24, Bowtie2 ≥2.5, Trim Galore ≥0.6
- methylKit ≥1.28 (Bioconductor 3.18+), R ≥4.3
- pandas ≥2.0, numpy ≥1.24, scipy ≥1.11, Python ≥3.10

## Prerequisites
- `bismark`, `bowtie2`, `trim_galore` on PATH for alignment; a reference genome FASTA
- R packages: `methylKit`, `BSmooth` (from `bsseq`) for DMR calling
- Python: `pandas`, `numpy`, `scipy`, `matplotlib`
- Familiarity with FASTQ/BAM basics (see `bio-read-qc-fastp-workflow`, `bio-read-alignment-bowtie2-alignment`)

## Quick Reference

| Step | Tool | Notes |
|------|------|-------|
| Genome prep | `bismark_genome_preparation` | Creates CT/GA converted index |
| Alignment | `bismark` | Bowtie2 backend |
| Deduplication | `deduplicate_bismark` | Remove PCR duplicates |
| Extraction | `bismark_methylation_extractor` | CpG/CHH/CHG contexts |
| Load in R | `methylKit::methRead()` | From bismark coverage files |
| DMP test | `methylKit::calculateDiffMeth()` | Logistic regression |
| DMR calling | `BSmooth::BSmooth()` | Smoothing + t-statistics |
| Age clock | Horvath coefficient dot product | 353 CpGs, anti-logit transform |

**Goal:** Turn raw bisulfite FASTQ reads into per-CpG methylation percentages.
**Approach:** Bismark converts reads and genome to a 3-letter alphabet (C→T), aligns with Bowtie2, then counts methylated (still C) vs unmethylated (converted to T) calls per cytosine.

```bash
# 1. Prepare bisulfite-converted genome index (one-time, per reference)
bismark_genome_preparation /path/to/genome/

# 2. Trim adapters + low-quality bases (WGBS/RRBS both benefit from this)
trim_galore --paired --fastqc sample_R1.fastq.gz sample_R2.fastq.gz

# 3. Align to bisulfite-converted genome
bismark --genome /path/to/genome/ -1 sample_R1_val_1.fq.gz -2 sample_R2_val_2.fq.gz

# 4. Remove PCR duplicates (skip for RRBS - fragments are enzyme-defined, not random)
deduplicate_bismark --paired sample_bismark_bt2_pe.bam

# 5. Extract per-cytosine methylation calls
bismark_methylation_extractor --paired-end --CpG --comprehensive \
    --cytosine_report --genome_folder /path/to/genome/ \
    sample_bismark_bt2_pe.deduplicated.bam

# 6. QC report (alignment rate, methylation bias, non-conversion rate)
bismark2report
```

**Goal:** Load Bismark coverage output and filter by read depth before downstream analysis.
**Approach:** the 6-column `.cov` file gives methylated/unmethylated counts per CpG; require a minimum depth to trust the beta value.

```python
import pandas as pd
import matplotlib.pyplot as plt


def load_bismark_cov(path: str, min_coverage: int = 10) -> pd.DataFrame:
    """Load a Bismark .cov(.gz) file and compute beta values, filtered by depth.

    Columns are: chrom, start, end, pct_methylated, count_methylated, count_unmethylated.
    """
    cov = pd.read_csv(
        path, sep="\t", compression="infer",
        names=["chrom", "start", "end", "pct_meth", "count_M", "count_U"],
    )
    cov["coverage"] = cov["count_M"] + cov["count_U"]
    cov["beta"] = cov["count_M"] / cov["coverage"]
    return cov[cov["coverage"] >= min_coverage].copy()


cpg = load_bismark_cov("sample.bismark.cov.gz", min_coverage=10)
plt.hist(cpg["beta"], bins=50, edgecolor="black")
plt.xlabel("Methylation level (beta)")
plt.ylabel("CpG count")
plt.title("CpG methylation distribution")
plt.show()
```

**Goal:** Call differentially methylated positions (DMPs) between two groups.
**Approach:** `methylKit` merges per-sample coverage into one object, filters/normalizes, then runs a logistic regression test per CpG.

```r
library(methylKit)

# Load samples (paths to Bismark .cov files)
methyl_list <- methRead(
  list("ctrl.cov", "treated.cov"),
  sample.id = list("ctrl", "treated"),
  assembly = "hg38", treatment = c(0, 1),
  context = "CpG", mincov = 10
)

# Filter low/extreme coverage, then normalize between samples
filtered  <- filterByCoverage(methyl_list, lo.count = 10, hi.perc = 99.9)
normalized <- normalizeCoverage(filtered)

# Merge samples at common CpGs and test each site (logistic regression)
united     <- unite(normalized, destrand = FALSE)
dm_results <- calculateDiffMeth(united)

# Keep DMPs with >=25 percentage-point change and q < 0.05
dmps <- getMethylDiff(dm_results, difference = 25, qvalue = 0.05)
```

**Goal:** Estimate DNAm (epigenetic) age from a beta-value matrix using the Horvath clock.
**Approach:** dot-product the sample's beta values at 353 clock CpGs with published coefficients, then apply the clock's inverse-logit age transform.

```python
import numpy as np
import pandas as pd


def anti_trafo(x: float, adult_age: int = 20) -> float:
    """Invert Horvath's age transform to recover years from the linear predictor."""
    if x < 0:
        return (1 + adult_age) * np.exp(x) - 1
    return (1 + adult_age) * x + adult_age


def predict_horvath_age(beta_matrix: pd.DataFrame, clock_coef: pd.DataFrame) -> pd.Series:
    """Predict DNAm age from a samples x CpGs beta matrix.

    clock_coef must have columns 'CpGmarker' and 'CoefficientTraining'.
    """
    cpgs = clock_coef["CpGmarker"]
    beta_clock = beta_matrix[cpgs]
    coef = clock_coef.set_index("CpGmarker")["CoefficientTraining"]
    linear_predictor = beta_clock.dot(coef)
    return linear_predictor.apply(anti_trafo)


clock = pd.read_csv("horvath_clock_cpgs.csv")  # CpGmarker, CoefficientTraining
predicted_age = predict_horvath_age(beta_matrix, clock)
```

## Pitfalls
- **Non-conversion rate > 1%**: poor bisulfite conversion; check the lambda/unmethylated spike-in control
- **Coverage asymmetry**: RRBS enriches CpG-dense regions (islands); WGBS is genome-wide but shallower per site
- **Strand merging**: use `destrand=TRUE` in `methRead`/`unite` to merge symmetric CpG strands for extra power
- **Deduplication on RRBS**: don't run `deduplicate_bismark` on RRBS — fragment ends are enzyme-defined, not random, so "duplicates" are real
- **Array vs WGBS**: Illumina EPIC covers ~850K CpGs; WGBS covers ~28M; only partial overlap, so clocks trained on arrays need array-covered CpGs
- **Clock tissue bias**: Horvath is pan-tissue; Hannum and PhenoAge are blood-optimized; GrimAge is mortality-tuned, not raw age

## See Also
- `bio-applied-dmr-analysis` — deeper DMR-calling workflows and statistics
- `bio-applied-epigenetic-clocks` — extended epigenetic clock methods (Hannum, PhenoAge, GrimAge)
- `chipseq-epigenomics` — related chromatin/epigenomic peak analysis
- `algo-sequence-alignment` — background on the alignment algorithms Bismark builds on