Bio Chipseq Differential Binding

Version Compatibility

Reference examples tested with: DiffBind 3.20+, DESeq2 1.42+, edgeR 4.0+, csaw 1.36+, PyDESeq2 0.5+, 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

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

Decision Tree: Choosing Normalization

Spike-In Scaling Factor Calculation

ChIP-Rx (Orlando 2014) uses Drosophila chromatin spike-in added at fixed concentration BEFORE IP. Egan 2016 adds a fixed MASS of Drosophila S2 chromatin (matched to target chromatin by the ~27:1 human:fly genome-size ratio), not a fixed cell count.

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=...):

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 (spikein = TRUE forces library = DBA_LIBSIZE_BACKGROUND internally)
dba_obj <- dba.normalize(dba_obj, spikein = TRUE,
                          normalize = DBA_NORM_LIB)

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

Internal-control sanity check: 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 saturated (always 100k+ reads); check titration linearity

Per the Patel et al 2024 Nat Biotechnol survey, improper spike-in normalization is common: of 53 datasets examined, only 27 (~51%) had adequate matched input controls across conditions.

DiffBind Workflow (BAMs + Peaks)

Goal: Run DiffBind from a sample sheet to consensus peaks, normalized counts, and tested differential binding.

Approach: Build sample sheet, count reads in consensus peaks (summit-centered for narrow marks, full peak width for broad), choose normalization based on the three-problem framework above, then test with DESeq2 or edgeR backend.

library(DiffBind)

# Sample sheet: SampleID, Condition, Replicate, bamReads, bamControl, Peaks, PeakCaller
dba_obj <- dba(sampleSheet = 'samples.csv')

# Counting: summits=250 (narrow); FALSE for broad histones (use full peak width)
dba_obj <- dba.count(dba_obj, summits = 250, minOverlap = 2, bParallel = TRUE)

# Normalization — choose per the three-problem framework
dba_obj <- dba.normalize(dba_obj,                      # default: full library size
                          method = DBA_DESEQ2,
                          normalize = DBA_NORM_LIB,
                          library = DBA_LIBSIZE_FULL)

# Background bin TMM for composition bias / broad marks
# dba_obj <- dba.normalize(dba_obj, background = TRUE)

# Spike-in scaling for global shifts
# dba_obj <- dba.normalize(dba_obj, spikein = TRUE)

dba_obj <- dba.contrast(dba_obj, design = '~ Condition')
dba_obj <- dba.analyze(dba_obj, method = DBA_DESEQ2)

# Always inspect what was actually applied
dba.normalize(dba_obj, bRetrieve = TRUE)

DiffBind 3.20+ defaults (verified via bRetrieve):

• summits = 200 — narrow recentering (set FALSE for broad histones)
• library = DBA_LIBSIZE_FULL — full library, conservative
• normalize = DBA_NORM_LIB — library size only, not reads-in-peaks RLE
• Blacklist filtering on
• background = FALSE — set TRUE for composition bias

csaw Workflow (Windows-Based)

Goal: Detect differential binding from sliding windows without committing to predefined peaks; robust to composition bias via background bin TMM.

Approach: Count reads in overlapping windows, normalize on 10 kb bins (composition bias) and/or loess (trended bias), test with edgeR QL-F, then merge significant windows into regions.

library(csaw)
library(edgeR)

bam_files <- c('ctrl_1.bam', 'ctrl_2.bam', 'treat_1.bam', 'treat_2.bam')
condition <- factor(c('ctrl', 'ctrl', 'treat', 'treat'))

# Window counts: 150 bp window, 50 bp spacing (sharp marks); 1-2 kb for broad
param <- readParam(minq = 30, pe = 'both', dedup = TRUE,
                    discard = import('hg38-blacklist.v2.bed'))
windows <- windowCounts(bam_files, width = 150, ext = 200, param = param)

# Composition bias via 10 kb bins (always applied for ChIP-seq)
bg_bins <- windowCounts(bam_files, bin = TRUE, width = 10000, param = param)
windows <- normFactors(bg_bins, se.out = windows)

# Optional: trended bias via non-linear loess (use cautiously)
# windows <- normOffsets(windows, se.out = TRUE)

# edgeR QL-F test
y <- asDGEList(windows)
design <- model.matrix(~condition)
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design, robust = TRUE)
results <- glmQLFTest(fit, coef = 2)

# Merge significant windows within 1 kb into regions
merged <- mergeResults(windows, results$table, tol = 1000, merge.args = list(max.width = 5000))

For broad marks, increase window width to 1-2 kb and merge tolerance to 5 kb. For TFs, 150 bp window + 50 bp spacing is standard.

DESeq2 from Count Matrix

library(DESeq2)

counts <- read.delim('counts.tsv', row.names = 1, check.names = FALSE)
coldata <- data.frame(
    condition = factor(c('ctrl', 'ctrl', 'ctrl', 'treat', 'treat', 'treat')),
    row.names = colnames(counts)
)

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~ condition)

# Pre-filtering for ChIP-seq is LESS aggressive than RNA-seq (peaks already enriched)
keep <- rowSums(counts(dds)) >= 1  # remove only all-zero
dds <- dds[keep, ]
dds$condition <- relevel(dds$condition, ref = 'ctrl')

dds <- DESeq(dds)
res <- results(dds, alpha = 0.05)  # match independent-filtering optimization

# Optional: LFC shrinkage for ranking/visualization only — does NOT change padj
library(apeglm)
resLFC <- lfcShrink(dds, coef = 'condition_treat_vs_ctrl', type = 'apeglm')

For spike-in normalization, set sizeFactors(dds) from scaling factors before DESeq(dds). For background-bin TMM (composition bias) use csaw to compute size factors, then transfer.

Per-Tool Failure Modes

DiffBind -- summits=200 default destroys broad marks

Trigger: Default DiffBind 3.x call with histone broad marks (H3K27me3, H3K9me3).

Mechanism: summits=200 re-centers peaks to summit ± 200 bp, throwing away most of a 10-100 kb broad domain.

Symptom: Differential count for broad marks much lower than expected; signal concentrated at narrow centers of broad regions.

Fix: dba.count(obj, summits = FALSE, ...) for broad marks; OR use full-width consensus peaks; OR switch to csaw with 1-2 kb windows.