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 -> sort + index -> MethylDackel -> bedMethyl
  |           |                 |               |                |             |             |
 FastQC   Trim Galore    Bismark (Bowtie2)  deduplicate_   samtools      mbias +       coverage
                                             bismark                     extract        stats

Paired-end only: TRIM_GALORE and BISMARK_ALIGN both take two mates, and --reads must match pairs.

ENCODE Repository

• GitHub: ENCODE-DCC/dna-me-pipeline
• Container: built from scripts/Dockerfile in this skill (docker build -t encode-toolkit/pipeline-wgbs:1.0.0 scripts/); override with --container
• WDL: Available for Cromwell execution
• This skill: Nextflow DSL2 reimplementation for portability

Core Tools and Versions

These are the versions in scripts/Dockerfile, which is what the workflow runs.

The conda environment in bioinformatics-installer (environments/wgbs-env.yml) is a separate manual route pinned to the same versions of these tools; it also carries bedtools, which is not in the image.

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

--reads and --genome_dir are both required; the workflow stops before the first task if either is missing. --genome_dir is a Bismark genome folder — the output of bismark_genome_preparation, which must still contain the genome .fa, because MethylDackel reads it from there.

Quick Start (local, Docker)

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

SLURM HPC

nextflow run scripts/main.nf -profile slurm \
    --container /path/to/pipeline-wgbs.sif \
    --slurm_queue normal \
    --reads '/data/fastq/*_R{1,2}.fastq.gz' \
    --genome_dir /ref/bismark_index \
    --outdir results/ \
    -resume

Cloud

# Google Cloud Batch
nextflow run scripts/main.nf -profile gcp \
    --container us-docker.pkg.dev/<project>/<repo>/pipeline-wgbs:1.0.0 \
    --gcp_project <project> \
    --gcp_workdir gs://<bucket>/work \
    --reads 'gs://<bucket>/fastq/*_R{1,2}.fastq.gz' \
    --genome_dir gs://<bucket>/ref/bismark_index \
    --outdir gs://<bucket>/results

# AWS Batch
nextflow run scripts/main.nf -profile aws \
    --container <account>.dkr.ecr.<region>.amazonaws.com/pipeline-wgbs:1.0.0 \
    --aws_queue <job-queue> \
    --aws_workdir s3://<bucket>/work \
    --reads 's3://<bucket>/fastq/*_R{1,2}.fastq.gz' \
    --genome_dir s3://<bucket>/ref/bismark_index \
    --outdir s3://<bucket>/results

--outdir only sets where results are published; Google Batch and AWS Batch stage every task through the work directory, and the workflow stops with an error if it or the project/queue is missing.

Resource Requirements

Every process asks for memory { N.GB * task.attempt }, so a task killed for running out of memory is retried with more: the second attempt gets twice the figure in the table, the third three times it, bounded by --max_memory. nextflow.config scales the time the same way for BISMARK_ALIGN, DEDUPLICATE and METHYLDACKEL_EXTRACT, bounded by --max_time; the other processes declare no time limit. A task is retried only for exit codes 130-145 and 104 (killed for exceeding a limit); any other failure stops the run.

Pipeline Parameters

Infrastructure parameters (nextflow.config)

Notes on what the workflow does and does not do:

• Aligner: Bismark with Bowtie2. bwa-meth is described in references/02-bismark-alignment.md

as a manual alternative; it is not an option of this workflow and is not in the image.

• Deduplication: deduplicate_bismark. Picard is not installed and is not used.
• BAM filtering: there is no MAPQ or flag filtering step. The BAM handed to

MethylDackel is the Bismark output, deduplicated and sorted. references/03-dedup-filtering.md gives the manual filtering commands if you want them.

• Overlapping mates: MethylDackel never counts both mates of an overlapping read pair, so

there is no switch for it. --merge_context is a separate choice about per-CpG versus per-cytosine output.

• Read-position trimming at extraction: none. MethylDackel extract is called without

