ENCODE WGBS Pipeline: FASTQ to Methylation Calls

When to Use

• User wants to run a WGBS/bisulfite sequencing pipeline from FASTQ to methylation calls
• User asks about "WGBS pipeline", "bisulfite sequencing", "methylation calling", "Bismark", or "bedMethyl"
• User needs to process whole-genome bisulfite sequencing data following ENCODE standards
• Example queries: "process my WGBS FASTQs", "call methylation levels from bisulfite-seq", "run Bismark on my WGBS data"

Execute the ENCODE DNA methylation pipeline for Whole Genome Bisulfite Sequencing data, producing per-CpG methylation levels in bedMethyl format.

Pipeline Overview

FASTQ -> Trim adapters -> Bismark align -> Deduplicate -> MethylDackel extract -> bedMethyl
  |           |               |                |                |                    |
  QC      Trim Galore    Bismark/bwa-meth   Picard         Per-CpG calls      Final output

ENCODE Repository

• GitHub: ENCODE-DCC/dna-me-pipeline
• Container: encodedcc/dna-me-pipeline
• WDL: Available for Cromwell execution
• This skill: Nextflow DSL2 reimplementation for portability

Core Tools and Versions

Key Literature

1. Krueger & Andrews 2011 - "Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications" (Bioinformatics, ~4,000 citations) DOI: 10.1093/bioinformatics/btr167

2. Lister et al. 2009 - "Human DNA methylomes at base resolution show widespread epigenomic differences" (Nature, ~5,000 citations) DOI: 10.1038/nature08514

3. Schultz et al. 2015 - "Human body epigenome maps reveal noncanonical DNA methylation variation" (Nature, ~1,500 citations) DOI: 10.1038/nature14248

4. Pedersen et al. 2014 - "Fast and accurate alignment of long bisulfite-seq reads" arXiv:1401.1129 (bwa-meth)

5. Amemiya et al. 2019 - "The ENCODE Blacklist" (Scientific Reports, ~1,372 citations) DOI: 10.1038/s41598-019-45839-z

Execution

Quick Start (Local)

nextflow run main.nf \
    -profile local \
    --reads '/data/fastq/*_R{1,2}.fastq.gz' \
    --genome_dir '/ref/bismark_index' \
    --outdir results/ \
    -resume

SLURM HPC

nextflow run main.nf \
    -profile slurm \
    --reads '/data/fastq/*_R{1,2}.fastq.gz' \
    --genome_dir '/ref/bismark_index' \
    --outdir results/ \
    -resume

Cloud (GCP / AWS)

nextflow run main.nf \
    -profile gcp \
    --reads 'gs://bucket/fastq/*_R{1,2}.fastq.gz' \
    --genome_dir 'gs://bucket/ref/bismark_index' \
    --outdir 'gs://bucket/results/' \
    -resume

Resource Requirements

Pipeline Parameters

Output Files

results/
  fastqc/                    # Raw read quality
  trim_galore/               # Trimmed reads + reports
  bismark/
    alignments/              # Sorted, deduplicated BAMs
    dedup_reports/           # Duplication metrics
    methylation/             # bedMethyl files (primary output)
      {sample}.CpG.bedMethyl.gz
      {sample}.CHG.bedMethyl.gz   # Non-CpG contexts
      {sample}.CHH.bedMethyl.gz
    conversion_rate/         # Lambda/pUC19 conversion QC
  coverage/
    {sample}.coverage_stats.txt
  multiqc/
    multiqc_report.html

bedMethyl Format

The primary output is per-CpG methylation in bedMethyl format:

chr1  10468  10470  .  1000  +  10468  10470  0,0,0  12  83.3

Columns: chr, start, end, name, score, strand, thickStart, thickEnd, color, coverage, methylation_percentage

QC Thresholds (ENCODE Standards)

Critical Pitfalls

RRBS vs WGBS

RRBS (Reduced Representation) uses MspI digestion and covers ~10% of CpGs. WGBS covers the full genome. These are DIFFERENT protocols:

• RRBS: Skip deduplication (--skip_dedup true), different trimming
• WGBS: Full dedup required, standard Trim Galore settings
• Never mix RRBS and WGBS data in the same analysis

Strand-Specific vs Merged CpG

Bismark reports methylation per strand by default. For most analyses, merge complementary CpG strands:

• Forward C at position N and reverse G at position N+1 are the same CpG
• MethylDackel --mergeContext handles this automatically
• Always specify --mergeContext unless you need strand-specific data

Incomplete Bisulfite Conversion

Conversion artifacts produce false methylation calls:

• Always include lambda phage or pUC19 spike-in DNA
• Unmethylated spike-in should show ≥98% conversion
• If conversion <98%, the library has systematic artifacts -- do NOT proceed

M-bias Plots

MethylDackel generates M-bias plots showing methylation level by read position:

• End-repair artifacts cause elevated methylation at read ends
• Use --OT and --OB flags to trim affected positions
• ENCODE typically trims 10 bp from 5' of read 2

Low Coverage Regions

Regions with <5x coverage have unreliable methylation estimates:

• Filter bedMethyl to --min_coverage 5 (default)
• For differential methylation analysis, consider --min_coverage 10
• Report the fraction of CpGs meeting coverage threshold

Provenance Integration

After pipeline completion, log all outputs:

# Log derived bedMethyl files
encode_log_derived_file(
    file_path="/results/bismark/methylation/sample1.CpG.bedMethyl.gz",
    source_accessions=["ENCSR...", "ENCFF..."],
    description="CpG methylation calls from ENCODE WGBS pipeline",
    file_type="bedMethyl",
    tool_used="Bismark 0.24.2 + MethylDackel 0.6.1",
    parameters="bismark --genome /ref -1 R1.fq.gz -2 R2.fq.gz; MethylDackel extract --mergeContext --minDepth 5"
)

Reference Files

Detailed step-by-step documentation is provided in the references/ directory:

1. 01-qc-trimming.md -- Bisulfite-specific adapter trimming with Trim Galore
2. 02-bismark-alignment.md -- Bismark alignment and bwa-meth alternative
3. 03-dedup-filtering.md -- Deduplication and BAM filtering
4. 04-methylation-calling.md -- MethylDackel extraction and bedMethyl generation
5. 05-qc-metrics.md -- Conversion rate QC, coverage stats, M-bias

Walkthrough: Processing ENCODE WGBS