dna-methylation

name	dna-methylation
description	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).
tool_type	bash
primary_tool	bismark

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.

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

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.

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  
dm_results  calculateDiffMethunited


dmps  getMethylDiffdm_results difference   qvalue  

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.

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