# QC and Trimming for WGBS Data

Bisulfite-converted reads require specialized trimming to handle adapter
contamination and conversion-induced sequence biases.

## Pre-Trimming QC with FastQC

```bash
fastqc --threads 4 --outdir fastqc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz
```

Key checks before trimming:
- Per-base sequence quality (expect Phred >28 across most positions)
- Adapter content (Illumina universal adapter is common)
- Per-base sequence content (bisulfite conversion causes C->T skew -- this is EXPECTED)
- Overrepresented sequences

**Important**: FastQC will flag per-base sequence content as FAIL for bisulfite data.
This is normal and expected because unmethylated C is converted to T.

## Adapter Trimming with Trim Galore

Trim Galore auto-detects Illumina adapters and applies bisulfite-aware trimming. This is
the command the workflow runs, with `--cores` set from the process CPU count:

```bash
trim_galore \
    --paired \
    --quality 20 \
    --phred33 \
    --length 36 \
    --cores 4 \
    --clip_R2 10 \
    --three_prime_clip_R1 1 \
    --fastqc \
    sample_R1.fastq.gz \
    sample_R2.fastq.gz
```

The trimmed reads, the trimming reports and the FastQC reports that `--fastqc` produces
for the trimmed reads are all published to `trim_galore/`, and the trimming and trimmed-read
FastQC reports are fed to MultiQC. FastQC on the *raw* reads is a separate process and
lands in `fastqc/`.

### Parameter Rationale

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--quality 20` | Phred 20 | Standard quality cutoff |
| `--length 36` | 36 bp | Minimum length after trimming |
| `--clip_R2 10` | 10 bp | Remove end-repair artifacts from read 2 (5' end) |
| `--three_prime_clip_R1 1` | 1 bp | Remove filled-in C from end repair |
| `--cores 4` | 4 | Parallel processing (Trim Galore uses 3x cores internally) |

### RRBS-Specific Trimming (not reachable from this workflow)

The trimming above is hardcoded and WGBS-specific: there is no `--rrbs` switch and no
parameter that changes the clip settings, and the workflow re-trims whatever it is given.
For RRBS data, run Trim Galore by hand with the `--rrbs` flag to handle MspI-digested
fragments, and then run Bismark and MethylDackel on the result yourself rather than through
this workflow:

```bash
trim_galore \
    --paired \
    --rrbs \
    --quality 20 \
    --length 36 \
    --fastqc \
    --cores 4 \
    --output_dir trim_galore/ \
    sample_R1.fastq.gz \
    sample_R2.fastq.gz
```

The `--rrbs` flag removes 2 bp from the 3' end of reads that were filled in
during end-repair of MspI-digested fragments (CCGG sites).

## Post-Trimming QC

Trim Galore produces a trimming report with:
- Total reads processed
- Reads with adapter detected
- Reads too short after trimming
- Quality-trimmed bases

Verify in the report:
- Adapter detection rate should be >0% (if 0%, adapters may not be standard)
- Reads passing filter should be >90%
- Read length distribution should be reasonable (median >50 bp for 150 bp reads)

# Bismark Alignment for WGBS

Bisulfite-converted reads require a specialized aligner that accounts for
C-to-T conversion. Bismark handles this by aligning to both C-to-T and G-to-A
converted genomes simultaneously.

## Genome Preparation (one-time prep, outside the workflow)

Build the bisulfite-converted genome index. The directory you prepare is what you pass as
`--genome_dir`, and it must still contain the genome `.fa`, because MethylDackel reads the
FASTA from the same directory later.

```bash
bismark_genome_preparation \
    --bowtie2 \
    --parallel 4 \
    --verbose \
    /ref/genome/
```

This creates two converted genomes:
- C-to-T converted (for original top strand and complementary to bottom)
- G-to-A converted (for original bottom strand and complementary to top)

Requires approximately 12 GB disk space for human genome.

## Bismark Alignment

This is what the workflow runs, in the task directory:

```bash
mkdir -p tmp
bismark \
    --genome /ref/genome/ \
    --bowtie2 \
    --parallel 4 \
    --score_min L,0,-0.2 \
    --no_mixed \
    --no_discordant \
    --maxins 1000 \
    --temp_dir $PWD/tmp \
    -1 sample_R1_val_1.fq.gz \
    -2 sample_R2_val_2.fq.gz
```

Of the two outputs, only the `*_PE_report.txt` is published (to
`bismark/alignments/`). The raw unsorted BAM stays in the work directory; what reaches
`bismark/alignments/<sample>.sorted.bam` is the deduplicated, coordinate-sorted BAM from
the later step.

### Key Parameters

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--parallel 4` | 4 instances | Each uses 2 Bowtie2 threads = 8 total threads |
| `--score_min L,0,-0.2` | Linear penalty | ENCODE default; tolerates bisulfite mismatches |
| `--no_mixed` | Discard | Both mates must align |
| `--no_discordant` | Discard | Mates must be properly paired |
| `--maxins 1000` | 1000 bp | Maximum insert size for paired reads |

