Bisulfite Sequencing Explained: From Raw Reads to Methylation Calls with Bismark

DNA methylation is a chemical mark, not a change to the sequence itself — a methyl group parked on a cytosine that a plain FASTQ file has no way to record. So how do you read an invisible modification with a sequencer that only reports A, C, G, and T? The answer is a clever chemical trick called bisulfite conversion, which rewrites the genome so that methylation status becomes a difference in the letters themselves. This post walks the full path: the conversion chemistry, why it breaks normal aligners, how Bismark solves the resulting three-letter genome problem, and how you go from raw reads to per-CpG methylation percentages and differentially methylated regions.

What Methylation Is and Why It Matters

In mammals, methylation almost always means 5-methylcytosine (5mC) at cytosines sitting in a CpG dinucleotide — a C immediately followed by a G. These marks silence transposons, imprint genes, and shut down promoters when CpG islands become hypermethylated, a hallmark of many cancers. Because the mark controls whether genes fire without touching the underlying DNA sequence, it is one of the central layers of epigenetic regulation.

The catch is measurement. A methyl group does not change base pairing, so an ordinary sequencing run reads a methylated C and an unmethylated C identically. To see the mark, you need a chemistry that treats the two differently — and that is exactly what sodium bisulfite provides.

The Bisulfite Conversion Trick

Treat DNA with sodium bisulfite and something selective happens: unmethylated cytosines are deaminated to uracil, which is then read as thymine (T) after PCR amplification. Methylated cytosines are protected and stay as C. The methyl group physically shields the base from the reaction.

That single asymmetry is the whole game. After conversion:

• A position that reads C in your data was methylated.
• A position that reads T (but is C in the reference) was unmethylated.

So methylation status is now encoded as a C-versus-T difference at every cytosine. If you sequence a region 30 times and 24 reads show C at a given CpG while 6 show T, that site is roughly 80% methylated. Methylation is no longer invisible — it is spelled out in the letters. The cost is that most of the cytosines in the genome (the unmethylated ones) have been erased and replaced with thymine.

The Three-Letter Genome Problem

Here is why you cannot just hand bisulfite reads to BWA or Bowtie2 and expect sensible results. After conversion, a read derived from an unmethylated region is essentially a three-letter sequence — nearly all cytosines have become thymines, so the read is written in an A/G/T alphabet with very few surviving Cs. A normal aligner scoring that read against a four-letter reference sees a mismatch at every converted C. Those pile up into huge penalties, and the read either fails to map or maps to the wrong place.

To make matters worse, the conversion is strand-asymmetric. The two strands of DNA are no longer complementary after treatment. The original top strand loses its Cs (becoming C→T reads); the original bottom strand independently loses its Cs, which on the top-strand coordinate system looks like G→A changes. A single genomic locus can therefore give rise to reads that look C→T-converted or G→A-converted depending on which strand and which PCR product they came from. You must account for all of them, or you throw away half your data and bias the half you keep.

If the mechanics of scoring mismatches and seeding are unfamiliar, the general primer in how read aligners work is worth a detour before continuing.

How Bismark Solves It

Bismark reduces the four-letter problem to a three-letter one deliberately, then does the alignment there. Its strategy has two phases.

First, a one-time genome preparation with bismark_genome_preparation builds fully C→T-converted and G→A-converted copies of the reference, indexed for Bowtie2 (or HISAT2). Second, at alignment time Bismark performs the same in-silico conversion on your reads: it makes a C→T version and a G→A version of every read, aligns each against the matching converted genome index, and then compares the resulting alignments to decide the best unique hit. This is the dual-strand mapping that handles both the original-top and original-bottom strand products at once.

The elegance is that once everything — reads and genome — lives in the same reduced three-letter alphabet, the converted Cs are no longer mismatches; they simply agree. Bismark then does something the reduced alphabet threw away: it walks each aligned read back against the original unconverted reference and, at every CpG, CHG, and CHH position, asks whether the read shows a C (methylated) or a T (unmethylated). That per-base comparison is where the actual methylation calls come from.

A Practical Bismark Workflow

