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

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

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.

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

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.

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.