--nOT/--nOB, because --clip_R2 10 at the trimming step already removes the end-repair bias. See the M-bias pitfall below.

• Bisulfite conversion rate: not computed by the workflow. Estimate it from the CHH

methylation percentage in the Bismark alignment report, or align a lambda/pUC19 spike-in separately when the sample carries real non-CpG methylation (ESCs, neurons). references/05-qc-metrics.md has both commands.

Output Files

results/
  fastqc/                      # FastQC on the raw reads
  trim_galore/                 # trimmed reads, trimming reports, and FastQC on the trimmed reads
  bismark/
    alignments/
      {sample}.sorted.bam      # final BAM: deduplicated unless --skip_dedup, sorted
      {sample}.sorted.bam.bai
      *_PE_report.txt          # Bismark alignment report (mapping rate, context methylation)
    dedup_reports/
      *.deduplication_report.txt
    mbias/
      {sample}_mbias_*.svg     # M-bias plots, one per strand
      {sample}_mbias_report.txt  # MethylDackel's suggested inclusion bounds
    methylation/
      {sample}_CpG.bedGraph    # MethylDackel output, every covered cytosine/CpG
      {sample}_CHG.bedGraph
      {sample}_CHH.bedGraph
      {sample}.CpG.bedMethyl.gz (+ .tbi)   # primary output, filtered by --min_coverage
      {sample}.CHG.bedMethyl.gz (+ .tbi)   # non-CpG contexts
      {sample}.CHH.bedMethyl.gz (+ .tbi)
  coverage/
    {sample}.coverage_stats.txt
  multiqc/
    multiqc_report.html
  pipeline_info/
    timeline.html, report.html, trace.txt  # Nextflow execution reports

Note the naming: the MethylDackel bedGraphs use an underscore before the context ({sample}_CpG.bedGraph), the bedMethyl files a dot ({sample}.CpG.bedMethyl.gz).

bedMethyl Format

