bio-atac-seq-consensus-peakset

What this skill does

## Version Compatibility

Reference examples tested with: bedtools 2.31+, samtools 1.19+, BEDOPS 2.4.41+, GenomicRanges 1.54+, DiffBind 3.12+, Subread 2.0+ (featureCounts), pybedtools 0.10+.

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws unexpected errors, introspect the installed package and adapt rather than retrying.

# Consensus Peakset Construction

**"Build a single peakset I can count reads against across all my samples"** -> Combine per-replicate or per-condition peak calls into a non-redundant, fixed-width set of regions. The strategy chosen drives FDR calibration, peak-width fairness, and reproducibility downstream.

- CLI: `bedtools merge` (simple union) and `bedtools multiinter -j` (per-sample membership)
- CLI: Corces 2018 iterative overlap removal (custom shell)
- R: `DiffBind::dba.count(summits=250)` (built-in fixed-width)
- Python: `pybedtools` for programmatic merging

The peakset choice is rarely default-correct. Wrong width or wrong overlap rule propagates to every downstream analysis (differential, motif, footprint, ML).

## Why a Consensus Peakset Matters

ATAC peaks vary in width across replicates: same regulatory element might be called 200 bp in rep1 and 800 bp in rep2 because of stochastic Tn5 cuts at edges. Counting reads in different-width intervals confounds peak width with biological signal. A fixed-width consensus avoids this.

For ENCODE-style differential analysis: ALL samples must be counted against the SAME peak coordinates; otherwise the count matrix is non-rectangular and statistical models are misspecified.

## Strategy Taxonomy

| Strategy | Implementation | Width | When to use | Fails when |
|----------|---------------|-------|-------------|------------|
| Naive union | `bedtools merge` of all peaks | Variable, tends wide | Quick exploratory; never for differential | Width inflation drives spurious differential |
| Naive intersection | `bedtools multiinter` requiring all samples | Variable | High-stringency reproducibility | Loses real condition-specific peaks |
| Majority-rule overlap | `multiinter` requiring >= n/2 samples | Variable | Balance; DiffBind default with `minOverlap` | Width still varies; counts are width-biased |
| Iterative overlap removal (Corces 2018) | Sort by significance, greedily keep non-overlapping at fixed width | 501 bp fixed | ML features; cross-study comparison; modern ATAC standard | Loses sub-501bp resolution; overweights high-significance peaks |
| Summit-centered fixed width (DiffBind) | `dba.count(summits=250)` re-centers all peaks on summit +/- 250 bp | 501 bp fixed | Default for DiffBind workflows; integrates with replicate counts | Requires summit info (MACS narrowPeak); broad peaks lose width info |
| IDR-filtered union | Union of IDR-passed peaks across rep pairs | Variable | ENCODE pipeline-compliant; reproducibility-aware | Requires running IDR per pair; computationally heavier |
| Per-condition union, then global union | Each group consensus separately, then merge | Variable | Different cell types / strong condition shift | Same width issues as naive union |
| Width-controlled extension | Extend each peak to median width centered on midpoint | User-set | Quick fixed-width without summit info | Midpoint != summit; can shift biology |

Methodology evolves; verify against current ENCODE 4 ATAC standards (encodeproject.org/atac-seq) and the Corces 2018 iterative overlap algorithm before locking pipelines.

## Iterative Overlap Removal (Corces 2018)

This is the modern standard for fixed-width consensus peaksets used in cross-study comparison and machine-learning feature matrices.

**Algorithm:**
1. Pool all peaks from all samples; re-center each on its summit; extend +/- 250 bp -> 501 bp fixed-width peaks.
2. Sort by significance (narrowPeak column 7, signalValue, descending).
3. Walk down the sorted list. Keep each peak if it does not overlap any previously kept peak. Drop if overlap.
4. Output the kept peaks as the consensus.

**Goal:** Produce a non-overlapping, fixed-width peakset weighted toward strongest evidence.