The end-to-end path from FASTQ to methylation calls is short but each step matters. Adapter and quality trimming come first — Trim Galore is the near-universal choice here because it wraps Cutadapt and FastQC-style QC and has a dedicated --rrbs mode for reduced-representation libraries. If you want the reasoning behind trimming decisions in general, the NGS read QC and trimming walkthrough covers it.

# 1. Trim adapters and low-quality tails
trim_galore --paired --cores 4 sample_R1.fastq.gz sample_R2.fastq.gz

# 2. Prepare the bisulfite genome (once per reference)
bismark_genome_preparation --bowtie2 /path/to/genome/

# 3. Align the trimmed reads
bismark --genome /path/to/genome/ -1 sample_R1_val_1.fq.gz -2 sample_R2_val_2.fq.gz

# 4. Remove PCR duplicates
deduplicate_bismark --paired sample_R1_val_1_bismark_bt2_pe.bam

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

A few notes that save real debugging time. Run deduplicate_bismark on whole-genome libraries but never on RRBS data — RRBS fragments legitimately share start coordinates because they are cut at the same MspI sites, so deduplication would delete true signal. The --cytosine_report flag produces the genome-wide per-CpG coverage file (CX_report) that downstream tools consume. Finally, always read the Bismark HTML report and the alignment rate; bisulfite libraries map at a lower rate than standard genomic DNA, but a rate under roughly 50% usually signals a contamination or conversion problem.

Calling Methylation Percentages and DMRs

The extractor gives you counts, not conclusions. For each cytosine you get the number of methylated (C) and unmethylated (T) reads; the methylation level is simply methylated / (methylated + unmethylated). The biologically interesting question is rarely a single site, though — it is whether a region differs between conditions. That is a differentially methylated region (DMR), and it is analogous to differential expression in RNA-seq. If you have done RNA-seq differential expression, the mental model of replicates, dispersion, and multiple testing carries straight over.

methylKit, an R/Bioconductor package, reads Bismark's coverage files, filters sites by minimum depth (commonly 10x), and tests for differential methylation with logistic regression or Fisher's exact test, correcting for multiple comparisons with a q-value. You typically threshold on both effect size and significance — for example, q < 0.01 and a methylation difference > 25% — to call a site or tiled window as differentially methylated. Because bisulfite coverage is uneven, filtering out sites below a depth floor before testing is not optional; low-coverage sites produce wildly noisy percentages.

Bismark vs bwa-meth: Choosing an Aligner

Bismark is the didactic, well-documented standard, but it is not the only route. bwa-meth wraps BWA-MEM and is popular for large whole-genome bisulfite sequencing (WGBS) projects where speed and downstream compatibility matter, usually paired with MethylDackel for extraction.

Neither is strictly better. Bismark's exhaustive reporting makes it the right tool to learn on and the safer choice for RRBS; bwa-meth plus MethylDackel scales more comfortably to deep WGBS. Both feed cleanly into methylKit or DSS, so the choice rarely locks you out of downstream analysis.

Practical Takeaways

• Bisulfite conversion turns unmethylated C into T while methylated C is protected, so methylation becomes a readable C-vs-T difference at every cytosine.
• The conversion creates a three-letter genome that normal aligners choke on; Bismark solves it by converting both reads and reference to the same reduced alphabet and doing dual-strand mapping.
• Trim first with Trim Galore (use --rrbs for RRBS), and never run deduplicate_bismark on RRBS libraries — shared cut sites are real signal, not duplicates.
• Always check conversion efficiency (via a lambda spike-in, target > 99%) and the M-bias plot before trusting any methylation percentage.
• Use methylKit with a depth filter (e.g. 10x) and joint thresholds like q < 0.01 plus a methylation difference > 25% to call DMRs robustly.
• For large WGBS, consider bwa-meth plus MethylDackel for speed; for RRBS and learning, Bismark's reporting is hard to beat.

Bisulfite sequencing rewards the same discipline as any NGS assay: clean input, honest QC, and an aligner that matches the chemistry. If you are setting up the toolchain, start by getting a reproducible environment in place with conda environments and Bioconda channels, since Bismark, Trim Galore, and methylKit all install cleanly from Bioconda and pin far more reliably than a hand-built stack.