### Memory Requirements

Bismark parallel mode runs several Bowtie2 instances, each holding its own copy of a
converted genome index. `nextflow.config` allocates 48 GB to `BISMARK_ALIGN` with
`--parallel 4` on a human genome, and multiplies it by the attempt number on a retry,
capped by `--max_memory` (64 GB by default, so a retry gets 64 GB). Budget for the 48 GB
figure rather than a per-instance estimate.

## Alternative: bwa-meth

bwa-meth is not part of the bundled workflow or container image. The commands below are for
running it by hand when Bismark is too slow or memory-hungry for your genome.

bwa-meth is faster than Bismark for large genomes and uses less memory:

```bash
# Index genome (one-time)
bwameth.py index /ref/genome/genome.fa

# Align
bwameth.py \
    --threads 8 \
    --reference /ref/genome/genome.fa \
    sample_R1_trimmed.fq.gz \
    sample_R2_trimmed.fq.gz \
    | samtools sort -@ 4 -o sample_bwameth.bam

samtools index sample_bwameth.bam
```

### Bismark vs bwa-meth Comparison

| Feature | Bismark | bwa-meth |
|---------|---------|----------|
| Speed | Slower (~2x) | Faster |
| RAM | 48 GB allocated (parallel) | ~16 GB |
| Accuracy | Gold standard | Comparable |
| This workflow | The only aligner it runs | Not available; manual only |
| Methylation calling | Built-in | Requires MethylDackel |

Recommendation: Use Bismark for ENCODE compatibility. Use bwa-meth outside this workflow
when processing many samples and speed is critical.

## Lambda/pUC19 Spike-in Alignment (not run by this workflow)

If spike-in DNA was used, align to the spike-in genome by hand to measure the conversion
rate. The workflow does not do this and reports no conversion rate:

```bash
bismark \
    --genome /ref/lambda/ \
    --bowtie2 \
    -1 sample_R1_trimmed.fq.gz \
    -2 sample_R2_trimmed.fq.gz \
    --output_dir lambda_out/ \
    --unmapped
```

The lambda genome is fully unmethylated, so any detected methylation represents
incomplete bisulfite conversion. Expect conversion rate ≥98%.

## Alignment QC Checks

After alignment, verify in the Bismark report (`bismark/alignments/*_PE_report.txt`,
also parsed into the MultiQC report):
- **Mapping efficiency**: >70% for WGBS (lower than standard WGS due to conversion)
- **Unique alignments**: Should dominate over multimappers
- **C methylated in CpG context**: Typically 70-85% for mammalian somatic tissue
- **C methylated in CHG context**: Should be <1% in somatic tissue (>5% in plants/ESCs)
- **C methylated in CHH context**: Should be <1% in somatic tissue

# Deduplication and Filtering for WGBS

PCR duplicates inflate coverage estimates and bias methylation calls.
Deduplication is essential for WGBS but must NOT be used for RRBS.

## Bismark Deduplication

Bismark provides its own deduplication tool optimized for bisulfite data, and this is what
the workflow runs (skipped with `--skip_dedup`):

```bash
deduplicate_bismark \
    --bam \
    --paired \
    bismark_out/sample_pe.bam
```

Bismark deduplication identifies duplicates by their alignment positions on
both strands, accounting for bisulfite conversion. This is preferred over
Picard for bisulfite data; Picard is not installed in the image.

The deduplicated BAM is then sorted and indexed with samtools and published as
`bismark/alignments/<sample>.sorted.bam` (+ `.bai`). The deduplication report goes to
`bismark/dedup_reports/`.

### Expected Duplication Rates

| Library Quality | Duplication Rate | Action |
|----------------|------------------|--------|
| Good | <20% | Proceed |
| Acceptable | 20-40% | Proceed with caution |
| Poor | 40-60% | Consider resequencing |
| Very poor | >60% | Library failed -- redo |

## Manual Alternative: Picard MarkDuplicates (not run by this workflow)

Picard is neither in the container image nor used by the workflow. If you align with
bwa-meth by hand (`02-bismark-alignment.md`), which has no built-in dedup, install Picard
separately and run:

```bash
picard MarkDuplicates \
    INPUT=sample_bwameth_sorted.bam \
    OUTPUT=sample_dedup.bam \
    METRICS_FILE=sample_dup_metrics.txt \
    REMOVE_DUPLICATES=true \
    VALIDATION_STRINGENCY=LENIENT \
    ASSUME_SORTED=true
```

## Optional Manual BAM Filtering (not run by this workflow)

The workflow hands MethylDackel the deduplicated, sorted BAM with no MAPQ or flag
filtering. If you want a filtered BAM, produce it yourself from the published
`<sample>.sorted.bam` and point MethylDackel at that instead:

