Bio Atac Seq Consensus Peakset

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.2+ (featureCounts; --countReadPairs requires >= 2.0.2), 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 to count reads against for 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 (per-sample membership columns emitted by default)
• 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

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.

#!/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+1, $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 a faster approximation at scale, bedtools cluster + per-cluster top-significance selection is tempting, but it is NOT equivalent to greedy iterative overlap: cluster groups transitively-overlapping peaks (a chain can span well beyond 501 bp), so keeping one peak per cluster discards non-overlapping peaks the greedy algorithm would retain. Use it only as an approximation.

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 ... (which emits per-sample membership columns by default) and 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 median width as approximation.

IDR-filtered union -- Computational cost; threshold mismatch

Trigger: Running IDR for each rep pair, then unioning IDR-passed peaks.

Mechanism: IDR thresholds differ (true reps 0.05; pseudoreps 0.10). Mixing both into a union biases the consensus toward looser pseudorep peaks unless filtering carefully.

Fix: Decide upfront -- use only true-rep IDR (stricter, fewer peaks) OR pseudorep IDR (looser, more peaks). Don't mix.

Decision Tree by Goal