pipeline-wgbs

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

| Tool | Version | Purpose | Citation | |------|---------|---------|----------| | Trim Galore | 0.6.10 | Adapter + quality trimming (bisulfite-aware) | Krueger (Babraham) | | Bismark | 0.24.2 | Bisulfite-aware alignment + methylation | Krueger & Andrews 2011 | | bwa-meth | 0.2.7 | Alternative bisulfite aligner (faster) | Pedersen 2014 | | MethylDackel | 0.6.1 | Methylation extraction from BAM | Ryan (GitHub) | | Picard | 3.1.1 | Duplicate marking | Broad Institute | | samtools | 1.19 | BAM operations | Li et al. 2009 | | FastQC | 0.12.1 | Read quality assessment | Andrews (Babraham) | | MultiQC | 1.21 | Aggregated QC reporting | Ewels et al. 2016 |

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

1. Lister et al. 2009 - "Human DNA methylomes at base resolution show widespread

epigenomic differences" (Nature, ~5,000 citations) DOI: 10.1038/nature08514

1. Schultz et al. 2015 - "Human body epigenome maps reveal noncanonical DNA

methylation variation" (Nature, ~1,500 citations) DOI: 10.1038/nature14248

1. Pedersen et al. 2014 - "Fast and accurate alignment of long bisulfite-seq reads"

arXiv:1401.1129 (bwa-meth)

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

| Step | CPUs | RAM | Time (30x human) | |------|------|-----|-------------------| | Trim Galore | 4 | 4 GB | 1-2 hours | | Bismark align | 8 | 48 GB | 8-16 hours | | Deduplication | 2 | 16 GB | 1-2 hours | | MethylDackel | 4 | 8 GB | 1-2 hours | | Total | 8 | 48 GB | 12-24 hours |

Pipeline Parameters

| Parameter | Default | Description | |-----------|---------|-------------| | --reads | required | Glob pattern to paired FASTQ files | | --genome_dir | required | Path to Bismark genome index directory | | --outdir | ./results | Output directory | | --aligner | bismark | Aligner: bismark or bwameth | | --min_coverage | 5 | Minimum coverage for CpG reporting | | --no_overlap | true | Remove overlapping PE reads (avoid double-counting) | | --lambda_genome | null | Lambda genome index for conversion rate QC | | --skip_dedup | false | Skip deduplication (for RRBS data) |

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)

| Metric | Pass | Warning | Fail | |--------|------|---------|------| | Bisulfite conversion rate | ≥98% | 95-98% | <95% | | CpG coverage (genome-wide) | >10x | 5-10x | <5x | | Mapping rate | >70% | 50-70% | <50% | | Duplication rate | <30% | 30-50% | >50% | | CpG sites covered (>=5x) | >80% | 60-80% | <60% | | Lambda spike-in conversion | ≥98% | 95-98% | <95% |

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