bio-chipseq-differential-binding

What this skill does

## Version Compatibility

Reference examples tested with: DiffBind 3.20+, DESeq2 1.42+, edgeR 4.0+, csaw 1.36+, PyDESeq2 0.4+, NormR 1.28+, MAnorm2 1.2+, ChIPseqSpikeInFree 1.6+.

DiffBind 3.0+ changed defaults: `summits=200` (was FALSE), `dba.normalize()` now required, blacklist filtering on by default, full library size normalization replaces reads-in-peaks. Always run `packageVersion('DiffBind')` and inspect `dba.normalize(obj, bRetrieve=TRUE)` to confirm what was applied.

# Differential ChIP-seq Binding

**"Compare protein-DNA binding between experimental conditions"** -> Identify regions where IP signal changes significantly, accounting for sequencing depth, composition bias, trended biases, and global shifts that confound naive normalization.

- R (BAM + peaks): `DiffBind::dba()` -> `dba.count()` -> `dba.normalize()` -> `dba.analyze()`
- R (count matrix): `DESeq2::DESeq()` or `edgeR::glmQLFTest()` on a peaks-by-samples matrix
- R (windows-based, global-shift-robust): `csaw::windowCounts()` -> `csaw::normFactors()` -> `edgeR::glmQLFTest()`
- R (control-aware): `normr::diffR(chip1.bam, chip2.bam, genome)` joint binomial mixture
- Python (count matrix): `pydeseq2.DeseqDataSet()`

Choice of normalization matters more than choice of test statistic (RLE vs TMM vs csaw bin-TMM on the same reference reads produce nearly identical results). Choose by which of the three normalization problems applies.

## The Three Distinct Normalization Problems

| Problem | Symptom on MA plot | Cause | Fix |
|---------|--------------------|-------|-----|
| **Composition bias** | Loess shifts off y=0 systematically | Few high-signal peaks dominate read counts; small fold changes look large or inverted | TMM on background 10 kb bins (csaw / DiffBind `background=TRUE`); NOT reads-in-peaks |
| **Trended bias (intensity-dependent)** | Loess curve sweeps from + to - across abundance | Library-prep efficiency varies with fragment abundance | Non-linear loess offsets (csaw `normalizeOffsets`); use cautiously — can over-normalize biology |
| **Global shift (treatment changes most peaks)** | Loess entirely shifted off y=0; mean log2FC ≠ 0 | Drug/perturbation changes the genome-wide level of binding (HDACi, BETi, EZH2i, target KD) | Spike-in scaling (ChIP-Rx); no algorithmic fix works |

**Why this matters:** Most published ChIP-seq differential analyses default to RLE/TMM on reads-in-peaks, which assumes "most peaks are unchanged." For HDAC inhibitors, BET inhibitors, EZH2 inhibitors, or any large dosage / target-knockdown experiment, this assumption is violated. The algorithm forces the median log2FC to zero, hides the real effect, and amplifies noise around the new "zero." The result can have the wrong sign.

**Diagnostic:** Plot MA loess on differential results. A loess curve that sweeps abundance indicates trended bias. A uniformly shifted loess indicates a global shift. Both patterns together indicate normalization is failing in two ways simultaneously.

## Algorithmic Taxonomy

| Tool | Treats | Statistical model | Strength | Fails when |
|------|--------|-------------------|----------|------------|
| **DiffBind 3.20+** | Consensus peaks summit ± 250 bp (default) | DESeq2 or edgeR backend | Mature; integrated counting/normalization; spike-in support; blacklist filter on | Default `summits=200` recenters peaks (wrong for broad marks); `DBA_NORM_LIB` is conservative but misses global shifts unless `background=TRUE` |
| **DESeq2 (direct)** | Predefined peaks | NB GLM | Familiar; transparent | RLE on reads-in-peaks fails for global shifts |
| **edgeR (direct)** | Predefined peaks | NB GLM with TMM; quasi-likelihood F-test | Cleaner small-sample inference; QL-F controls type-I error | TMM on peak counts unstable if peaks globally shifting |
| **csaw** (Lun & Smyth 2016) | Sliding windows (typically 150 bp width, 50 bp shift) | edgeR QL-F | Gold-standard for global shifts; bin-TMM composition bias; loess for trended biases | Slower; requires BAMs not count matrix; window-merge step adds complexity |
| **NormR** (Helmuth 2016) | Genomic bins | Binomial mixture (background + enriched) | Control-aware; identifies enrichment/depletion/background simultaneously | Bin-level not peak-level; older codebase; less integration with downstream tools |
| **MAnorm2** (Tu 2021) | Peaks | Hierarchical model with mean-variance trend | Designed for cross-condition with replicates | Less widely adopted; sparse maintenance |
| **SpikChIP** (Blanco 2021) | Peaks | Spike-in-aware | Multi-sample spike-in comparison | Niche; specific spike-in protocol assumed |
| **SpikeFlow** (2024) | End-to-end | Snakemake wrapper around DiffBind + spike-in | Automated; multiple normalization options | Wrapper; inherits experimental-design errors upstream |
| **ChIPComp** (Chen 2015) | Peaks | Joint Poisson with input | Control-aware | Older; less maintained |
| **ChIPseqSpikeInFree** (Jin 2020) | Peaks (post-hoc) | Distribution-shape inference | Detects global shift WITHOUT spike-in | Post-hoc heuristic only; not definitive; sanity check |

