bio-applied-wgbs-bismark

name	bio-applied-wgbs-bismark
description	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.
tool_type	bash
primary_tool	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.

# 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, )
     np.log2((beta + eps) / ( - beta + eps))


 ():
    
    lo = stats.binom.ppf(alpha / , coverage, beta_true) / coverage
    hi = stats.binom.ppf( - alpha / , coverage, beta_true) / coverage
     hi - lo


 __name__ == :
    
    rng = np.random.default_rng()
    n = 
    cov = rng.negative_binomial(,  / , n).clip(=)
    m = rng.binomial(cov, )
    demo_df = pd.DataFrame({
        : , : np.arange(n), : np.arange(n) + ,
        : , : m, : cov - m,
    })
    demo_df[] = cov
    demo_df[] = demo_df[] / demo_df[]

    filtered = filter_by_coverage(demo_df, min_cov=)
     filtered[].() >= 
     filtered[].between(, ).()

    mvals = beta_to_mvalue(filtered[])
     np.isfinite(mvals).()

     binomial_ci_width() > binomial_ci_width()  
    ()
    ()

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.