**Approach:** Greedy non-overlap on summit-centered fixed-width peaks ranked by signalValue.

```bash
#!/bin/bash
# Corces 2018 iterative overlap removal
PEAKS_DIR=peaks/per_sample
GENOME_SIZES=hg38.chrom.sizes
WIDTH_HALF=250                                       # 501 bp total width

# 1. Pool all peaks, re-center on summit, extend
awk -v w=$WIDTH_HALF 'BEGIN{OFS="\t"}
    {summit=$2+$10; print $1, summit-w, summit+w, $4, $7, $6}' \
    $PEAKS_DIR/*.narrowPeak | \
    awk '$2 >= 0' | \
    bedtools slop -i - -g $GENOME_SIZES -b 0 | \
    sort -k1,1 -k2,2n > pooled_recentered.bed

# 2. Sort by signalValue descending (column 5 in our BED -- which was column 7 of narrowPeak)
sort -k5,5gr pooled_recentered.bed > pooled_by_sig.bed

# 3. Iterative greedy non-overlap (process in significance order)
python3 - <<'EOF'
import sys
kept = []
with open('pooled_by_sig.bed') as f:
    for line in f:
        chrom, start, end, name, sig, strand = line.strip().split('\t')[:6]
        start, end = int(start), int(end)
        overlap = any(c == chrom and not (end <= s or start >= e) for c, s, e in kept)
        if not overlap:
            kept.append((chrom, start, end))

with open('consensus_iterative.bed', 'w') as f:
    for c, s, e in sorted(kept):
        f.write(f'{c}\t{s}\t{e}\n')
EOF

sort -k1,1 -k2,2n consensus_iterative.bed > consensus_final.bed
echo "Consensus peakset: $(wc -l < consensus_final.bed) fixed-width 501bp peaks"
```

For better performance at scale, replace the Python loop with `bedtools cluster` followed by per-cluster top-significance selection.

## Per-Strategy Failure Modes

### Naive union -- Width-driven differential

**Trigger:** Using `bedtools merge` of all per-rep peaks; counting reads in merged intervals.

**Mechanism:** A peak appearing as 200 bp in one rep but 800 bp in another merges to 800 bp in the union. Read count in 800 bp interval is biased high; differential analysis flags it as condition-specific even when it's just width difference.

**Symptom:** Top differential peaks track peak width (mean width different by 100+ bp between conditions).

**Fix:** Use fixed-width strategy (Corces iterative or DiffBind `summits=250`). Never use merged variable-width peaks for differential.

### Intersection (peak in all reps) -- Loses condition-specific biology

**Trigger:** Using `bedtools multiinter -i ... -j -intervals` requiring all samples.

**Mechanism:** Peaks present only in one condition fail intersection requirement and are excluded from the consensus, even though they are the biology of interest.

**Symptom:** Differential analysis returns near-zero significant peaks; the condition-specific peaks were filtered out before testing.

**Fix:** Use per-condition consensus, then union of consensus. Or majority rule with per-condition minOverlap.

### Majority rule (DiffBind default-ish) -- Borderline peaks dropped

**Trigger:** `dba.count(minOverlap=ceiling(N/2))` where N = total replicates and one condition has fewer reps.

**Mechanism:** A peak in 2/2 reps of cond1 but 0/3 reps of cond2 is in 2/5 = 40% < 50% -> dropped.

**Fix:** Compute consensus per condition first (e.g. peak in >= 2/3 reps), then union across conditions. This preserves condition-specific peaks.

### Width-controlled extension -- Shifts off summit

**Trigger:** Extending peak to fixed width using midpoint (`(start+end)/2`).

**Mechanism:** Midpoint of the called peak is rarely at the actual summit; particularly for asymmetric peaks (TSS-flanking) or wide broadPeaks.

**Fix:** Use summit position from narrowPeak column 10 (`start + summit_offset`). If summit unavailable (broadPeak), use peak start + half me