The primary output is per-CpG methylation in the ENCODE bedMethyl layout (https://www.encodeproject.org/data-standards/wgbs/):

chr1  10468  10470  .  12  .  10468  10470  0,0,0  12  83

Columns: chr, start, end, name (.), score, strand, thickStart, thickEnd, colour, coverage, percent methylated. The score in column 5 is the read count capped at 1000, the strand in column 6 is . because MethylDackel's bedGraph carries no strand, and column 11 is an integer percentage. Sites below --min_coverage are dropped.

Coverage Statistics

coverage/{sample}.coverage_stats.txt is computed from the unfiltered {sample}_CpG.bedGraph, so it describes every covered site, not just the reported ones. It has five lines — this is a run at the default --min_coverage 5:

Covered CpGs (>=1x): 27184023
Mean coverage of covered CpGs: 12.4
Covered CpGs >=5x: 22903511 (84.3%)
Covered CpGs >=10x: 16992841 (62.5%)
Covered CpGs >=5x (--min_coverage, kept in bedMethyl): 22903511 (84.3%)

The 5x and 10x lines are fixed thresholds and do not follow --min_coverage. The last line does, and it is the one that describes what actually reached the bedMethyl files: at the default it repeats the >=5x line, with --min_coverage 10 it would repeat the >=10x line, and at any other value it stands on its own. All the percentages are of covered CpGs, not of all CpGs in the genome — the workflow never counts genomic CpGs that got zero reads. With --merge_context false the records are per cytosine rather than per CpG.

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: this workflow is not set up for it. Deduplication can be skipped

(--skip_dedup true), but the trimming is WGBS-specific and hardcoded (--clip_R2 10 --three_prime_clip_R1 1), there is no --rrbs switch, and the workflow re-trims whatever FASTQs it is given. Run RRBS trimming and alignment by hand — references/01-qc-trimming.md has the command

• 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
• --merge_context true is the workflow default; set it to false only when 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
• The workflow does not measure this. Align the spike-in separately, or read the CHH

methylation percentage out of the Bismark report, before trusting the calls

M-bias Plots

MethylDackel mbias runs on every sample and writes plots plus bismark/mbias/{sample}_mbias_report.txt, which captures MethylDackel's suggested inclusion bounds. Extraction itself excludes no read positions:

• End-repair artifacts cause elevated methylation at read ends. The 5' end-repair bias of

read 2 is already removed at the trimming step by --clip_R2 10

• If the plots still show bias, re-run MethylDackel extract manually with --OT/--OB

set to the suggested bounds. Note the argument order: --OT A,B,C,D includes positions A-B on read 1 and C-D on read 2, and the matching --nOT a,b,c,d excludes a bases from the start of read 1, b from its end, c from the start of read 2 and d from its end

• --maxDepth is not an option of MethylDackel 0.6.1 and must not be passed

Low Coverage Regions

Regions with <5x coverage have unreliable methylation estimates:

• The bedMethyl files are filtered to --min_coverage (default 5); the bedGraphs are not
• For differential methylation analysis, consider --min_coverage 10
• Report the >=5x and >=10x fractions from coverage/{sample}.coverage_stats.txt. Those

two thresholds are fixed in the workflow and do not follow --min_coverage; the last line of the file is the one that does, so quote it too whenever the run used a value other than 5

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 --CHG --CHH; bedMethyl filtered at --min_coverage 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 the bwa-meth manual alternative
3. 03-dedup-filtering.md -- Deduplication, and manual BAM filtering the workflow skips
4. 04-methylation-calling.md -- MethylDackel extraction and bedMethyl generation
5. 05-qc-metrics.md -- Coverage stats, M-bias, and manual conversion-rate checks

Walkthrough: Processing ENCODE WGBS from FASTQ to Methylation Calls

Goal: Process whole-genome bisulfite sequencing FASTQ files through the ENCODE pipeline to generate per-CpG methylation calls for epigenomic analysis. Context: WGBS requires bisulfite-aware alignment (Bismark) and per-CpG methylation extraction (MethylDackel), with ≥98% bisulfite conversion expected of the library.

Step 1: Find WGBS experiment

encode_get_experiment(accession="ENCSR765JPC")

Expected output:

{
  "accession": "ENCSR765JPC",
  "assay_title": "WGBS",
  "biosample_summary": "liver tissue male adult (54 years)",
  "assembly": ["GRCh38"],
  "bio_replicate_count": 2,
  "tech_replicate_count": 2,
  "status": "released"
}

Step 2: List FASTQ files

encode_list_files(experiment_accession="ENCSR765JPC", file_format="fastq")

Expected output (a JSON array of file records; fields abridged):

[
  {"accession": "ENCFF300BS1", "file_format": "fastq", "output_type": "reads", "file_size_human": "45.0 GB", "biological_replicates": [1], "status": "released"},
  {"accession": "ENCFF301BS2", "file_format": "fastq", "output_type": "reads", "file_size_human": "46.0 GB", "biological_replicates": [1], "status": "released"}
]

Interpretation: WGBS files are very large (~45GB per read file). Ensure adequate storage (>500GB for processing).

Step 3: Download and name the FASTQs so a read-pair glob can find them

encode_download_files(file_accessions=["ENCFF300BS1", "ENCFF301BS2"], download_dir="/data/wgbs/fastq")

ENCODE names every FASTQ after its accession (ENCFF300BS1.fastq.gz), with no _R1/_R2 in the name, so the two files of a pair share no prefix and the --reads glob cannot pair them. Link them into the shape the glob expects. Which mate an accession is comes from the ENCODE file record on encodeproject.org, which carries paired_end (1 or 2) and paired_with; the MCP file tools do not return those two fields:

cd /data/wgbs/fastq
ln -s ENCFF300BS1.fastq.gz liver_rep1_R1.fastq.gz
ln -s ENCFF301BS2.fastq.gz liver_rep1_R2.fastq.gz

Step 4: Run the WGBS pipeline

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

Key pipeline steps:

1. FastQC on the raw reads, then adapter/quality trimming with Trim Galore (which also runs FastQC on the trimmed reads)
2. Bisulfite-aware alignment (Bismark with Bowtie2)
3. Deduplication (deduplicate_bismark), then sort and index
4. M-bias assessment and methylation extraction (MethylDackel), then conversion to bedMethyl
5. CpG coverage statistics
6. MultiQC aggregation

Step 5: Validate output quality

Step 6: Identify differentially methylated regions

Feed per-CpG methylation into -> methylation-aggregation for cross-tissue comparison and HMR/UMR/PMD identification.

Integration with downstream skills

• Per-CpG methylation files feed into -> methylation-aggregation for cross-tissue atlas
• DMRs feed into -> peak-annotation for nearest gene assignment
• Methylation at regulatory elements connects to -> regulatory-elements
• CpG variant methylation integrates with -> variant-annotation
• Pipeline provenance logged by -> data-provenance

Code Examples

1. Find WGBS data for methylation analysis

encode_search_experiments(
  assay_title="WGBS",
  organ="brain"
)

Expected output:

{
  "results": [
    {
      "accession": "ENCSR321BRN",
      "assay_title": "WGBS",
      "biosample_summary": "brain tissue female adult (53 years)",
      "organ": "brain",
      "status": "released"
    }
  ],
  "total": 6,
  "limit": 25,
  "offset": 0,
  "has_more": false,
  "next_offset": null
}

2. Download processed methylation files

encode_search_files(
  assay_title="WGBS",
  organ="brain",
  file_format="bed",
  output_type="methylation state at CpG",
  assembly="GRCh38"
)

Expected output (fields abridged):

{
  "results": [
    {
      "accession": "ENCFF567MET",
      "file_format": "bed",
      "output_type": "methylation state at CpG",
      "assembly": "GRCh38",
      "file_size": 886144860,
      "file_size_human": "845.1 MB",
      "experiment_accession": "ENCSR321BRN"
    }
  ],
  "total": 3,
  "limit": 25,
  "offset": 0,
  "has_more": false,
  "next_offset": null
}

Integration

Related Skills

• pipeline-guide -- Parent skill with compute resource assessment and cloud setup
• methylation-aggregation -- Aggregate methylation data across samples/tissues
• quality-assessment -- Evaluate pipeline output quality metrics
• data-provenance -- Track all pipeline inputs, outputs, and parameters
• download-encode -- Download ENCODE WGBS FASTQ files for pipeline input
• publication-trust -- Verify literature claims backing analytical decisions

Presenting Results

When reporting WGBS pipeline results:

• Mapping and duplication: Report the mapping rate from bismark/alignments/*_PE_report.txt (>70% expected) and the duplication rate from bismark/dedup_reports/ (<30% expected)
• CpG coverage depth: Report the mean coverage of covered CpGs, the percentage reaching >=5x and >=10x, and the last line's --min_coverage count from coverage/{sample}.coverage_stats.txt, and state that these are percentages of covered CpGs rather than of all genomic CpGs
• Global methylation level: Report the genome-wide average CpG methylation percentage and note any non-CpG (CHG/CHH) methylation if relevant to the tissue. Both come from the Bismark report or from the bedGraphs
• bedMethyl output paths: Provide paths to the CpG bedMethyl file and the CHG/CHH context files, and state the --min_coverage value they were filtered at
• M-bias assessment: No read positions are excluded at extraction. Report MethylDackel's suggested bounds from bismark/mbias/{sample}_mbias_report.txt, and say whether a manual re-extraction with --OT/--OB was done
• Bisulfite conversion rate: This workflow does not compute it. If a lambda/pUC19 spike-in was aligned separately, or the CHH percentage was used as a proxy, report the number and say which method produced it (≥98% pass, 95-98% warning, <95% fail)
• Output summary: Report covered CpG sites, sites passing the coverage filter, and the MultiQC report at multiqc/multiqc_report.html
• Next steps: Suggest methylation-aggregation for cross-sample averaging and HMR/UMR/PMD identification

For the request: "$ARGUMENTS"