# Bio Applied Wgbs Bismark

> Align WGBS/RRBS bisulfite FASTQ with Bismark, extract per-CpG methylation into beta/M-values. Use for bismark_genome_preparation, deduplicate_bismark, bismark_methylation_extractor, or bismark.cov/CpG_report analysis.

- Skill: `pavel-kravchenko/bio-applied-wgbs-bismark` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-wgbs-bismark`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-wgbs-bismark/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/bio-applied-wgbs-bismark

---


# WGBS/RRBS Processing with Bismark

## When to Use

- Aligning bisulfite-converted FASTQ reads (WGBS or RRBS) to a reference genome.
- Building a Bismark bisulfite genome index and running the full trim → align → dedup → extract pipeline.
- Computing per-CpG beta values / M-values from Bismark coverage output, with coverage-based QC filtering.
- Deciding between WGBS, RRBS, and EPIC array for a methylation study design.
- Producing a genome-wide or regional methylation landscape plot (CpG island vs gene body vs intergenic).

## Version Compatibility

Bismark ≥0.24 (Bowtie2 ≥2.5 or HISAT2 backend), Trim Galore ≥0.6, Python ≥3.10 with pandas ≥2.0 / numpy ≥1.24 / scipy ≥1.10 / matplotlib ≥3.7 for downstream analysis. R ≥4.3 with minfi ≥1.48 (Bioconductor ≥3.18) for array comparison.

## Prerequisites

- `bismark`, `bowtie2` (or `hisat2`), `samtools`, `trim_galore` on PATH.
- A reference genome FASTA (one genome dir per assembly; index is genome-specific and reused across samples).
- Python: `pandas`, `numpy`, `scipy`, `matplotlib`.
- Prior concept: bisulfite chemistry (unmethylated C → U → reads as T; methylated C protected, reads as C) and the resulting 4-strand alignment problem (OT/CTOT/OB/CTOB) that ordinary aligners cannot resolve.
- Related skill: `dna-methylation` (DMR calling, methylKit) and `bio-applied-epigenetic-clocks` (downstream clock models) for what comes after this pipeline.

## Biology Quick Reference

5mC is deposited by DNMT3A/3B (de novo) and DNMT1 (maintenance), removed via TET1/2/3 oxidation to 5hmC. Occurs almost exclusively at CpG dinucleotides in mammals (~70-80% genome-wide methylated).

| Feature | Methylation | Function |
|---|---|---|
| CpG islands (CGIs, ≥200bp, obs/exp CpG >0.6, GC >50%) | Mostly unmethylated | Protect promoters (~70% overlap) |
| CGI shores (±2kb) | Variable, tissue-specific | Major differential methylation site |
| Gene bodies | Moderate | Associated with active transcription |
| Repeats/transposons | Heavy | Silences parasitic elements |

Disease: cancer → CGI promoter hypermethylation silences tumor suppressors (BRCA1, MLH1, CDKN2A) plus global hypomethylation; aging → gradual drift (basis of epigenetic clocks).

| Method | CpG coverage | Cost | Best for |
|---|---|---|---|
| WGBS | ~28M (all) | High (>$500/sample) | Comprehensive/novel DMR discovery |
| RRBS (MspI digest, cuts CCGG) | ~5M (CpG-enriched) | Medium (~6x cheaper) | Cost-effective, most CGIs covered |
| EPIC array | 850K (fixed probes) | Low (~$200) | Large cohorts, no CGI-shore/enhancer coverage |

## Core Pipeline

**Goal:** go from paired-end bisulfite FASTQ to a per-CpG methylation call table.
**Approach:** one-time genome prep, then per-sample trim → align → dedup → extract.

```bash
# Step 1: build bisulfite genome index (once per genome; makes CT_conversion/ and GA_conversion/)
bismark_genome_preparation /path/to/genome/hg38/

# Step 2: adapter/quality trimming (bisulfite libraries need aggressive trimming)
trim_galore --paired --fastqc sample_R1.fastq.gz sample_R2.fastq.gz

# Step 3: bisulfite-aware alignment (expect 60-80% mapping rate for WGBS;
# lower rates indicate conversion failure or contamination)
bismark --genome /path/to/genome/hg38/ \
    -1 sample_R1_val_1.fq.gz -2 sample_R2_val_2.fq.gz \
    --output_dir bismark_output/

# Step 4: remove PCR duplicates (essential for WGBS; less critical for RRBS,
# where genuine read overlap at MspI cut sites is expected)
deduplicate_bismark -p bismark_output/sample_R1_val_1_bismark_bt2_pe.bam

# Step 5: extract per-cytosine methylation calls in all three contexts
bismark_methylation_extractor \
    --paired-end --CpG --CHG --CHH --comprehensive \
    --cytosine_report --genome_folder /path/to/genome/hg38/ \
    bismark_output/sample_R1_val_1_bismark_bt2_pe.deduplicated.bam