```bash
samtools view -b -h \
    -q 10 \
    -F 1804 \
    -f 2 \
    bismark/alignments/sample.sorted.bam \
    | samtools sort -@ 4 -o sample_filtered.bam

samtools index sample_filtered.bam
```

### Filter Flag Explanation

| Flag | Binary | Meaning |
|------|--------|---------|
| `-q 10` | MAPQ >= 10 | Minimum mapping quality |
| `-F 4` | 0x4 | Remove unmapped reads |
| `-F 8` | 0x8 | Remove reads whose mate is unmapped |
| `-F 256` | 0x100 | Remove secondary alignments |
| `-F 512` | 0x200 | Remove reads failing QC |
| `-F 1024` | 0x400 | Remove PCR duplicates |
| `-f 2` | 0x2 | Keep only properly paired |

Combined: `-F 1804` = 4 + 8 + 256 + 512 + 1024, removing unmapped, mate-unmapped,
secondary, QC-fail and duplicate reads.

## RRBS: Skip Deduplication

For RRBS libraries, MspI digestion creates identical fragment starts at cut
sites. These are NOT PCR duplicates and must be retained. The trimming in this workflow is
WGBS-specific (`01-qc-trimming.md`), so RRBS is better processed by hand; if you run it
through the workflow anyway, at least set `--skip_dedup true`.

```bash
# Do NOT run deduplication for RRBS.
# To filter by hand instead, drop the duplicate bit from the mask (780 = 4+8+256+512):
samtools view -b -h -q 10 -F 780 -f 2 sample.bam \
    | samtools sort -@ 4 -o sample_filtered.bam
```

## Coverage Statistics

The workflow computes its coverage statistics from the CpG bedGraph, not from the BAM;
see `05-qc-metrics.md`. The commands below are alternatives you can run by hand on the
published BAM. `bedtools` is not in the container image — it ships with the conda
environment instead.

```bash
samtools depth -a bismark/alignments/sample.sorted.bam \
    | awk '{sum+=$3; n++} END {print "Mean coverage:", sum/n}'

samtools flagstat bismark/alignments/sample.sorted.bam > sample_flagstat.txt

# CpG-specific coverage: reads overlapping each CpG (last column), as a histogram.
# The BED goes in -a so the output is text; with a BAM in -a, bedtools writes BAM.
bedtools coverage \
    -a /ref/CpG_sites.bed \
    -b bismark/alignments/sample.sorted.bam \
    -counts \
    | awk '{print $NF}' | sort -n | uniq -c
```

# Methylation Calling with MethylDackel

MethylDackel (formerly PileOMeth) extracts per-CpG methylation levels from
bisulfite-aligned BAM files. This workflow then converts its bedGraph output to the
ENCODE bedMethyl layout (https://www.encodeproject.org/data-standards/wgbs/).

## M-bias Assessment (Run First)

The workflow runs `MethylDackel mbias` on every sample before extraction:

```bash
MethylDackel mbias \
    --CHG --CHH \
    /ref/genome/genome.fa \
    sample.sorted.bam \
    sample_mbias > sample_mbias_report.txt 2>&1
```

This produces SVG plots showing methylation level by read position, and the redirected
text output carries MethylDackel's suggested inclusion bounds. Both land in
`bismark/mbias/`. Look for:
- Elevated methylation at the 5' end of read 2 (end-repair artifact). The workflow already
  removes this at trimming with `--clip_R2 10`, so it should be gone
- Irregular methylation at read ends (adapter contamination)

## Methylation Extraction

This is exactly what the workflow runs (`--mergeContext` is dropped when
`--merge_context false`):

```bash
MethylDackel extract \
    --mergeContext \
    --CHG --CHH \
    --opref sample \
    /ref/genome/genome.fa \
    sample.sorted.bam
```

### Key Parameters

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--mergeContext` | Enabled by `--merge_context true` | Merge forward/reverse strand CpG data into one record |
| `--CHG --CHH` | Enabled | Also extract non-CpG methylation (useful for ESCs, neurons) |
| `--opref sample` | Output prefix | Names the three bedGraphs |

### Parameters the Workflow Deliberately Does Not Pass

| Parameter | Why not |
|-----------|---------|
| `--minDepth` | Extraction is unfiltered so the coverage statistics can see every covered site. The `--min_coverage` cut is applied when the bedMethyl files are written |
| `--maxDepth` | Not an option in MethylDackel 0.6.1. `getopt_long` rejects the long form and `extract` exits with a usage error. Do not pass it |
| `--nOT` / `--nOB` | The 5' end-repair bias of read 2 is removed at trimming (`--clip_R2 10`), so no read positions are excluded here |

### How `--OT`/`--nOT` Arguments Are Ordered

Getting this wrong silently trims the wrong end. From the MethylDackel 0.6.1 help,
`--OT A,B,C,D` means "include calls at positions from A through B on read #1 and C
through D on read #2". The `--nOT a,b,c,d` form is the complement and *excludes*:

| Slot | Meaning |
|------|---------|
| 1 | bases from the start (5') of read 1 |
| 2 | bases from the end (3') of read 1 |
| 3 | bases from the start (5') of read 2 |
| 4 | bases from the end (3') of read 2 |

So `--nOT 0,0,0,10` trims read 2's **3'** end, not its 5' end; read 2's 5' end is
`--nOT 0,0,10,0`. `--nOB 0,10,0,0` trims read 1's 3' end. If the M-bias plots still show
a problem, re-extract by hand with the bounds the mbias report suggests.

### Context-Specific Output Files

MethylDackel produces one bedGraph per context, named from `--opref` with an underscore.
They are published to `bismark/methylation/` and hold every covered site:

- `<sample>_CpG.bedGraph` -- CpG methylation (primary)
- `<sample>_CHG.bedGraph` -- CHG methylation (non-CpG)
- `<sample>_CHH.bedGraph` -- CHH methylation (non-CpG)

Columns are chrom, start, end, methylation percentage, methylated read count,
unmethylated read count, and the file starts with a `track` header line.

## Convert to bedMethyl Format

The workflow converts each context with this awk, then sorts, compresses and indexes:

```bash
for context in CpG CHG CHH; do
    awk -v min_cov=5 'BEGIN {OFS="\t"} !/^track/ {
        cov = $5 + $6
        if (cov == 0 || cov < min_cov) next
        score = (cov > 1000) ? 1000 : cov
        print $1, $2, $3, ".", score, ".", $2, $3, "0,0,0", cov, int(($5 / cov) * 100 + 0.5)
    }' sample_${context}.bedGraph \
        | sort -k1,1 -k2,2n \
        | bgzip > sample.${context}.bedMethyl.gz

    tabix -p bed sample.${context}.bedMethyl.gz
