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

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	Matches the Corces 501 bp convention (note: dba.count default is summits=200 -> 401 bp); 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