## Decision Tree: Choosing Normalization

| Scenario | Recommended | Tool / parameter |
|----------|-------------|------------------|
| Standard TF ChIP, local changes expected, balanced gain/loss | RLE on reads-in-peaks | DESeq2 default; DiffBind `DBA_NORM_RLE` |
| Histone marks, broad domains, local changes | TMM on background 10 kb bins | csaw `normFactors`; DiffBind `background=TRUE` |
| HDAC / BET / EZH2 inhibitor (global change) | Spike-in scaling | DiffBind `spikein=TRUE`; SpikeFlow; manual `sizeFactors()` from spike reads |
| Dosage titration, cell-cycle synchronization | Spike-in scaling | Same |
| ChIP target knockdown / degron | Spike-in or matched-input control subtraction | Spike-in preferred; bamCompare log2 ratio next |
| CUT&RUN/CUT&Tag standard | E. coli spike-in (carryover) | DiffBind custom scaling; see cut-and-run-tag |
| No spike-in available; suspect global shift | ChIPseqSpikeInFree | Post-hoc distribution-shape inference |
| Suspected trended (abundance-dependent) bias | Non-linear loess | csaw `normalizeOffsets` |
| Genome-wide enrichment/depletion analysis | NormR | Binomial mixture |
| Many conditions, large peak set | edgeR QL-F or DiffBind+edgeR | Better type-I control than DESeq2 Wald |

## Spike-In Scaling Factor Calculation

ChIP-Rx (Orlando 2014) uses Drosophila chromatin spike-in added at fixed concentration BEFORE IP. Egan 2016 protocol: 50,000 Drosophila S2 nuclei per 5M target cells.

**RRPM (reference-adjusted reads per million):**

```
scale_factor_i = min(N_spike_sample) / N_spike_sample_i
```

Apply to read counts pre-test, OR pass as `sizeFactors()` to DESeq2 / `normFactors()` to edgeR / a numeric library-size vector to DiffBind's `dba.normalize(..., library=...)`:

```r
spike_in_reads <- c(120000, 145000, 110000, 95000)  # per-sample Drosophila read count
sample_names <- c('ctrl_1', 'ctrl_2', 'treat_1', 'treat_2')
scale_factors <- min(spike_in_reads) / spike_in_reads
names(scale_factors) <- sample_names

# DESeq2 with spike-in size factors
sizeFactors(dds) <- 1 / scale_factors  # DESeq2 expects inverse convention

# Or directly via DiffBind 3.x
dba_obj <- dba.normalize(dba_obj, spikein = TRUE,
                          library = DBA_LIBSIZE_FULL,
                          normalize = DBA_NORM_LIB)
```

**Rx-Input variant** (Fursova 2019): additionally scale by input spike-in to correct IP efficiency variation.

**Internal-control sanity check (Patel L, Cao Y, Mendenhall EM, Benner C, Goren A 2024 *Nat Biotechnol* 42:1343):** After spike-in normalization, blacklist regions and constitutive housekeeping sites (U6 promoter, rRNA processing factors that are stable) should show no signal change. If they do, the normalization is broken — common causes:
- Spike-in scaling applied to peak counts instead of read counts
- Spike-in reads not deduplicated before scaling
- Spike-in genome not filtered for high-mapq before scaling
- Spike-in satur