done
```

Four details matter and are easy to get wrong:

- **`!/^track/`** skips MethylDackel's header line. Without it the header is parsed as data.
- **Column 5 (score)** is the read count capped at 1000, per the ENCODE definition — not
  the methylation percentage scaled by 10.
- **Column 6 (strand)** is `.`. The MethylDackel bedGraph carries no strand, and with
  `--mergeContext` a record covers both strands of the CpG anyway.
- **Column 11** is an integer percentage (rounded), and column 10 is the coverage.

Sites below `min_cov` are dropped, which is the only place `--min_coverage` is applied.

## Per-Chromosome Extraction (Parallel, not run by this workflow)

For very large BAMs you can parallelize extraction by chromosome by hand:

```bash
for chr in $(samtools idxstats sample.sorted.bam | cut -f1 | grep -v '*'); do
    MethylDackel extract \
        --mergeContext \
        --opref "perchr/${chr}" \
        -r "${chr}" \
        /ref/genome/genome.fa \
        sample.sorted.bam &
done
wait

# Concatenate results, dropping the per-file track headers
cat perchr/*_CpG.bedGraph | awk '!/^track/' | sort -k1,1 -k2,2n > sample_CpG.bedGraph
```

## Global Methylation Summary (manual)

Compute genome-wide methylation statistics from a published bedGraph:

```bash
awk '!/^track/ {
    meth += $5; unmeth += $6; n++
} END {
    total = meth + unmeth;
    print "Covered sites:", n;
    print "Mean methylation:", (meth/total)*100 "%";
    print "Methylated reads:", meth;
    print "Unmethylated reads:", unmeth
}' sample_CpG.bedGraph
```

Expected values for mammalian somatic tissue:
- Global CpG methylation: 70-85%
- CpG islands: 5-15% (mostly unmethylated)
- Gene bodies: 60-80%
- Intergenic: 75-90%

# WGBS QC Metrics and Conversion Rate Assessment

Quality control for WGBS requires bisulfite-specific metrics beyond
standard alignment QC. The most critical metric, the bisulfite conversion rate, is not
computed by this workflow — the sections below say which commands you have to run
yourself.

## What the Workflow Produces

| Output | Contents |
|--------|----------|
| `bismark/alignments/*_PE_report.txt` | Mapping efficiency, and the percentage of C methylated in CpG/CHG/CHH context |
| `bismark/dedup_reports/*.deduplication_report.txt` | Duplicate count and rate |
| `bismark/mbias/<sample>_mbias_*.svg`, `<sample>_mbias_report.txt` | M-bias plots and MethylDackel's suggested inclusion bounds |
| `coverage/<sample>.coverage_stats.txt` | Covered CpGs, their mean coverage, and the fraction reaching >=5x, >=10x and `--min_coverage`x |
| `multiqc/multiqc_report.html` | Aggregated report |

## Coverage Statistics

`COVERAGE_STATS` parses the unfiltered `<sample>_CpG.bedGraph`, so it sees every covered
site rather than only the ones that pass `--min_coverage`:

```bash
awk -v min_cov=5 '!/^track/ {
    cov = $5 + $6; sum += cov; n++;
    if (cov >= 5)       c5++;
    if (cov >= 10)      c10++;
    if (cov >= min_cov) cmin++
} END {
    if (n == 0) { print "Covered CpGs: 0"; exit }
    printf "Covered CpGs (>=1x): %d\n", n;
    printf "Mean coverage of covered CpGs: %.1f\n", sum/n;
    printf "Covered CpGs >=5x: %d (%.1f%%)\n", c5, c5/n*100;
    printf "Covered CpGs >=10x: %d (%.1f%%)\n", c10, c10/n*100;
    printf "Covered CpGs >=%dx (--min_coverage, kept in bedMethyl): %d (%.1f%%)\n", min_cov, cmin, cmin/n*100
}' sample_CpG.bedGraph > sample.coverage_stats.txt
```

`min_cov` is `--min_coverage`, so the five-line file looks like this at the default of 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%)
```

Three things to keep straight when reporting these numbers:

- The 5x and 10x thresholds are fixed in the workflow. They do not follow `--min_coverage`.
- The last line is the one that does, and it is the count that reached the bedMethyl
  files. At the default it duplicates the >=5x line; at any other `--min_coverage` it does
  not.
- The percentages are of *covered* CpGs. A CpG that received zero reads is not in the
  bedGraph and is not counted, so this is not the genome-wide CpG completeness.

With `--merge_context false` the records are per cytosine, not per CpG, and the counts
roughly double.

## Bisulfite Conversion Rate (Lambda Spike-in, not run by this workflow)

Lambda phage DNA is fully unmethylated. Any methylation detected on lambda
represents incomplete bisulfite conversion. Run this yourself, against the trimmed reads,
if a spike-in was included:

```bash
# Align to lambda genome
bismark \
    --genome /ref/lambda/ \
    --bowtie2 \
    --parallel 2 \
    -1 sample_R1_val_1.fq.gz \
    -2 sample_R2_val_2.fq.gz \
    --output_dir lambda_out/ \
    --unmapped

# Extract methylation from lambda alignments
MethylDackel extract \
    --mergeContext \
    --minDepth 1 \
    /ref/lambda/genome.fa \
    lambda_out/sample_pe.bam

# Calculate conversion rate
awk '!/^track/ {meth+=$5; unmeth+=$6} END {
    total=meth+unmeth;
    conv=(unmeth/total)*100;
    print "Conversion rate: " conv "%";
    print "Unconverted (false methylation): " (meth/total)*100 "%"
}' lambda_CpG.bedGraph
```

### Conversion Rate Thresholds

| Rate | Status | Action |
|------|--------|--------|
| >99.5% | Excellent | Proceed |
| 99.0-99.5% | Acceptable | Proceed with note |
| 98.0-99.0% | Warning | May inflate methylation estimates |
| <98.0% | Fail | Do NOT use this library |

## Non-CpG Methylation as Conversion Proxy

If no spike-in is available, use CHH methylation as a proxy. The Bismark alignment report
already carries "C methylated in CHH context"; the same number can be recomputed from the
published CHH bedGraph:

```bash
awk '!/^track/ {meth+=$5; unmeth+=$6} END {
    print "CHH methylation: " (meth/(meth+unmeth))*100 "%"
}' sample_CHH.bedGraph
```

In somatic tissue, CHH methylation should be <1%. Higher values suggest
incomplete conversion. Exception: embryonic stem cells and neurons can have
genuine non-CpG methylation (2-5%), which is why this is a proxy rather than a
measurement.

## Additional Manual Coverage Checks

```bash
# Genome-wide coverage distribution from the published BAM
samtools depth -a bismark/alignments/sample.sorted.bam | \
    awk '{cov[$3]++} END {for (c in cov) print c, cov[c]}' | \
    sort -k1,1n > coverage_distribution.txt

# Mean and median coverage (needs gawk for asort)
samtools depth -a bismark/alignments/sample.sorted.bam | \
    awk '{sum+=$3; n++; a[n]=$3; if ($3>=5) sum5++} END {
        asort(a);
        print "Mean:", sum/n;
        print "Median:", a[int(n/2)];
        print "Total bases:", n;
        print "Bases >=5x:", sum5/n*100 "%"
    }'
```

## Mapping Statistics

The mapping efficiency in the Bismark `*_PE_report.txt` is the primary number. For a
flag-level breakdown of the final BAM:

```bash
samtools flagstat bismark/alignments/sample.sorted.bam > flagstat.txt
```

Key values to extract:
- Total reads in the final (deduplicated) BAM
- Mapped reads (Bismark mapping efficiency >70% expected)
- Properly paired (expect >95% of mapped)
- Duplication rate (from the deduplication report)

## MultiQC Report

The workflow runs:

```bash
multiqc --title "ENCODE WGBS Pipeline" --filename multiqc_report --force .
```

`--filename` is required: without it MultiQC 1.21 derives the file name from `--title`
and writes `ENCODE-WGBS-Pipeline_multiqc_report.html`, which would not match the declared
output.

The inputs collected are the raw-read FastQC reports, the Trim Galore trimming reports,
the FastQC reports for the trimmed reads, the Bismark `*_PE_report.txt` files and the
`*.deduplication_report.txt` files, so the duplication rate appears in the report as well
as in `bismark/dedup_reports/`. With `--skip_dedup` there are no deduplication reports to
collect. Picard metrics do not exist for this pipeline.

## Summary QC Table Format

Generate a per-sample summary for reporting:

```bash
echo -e "Sample\tTotal_Reads\tMapping_Rate\tDedup_Rate\tConversion\tMean_CpG_Cov\tCpGs_5x"
echo -e "${SAMPLE}\t${TOTAL}\t${MAP_RATE}\t${DEDUP_RATE}\t${CONV_RATE}\t${MEAN_COV}\t${CPGS_5X}"
```

`CONV_RATE` has to come from the manual check above; everything else is in the published
reports.

# WGBS Pipeline — Literature References

**Last updated:** 2026-03-07
**Purpose:** Reference catalog for the pipeline-wgbs skill — papers defining ENCODE WGBS processing standards, bisulfite-aware alignment, methylation extraction, and quality assessment for whole-genome bisulfite sequencing.

---

## DNA Methylation Method & Standards

---

### Lister et al. 2009 — Human DNA methylomes at base resolution

- **Citation:** Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nott JR, Lee L, Ye Z, Ngo QM, Edsall L, Antosiewicz-Bourget J, Stewart R, Ruotti V, Millar AH, Thomson JA, Ren B, Ecker JR. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature, 462(7271):315-322, 2009.
- **DOI:** [10.1038/nature08514](https://doi.org/10.1038/nature08514)
- **PMID:** 19829295 | **PMC:** PMC2857523
- **Citations:** ~5,000
- **Key findings:** First base-resolution whole-genome DNA methylation maps for human (H1 ESCs and IMR90 fibroblasts). Demonstrated that ~25% of all methylated cytosines in embryonic stem cells are in non-CpG context (CHG, CHH), a phenomenon largely absent in differentiated cells. Identified widespread methylation differences between cell types, including large partially methylated domains (PMDs) covering ~40% of the genome in fibroblasts that correlate with gene repression and late-replicating chromatin. Established WGBS as the gold standard for comprehensive methylation profiling at single-nucleotide resolution, requiring ~30x genome coverage for reliable quantification.

---

### Schultz et al. 2015 — Human body epigenome maps

- **Citation:** Schultz MD, He Y, Whitaker JW, Hariharan M, Mukamel EA, Leung D, Rajagopal N, Nery JR, Urich MA, Chen H, Lin S, Lin Y, Jung I, Schmitt AD, Selvaraj S, Ren B, Sejnowski TJ, Wang W, Ecker JR. Human body epigenome maps reveal noncanonical DNA methylation variation. Nature, 523(7559):212-216, 2015.
- **DOI:** [10.1038/nature14248](https://doi.org/10.1038/nature14248)
- **PMID:** 26030523 | **PMC:** PMC4539777
- **Citations:** ~1,500
- **Key findings:** Generated WGBS methylomes for 18 human tissues as part of the Roadmap Epigenomics Project. Identified tissue-specific differentially methylated regions (DMRs) and described three classes of genomic methylation: highly methylated regions (HMRs, >80%), partially methylated domains (PMDs, 50-80%), and unmethylated regions (UMRs, <20%). PMDs are associated with late-replicating heterochromatin and expand during aging and in cancer. UMRs at CpG islands mark active promoters, while tissue-specific UMRs often mark active enhancers. Established the analytical framework for WGBS data interpretation used in this pipeline.

---

### ENCODE Project Consortium 2020 — Expanded encyclopaedias of DNA elements

- **Citation:** ENCODE Project Consortium et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature, 583(7818):699-710, 2020.
- **DOI:** [10.1038/s41586-020-2493-4](https://doi.org/10.1038/s41586-020-2493-4)
- **PMID:** 32728249 | **PMC:** PMC7410828
- **Citations:** ~2,500
- **Key findings:** ENCODE Phase 3 paper establishing DNA methylation as a core epigenomic layer alongside histone modifications and chromatin accessibility. WGBS data is used to identify methylation valleys (large UMRs) at key developmental genes and to classify candidate cis-Regulatory Elements (cCREs) by their methylation status. Established the ENCODE WGBS pipeline standards: bisulfite conversion rate ≥98%, minimum 30x CpG coverage for reliable DMR calling, and GRCh38/mm10 as reference assemblies.

---

### Hitz et al. 2023 — The ENCODE Uniform Analysis Pipelines

- **Citation:** Hitz BC, Lee JW, Jolanki O, et al. The ENCODE Uniform Analysis Pipelines. Research Square (preprint), 2023.
- **DOI:** [10.21203/rs.3.rs-311127/v1](https://doi.org/10.21203/rs.3.rs-311127/v1)
- **Citations:** ~84
- **Key findings:** Documents the official ENCODE WGBS pipeline specification: Bismark for bisulfite-aware alignment, deduplication with Picard MarkDuplicates, MethylDackel for methylation extraction, and lambda phage spike-in for conversion rate quality control. Specifies bedMethyl as the standard output format and defines QC thresholds: bisulfite conversion ≥98%, mapping rate >70%, and sufficient CpG coverage depth for downstream DMR analysis. Pipeline infrastructure uses Docker containers and WDL for reproducible execution.

---

### Foox et al. 2021 — Performance assessment of DNA sequencing platforms

- **Citation:** Foox J, Nordlund J, Lalonde E, Lim S, Deschamps-Francoeur G, Trognitz F, Holmes M, Close MC, et al. Performance assessment of DNA sequencing platforms in the ABRF Next-Generation Sequencing Study. Nature Biotechnology, 39(9):1129-1140, 2021.
- **DOI:** [10.1038/s41587-021-01049-5](https://doi.org/10.1038/s41587-021-01049-5)
- **PMID:** 34504353
- **Citations:** ~200
- **Key findings:** Multi-site benchmarking of sequencing platforms for WGBS applications. Demonstrated that bisulfite conversion rate is the most critical QC parameter, with rates below 99% introducing thousands of false-positive methylation calls across the genome. Platform-specific biases in GC-rich regions affect CpG island coverage. Supports the ENCODE requirement for ≥98% conversion rate and provides the evidence base for this pipeline's QC thresholds.

---

## Core Pipeline Tools

---

### Krueger & Andrews 2011 — Bismark

- **Citation:** Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics, 27(11):1571-1572, 2011.
- **DOI:** [10.1093/bioinformatics/btr167](https://doi.org/10.1093/bioinformatics/btr167)
- **PMID:** 21493656 | **PMC:** PMC3102221
- **Citations:** ~4,000
- **Key findings:** Introduced Bismark, the most widely used bisulfite-aware aligner. Bismark performs in silico C-to-T conversion of both the reference genome and sequencing reads, then aligns using Bowtie2 as the backend. After alignment, it determines the original methylation state of each cytosine by comparing the unconverted read to the reference. Handles all three cytosine contexts (CpG, CHG, CHH) and generates comprehensive methylation reports. The ENCODE pipeline uses Bismark as the default aligner for WGBS data, leveraging its --paired mode with --no_overlap to avoid double-counting methylation in overlapping paired-end reads.

---

### Pedersen et al. 2014 — bwa-meth

- **Citation:** Pedersen BS, Eyring K, De S, Yang IV, Schwartz DA. Fast and accurate alignment of long bisulfite-seq reads. arXiv:1401.1129, 2014.
- **DOI:** [10.48550/arXiv.1401.1129](https://doi.org/10.48550/arXiv.1401.1129)
- **Citations:** ~200
- **Key findings:** Introduced bwa-meth, an alternative bisulfite-aware aligner that wraps BWA-MEM for bisulfite-converted reads. Rather than creating four converted reference copies (like Bismark), bwa-meth converts only the C-to-T strand and the G-to-A strand, achieving faster alignment while maintaining accuracy. Particularly advantageous for longer reads (>100 bp) where BWA-MEM's local alignment outperforms Bowtie2's end-to-end mode. Included in this pipeline as an alternative aligner option (--aligner bwameth) for users who prefer faster processing or have long-read bisulfite data.

---

### MethylDackel — Methylation extraction from BAM files

- **Citation:** Ryan DP. MethylDackel: A (mostly) universal methylation extractor for BS-seq experiments. GitHub, 2023.
- **URL:** [https://github.com/dpryan79/MethylDackel](https://github.com/dpryan79/MethylDackel)
- **Key findings:** MethylDackel (formerly PileOMeth) extracts per-base methylation metrics from BAM files produced by any bisulfite-aware aligner. Key features used in this pipeline: --mergeContext flag to merge complementary CpG dinucleotide strands, --minDepth to filter low-coverage positions, M-bias plot generation for identifying end-repair artifacts, and --OT/--OB position-based trimming to remove biased read positions. MethylDackel outputs in bedMethyl format, the ENCODE standard for methylation data. Replaces Bismark's built-in methylation extractor with improved speed and flexibility.

---

## Quality Control

---

### Ewels et al. 2016 — MultiQC

- **Citation:** Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32(19):3047-3048, 2016.
- **DOI:** [10.1093/bioinformatics/btw354](https://doi.org/10.1093/bioinformatics/btw354)
- **PMID:** 27312411 | **Citations:** ~6,800
- **WGBS role:** Aggregates QC metrics from FastQC, Bismark alignment reports, Bismark deduplication reports, and MethylDackel M-bias plots into a unified HTML report. Provides batch-level assessment of bisulfite conversion rates, mapping rates, and duplication across all samples.

---

## Supplementary Tools (Non-WGBS-Specific)

See pipeline-chipseq/references/literature.md for detailed descriptions of shared tools.

---

### Li et al. 2009 — SAMtools

- **DOI:** [10.1093/bioinformatics/btp352](https://doi.org/10.1093/bioinformatics/btp352) | **PMID:** 19505943 | **Citations:** ~53,700
- **WGBS role:** BAM sorting, indexing, and filtering. Used for coordinate sorting after Bismark alignment and computing alignment statistics. Provides samtools flagstat for mapping rate assessment.

---

### Broad Institute — Picard MarkDuplicates

- **URL:** [https://broadinstitute.github.io/picard/](https://broadinstitute.github.io/picard/)
- **WGBS role:** PCR duplicate marking and removal. Critical for WGBS because PCR amplification of bisulfite-converted DNA is biased toward unmethylated molecules (which align better after conversion), making duplicates non-random. High duplication rates (>30%) indicate insufficient library complexity and can bias methylation estimates.

---

### Andrews 2010 — FastQC

- **URL:** [https://www.bioinformatics.babraham.ac.uk/projects/fastqc/](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
- **WGBS role:** Raw read quality assessment. Bisulfite-converted reads have characteristic quality profiles: reduced base diversity (most Cs converted to Ts), asymmetric per-base composition, and lower quality scores than standard WGS. These patterns are expected and not cause for concern in WGBS data.

---

### Martin 2011 — Cutadapt (basis for Trim Galore)

- **DOI:** [10.14806/ej.17.1.200](https://doi.org/10.14806/ej.17.1.200) | **Citations:** ~13,000
- **WGBS role:** Adapter trimming via Trim Galore wrapper. Particularly important for WGBS because bisulfite conversion reduces sequence complexity, making adapter contamination harder to detect. Trim Galore includes bisulfite-specific adapter detection and optional RRBS-aware trimming (--rrbs flag) for reduced representation data.

---

### Amemiya et al. 2019 — ENCODE Blacklist

- **DOI:** [10.1038/s41598-019-45839-z](https://doi.org/10.1038/s41598-019-45839-z) | **PMID:** 31249361 | **Citations:** ~1,372
- **WGBS role:** Blacklist v2 filtering removes artifact-prone regions before downstream analysis. In WGBS, blacklisted regions can show aberrant methylation estimates due to alignment artifacts in repetitive sequences. Applied after methylation calling to filter bedMethyl output.

---

## Additional Relevant Papers

---

### Hansen et al. 2012 — BSmooth: large-scale WGBS smoothing

- **Citation:** Hansen KD, Langmead B, Irizarry RA. BSmooth: from whole genome bisulfite sequencing reads to differentially methylated regions. Genome Biology, 13(10):R83, 2012.
- **DOI:** [10.1186/gb-2012-13-10-r83](https://doi.org/10.1186/gb-2012-13-10-r83)
- **PMID:** 23034175 | **PMC:** PMC3491411
- **Citations:** ~1,200
- **Key findings:** Introduced local-likelihood smoothing for WGBS data, demonstrating that single-CpG methylation estimates from WGBS are inherently noisy even at high coverage, and that smoothing across neighboring CpGs substantially improves differential methylation detection. Established the analytical framework for identifying differentially methylated regions (DMRs) from WGBS data, implemented in the bsseq R/Bioconductor package. Relevant for downstream analysis of pipeline output.

---

### Peters et al. 2015 — DMRcate: differentially methylated region analysis

- **Citation:** Peters TJ, Buckley MJ, Statham AL, Pidsley R, Samaras K, Lord RV, Clark SJ, Molloy PL. De novo identification of differentially methylated regions in the human genome. Epigenetics & Chromatin, 8:6, 2015.
- **DOI:** [10.1186/1756-8935-8-6](https://doi.org/10.1186/1756-8935-8-6)
- **PMID:** 25972926 | **PMC:** PMC4429357
- **Citations:** ~900
- **Key findings:** Introduced DMRcate, a method for de novo identification of differentially methylated regions using kernel smoothing and bump-hunting. Works with both array (450K/EPIC) and WGBS data, providing a unified framework for DMR detection. Relevant for downstream analysis of bedMethyl output from this pipeline, particularly when comparing methylation across conditions or tissues.