# QC summary (alignment rate, non-conversion rate, strand balance) as HTML
bismark2summary bismark_output/*.bam
```

Output formats: `*.bismark.cov.gz` (chrom, start, end, %methylated, count_M, count_U — only covered sites) and `*.CpG_report.txt.gz` (every CpG in the genome, including zero-coverage).

## Beta Values, M-Values, and Coverage Filtering

**Goal:** turn a Bismark coverage table into filtered beta/M-values ready for downstream stats.
**Approach:** compute β = M/(M+U) per CpG, drop low/very-high coverage sites, logit-transform for testing.

```python
import numpy as np
import pandas as pd
from scipy import stats


def load_bismark_cov(path):
    """Parse a *.bismark.cov(.gz) file into a tidy DataFrame with a beta column.

    Bismark .cov format (tab-separated, no header):
    chrom, start, end, pct_methylated, count_methylated, count_unmethylated
    """
    cols = ["chrom", "start", "end", "pct_meth", "count_M", "count_U"]
    df = pd.read_csv(path, sep="\t", header=None, names=cols)
    df["coverage"] = df["count_M"] + df["count_U"]
    df["beta"] = df["count_M"] / df["coverage"]
    return df


def filter_by_coverage(df, min_cov=10, max_pct=99):
    """Drop CpGs below min_cov reads or above the max_pct coverage percentile.

    Low coverage -> unreliable beta (binomial sampling noise).
    Very high coverage -> possible PCR duplicates that survived dedup.
    """
    max_cov = np.percentile(df["coverage"], max_pct)
    mask = (df["coverage"] >= min_cov) & (df["coverage"] <= max_cov)
    return df.loc[mask].copy()


def beta_to_mvalue(beta, eps=0.01):
    """Logit-transform beta values to M-values (log2 scale, more homoscedastic).

    M = log2((beta + eps) / (1 - beta + eps)); preferred for linear-model testing.
    """
    beta = np.clip(beta, 0, 1)
    return np.log2((beta + eps) / (1 - beta + eps))


def binomial_ci_width(coverage, beta_true=0.5, alpha=0.05):
    """95% binomial CI width for a beta estimate at a given coverage (worst case beta=0.5).

    Demonstrates why low-coverage CpGs are unreliable: at 5x coverage the CI
    spans most of [0,1]; by 30x it narrows substantially.
    """
    lo = stats.binom.ppf(alpha / 2, coverage, beta_true) / coverage
    hi = stats.binom.ppf(1 - alpha / 2, coverage, beta_true) / coverage
    return hi - lo


if __name__ == "__main__":
    # demo() self-check with synthetic data standing in for a real .cov file
    rng = np.random.default_rng(42)
    n = 5000
    cov = rng.negative_binomial(5, 5 / 25, n).clip(min=1)
    m = rng.binomial(cov, 0.4)
    demo_df = pd.DataFrame({
        "chrom": "chr1", "start": np.arange(n), "end": np.arange(n) + 1,
        "pct_meth": 0, "count_M": m, "count_U": cov - m,
    })
    demo_df["coverage"] = cov
    demo_df["beta"] = demo_df["count_M"] / demo_df["coverage"]

    filtered = filter_by_coverage(demo_df, min_cov=10)
    assert filtered["coverage"].min() >= 10
    assert filtered["beta"].between(0, 1).all()

    mvals = beta_to_mvalue(filtered["beta"])
    assert np.isfinite(mvals).all()

    assert binomial_ci_width(5) > binomial_ci_width(30)  # more reads -> tighter CI
    print(f"Retained {len(filtered)}/{len(demo_df)} CpGs after coverage filtering")
    print("All checks passed.")
```

## EPIC Array Alternative (R)

**Goal:** compare against Illumina EPIC 850K array data when WGBS is not feasible.
**Approach:** load IDATs with minfi, normalize, extract a beta matrix.

```r
library(minfi)

# Load raw intensity data (one IDAT pair per sample)
RGSet <- read.metharray.exp("idat_directory/")

# Noob background/dye-bias correction (BMIQ can follow for type I/II probe bias)
MSet <- preprocessNoob(RGSet)

# beta: CpGs x samples matrix, values in [0, 1]
beta <- getBeta(MSet)
```

Arrays cover only ~55% of CpG islands and miss CGI shores, enhancers, and non-CpG methylation entirely — use WGBS for novel DMR discovery or cancer epigenome studies.

## Pitfalls

- **Coordinate systems**: `.bismark.cov` is 0-based half-open like BED; mixing with 1-based VCF/GFF coordinates causes off-by-one errors when joining with annotations.
- **Low-coverage beta is noise, not signal**: a CpG with 3 reads showing 2 methylated could plausibly be 33-100% methylated — always apply a ≥10x (ideally ≥20x for strict work) coverage filter before comparing samples.
- **Skipping deduplication on WGBS**: bisulfite conversion reduces sequence complexity, so PCR duplicates map identically far more often than in standard DNA-seq; `deduplicate_bismark` before `bismark_methylation_extractor`, not after.
- **Non-conversion contamination**: check CHH/CHG "methylation" as a proxy for incomplete bisulfite conversion (target >99.5% conversion); high CHH signal genome-wide usually means a QC failure, not real biology.
- **Batch effects and multiple testing**: check for batch confounding before interpreting differential methylation, and apply FDR (Benjamini-Hochberg) when testing thousands of CpGs/regions.

## See Also

- `dna-methylation` — differentially methylated region (DMR) calling downstream of this pipeline.
- `bio-applied-epigenetic-clocks` — age-prediction models built on beta-value matrices.