Bio Applied Differential Binding

Differential Binding & Peak Annotation

When to Use

• Comparing ChIP-seq or ATAC-seq peak signal between two or more conditions (treatment vs. control, knockdown vs. WT) with replicates
• Building a consensus peak set from multiple MACS2 narrowPeak/broadPeak calls before quantitative comparison
• Running DESeq2- or edgeR-based differential testing on peak read counts (same statistical machinery as RNA-seq, applied to genomic intervals)
• Annotating differential (or all) peaks to genomic features — promoter, exon, intron, distal intergenic — and to nearest gene/TSS
• Generating volcano/MA plots and PCA of binding affinity to QC replicate consistency before biological interpretation

Version Compatibility

• R ≥ 4.3, Bioconductor ≥ 3.18
• DiffBind ≥ 3.12 (wraps DESeq2 ≥ 1.42 or edgeR ≥ 4.0)
• ChIPseeker ≥ 1.38, with a matching TxDb.* package (e.g. TxDb.Hsapiens.UCSC.hg38.knownGene) and org.Hs.eg.db
• Python ≥ 3.10, pandas ≥ 2.0, matplotlib ≥ 3.8 (for plotting exported results only — DiffBind/ChIPseeker have no Python equivalent)

Prerequisites

• R packages: BiocManager::install(c("DiffBind", "ChIPseeker", "clusterProfiler", "org.Hs.eg.db", "TxDb.Hsapiens.UCSC.hg38.knownGene"))
• Deduplicated, indexed BAMs per sample plus matching Input/IgG control BAMs
• Peak calls per sample (MACS2 narrowPeak/broadPeak) — see bio-chip-seq-peak-calling
• Familiarity with DESeq2-style count-based differential testing (bio-differential-expression-deseq2-basics)

Differential Binding with DiffBind

Goal: identify peaks whose read counts differ significantly between two conditions, correcting for library size and using replicate information.

Approach: build a sample sheet, count reads across a consensus peak set (peaks reproducible in ≥2 samples), normalize, then run DESeq2 (or edgeR) via dba.analyze().

Sample sheet (samplesheet.csv) — one row per sample:

SampleID,Condition,Replicate,bamReads,ControlID,bamControl,Peaks,PeakCaller
CTCF_Ctrl1,Control,1,dedup/ctcf_ctrl1.bam,Input1,dedup/input1.bam,peaks/ctcf_ctrl1_peaks.narrowPeak,narrow
CTCF_Ctrl2,Control,2,dedup/ctcf_ctrl2.bam,Input2,dedup/input2.bam,peaks/ctcf_ctrl2_peaks.narrowPeak,narrow
CTCF_Trt1,Treatment,1,dedup/ctcf_trt1.bam,Input3,dedup/input3.bam,peaks/ctcf_trt1_peaks.narrowPeak,narrow
CTCF_Trt2,Treatment,2,dedup/ctcf_trt2.bam,Input4,dedup/input4.bam,peaks/ctcf_trt2_peaks.narrowPeak,narrow

library(DiffBind)

# 1. Load sample sheet -> DBA object; peaks in >=2 samples form the consensus set
dba_obj <- dba(sampleSheet = "samplesheet.csv")
print(dba_obj)

# 2. Count reads in the consensus peak set (summarizeOverlaps backend is exact but slower)
dba_obj <- dba.count(dba_obj, bUseSummarizeOverlaps = TRUE, minOverlap = 2)

# 3. Normalize for library size / composition (RLE matches DESeq2's default)
dba_obj <- dba.normalize(dba_obj, normalize = DBA_NORM_RLE)

# 4. Define the contrast: Treatment vs. Control, requiring >=2 samples per group
dba_obj <- dba.contrast(dba_obj, categories = DBA_CONDITION, minMembers = 2)

# 5. Run the differential test (DESeq2 by default; DBA_EDGER for edgeR/TMM)
dba_obj <- dba.analyze(dba_obj, method = DBA_DESEQ2)

# 6. Extract significant peaks: FDR < 0.05 and |log2FC| > log2(1.5)
db_peaks <- dba.report(dba_obj, th = 0.05, fold = log2(1.5))
print(db_peaks)

# 7. QC plots: do replicates cluster by condition? is signal driven by a few peaks?
dba.plotPCA(dba_obj, DBA_CONDITION, label = DBA_ID)
dba.plotVolcano(dba_obj)

# Export for downstream annotation / sharing
rtracklayer::export(db_peaks, "diffbind_results.bed")
write.csv(as.data.frame(db_peaks), "diffbind_results.csv", row.names = FALSE)

Peak Annotation with ChIPseeker

Goal: map each (differential) peak to the nearest gene/TSS and a genomic feature category, for downstream GO/KEGG enrichment.

Approach: annotatePeak() against a TxDb for the genome build; tssRegion sets the promoter window; annoDb adds gene symbols/Entrez IDs.

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

peaks <- readPeakFile("diffbind_results.bed", as = "GRanges")

anno <- annotatePeak(
  peaks,
  tssRegion = c(-2000, 200),   # 2kb upstream to 200bp downstream of TSS
  TxDb      = txdb,
  annoDb    = "org.Hs.eg.db"
)

plotAnnoPie(anno)             # proportion of peaks per genomic feature
plotAnnoBar(anno)
plotDistToTSS(anno, title = "Distribution of peaks relative to TSS")

anno_df <- as.data.frame(anno)
head(anno_df[, c("seqnames", "start", "end", "annotation", "SYMBOL", "distanceToTSS")])

promoter_peaks <- anno_df[abs(anno_df$distanceToTSS) < 2000, ]
cat("Promoter-associated peaks:", nrow(promoter_peaks), "\n")

write.csv(anno_df, "diffbind_annotated.csv", row.names = FALSE)

Plotting DiffBind Results in Python

Goal: rebuild volcano and MA plots from an exported diffbind_results.csv without needing R installed on the plotting machine.

Approach: read the DESeq2-style columns (log2FoldChange/log2FC, padj/FDR) and classify by direction before plotting.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt


def classify_binding(df, fdr_col="FDR", lfc_col="log2FC", fdr_thresh=0.05, lfc_thresh=1.5):
    """Label each peak as Gained/Lost/Unchanged from DiffBind dba.report() output.

    df must have numeric fdr_col and lfc_col columns (e.g. loaded from
    diffbind_results.csv). Returns df with 'Significant' and 'Direction' added.
    """
    df = df.copy()
    df["Significant"] = (df[fdr_col] < fdr_thresh) & (df[lfc_col].abs() > np.log2(lfc_thresh))
    df["Direction"] = "Unchanged"
    df.loc[df["Significant"] & (df[lfc_col] > 0), "Direction"] = "Gained"
    df.loc[df["Significant"] & (df[lfc_col] < 0), "Direction"] = "Lost"
    return df


def plot_volcano(df, lfc_col="log2FC", fdr_col="FDR", fdr_thresh=0.05, lfc_thresh=1.5, ax=None):
    """Volcano plot of differential binding results, colored by direction."""
    ax = ax or plt.gca()
    colors = {"Gained": "coral", "Lost": "steelblue", "Unchanged": "lightgray"}
    neg_log10_fdr = -np.log10(df[fdr_col].clip(lower=1e-50))
    for direction, grp in df.groupby("Direction"):
        idx = grp.index
        ax.scatter(grp[lfc_col], neg_log10_fdr.loc[idx], c=colors[direction],
                   s=6 if direction != "Unchanged" else 3,
                   alpha=0.8 if direction != "Unchanged" else 0.4,
                   label=f"{direction} (n={len(grp):,})")
    ax.axhline(-np.log10(fdr_thresh), color="black", ls="--", lw=1, label=f"FDR={fdr_thresh}")
    ax.axvline(np.log2(lfc_thresh), color="gray", ls=":", lw=1)
    ax.axvline(-np.log2(lfc_thresh), color="gray", ls=":", lw=1)
    ax.set_xlabel("log2 Fold Change")
    ax.set_ylabel("-log10(FDR)")
    ax.legend(markerscale=2)
    return ax


# Example: results = pd.read_csv("diffbind_results.csv")  # from dba.report() export
# results = classify_binding(results)
# plot_volcano(results)
# plt.show()

Pitfalls

• Coordinate systems: BED is 0-based half-open; VCF/GFF/GRanges are 1-based inclusive — mixing them causes off-by-one errors when exporting DiffBind results and re-importing to ChIPseeker
• Consensus peak definition: minOverlap/minMembers defaults to peaks in ≥2 samples; too permissive inflates the peak set with singleton noise, too strict drops condition-specific real sites
• Normalization choice: for TF ChIP-seq use RLE/TMM (library-size driven); for broad marks with global signal shifts (e.g. total H3K27ac loss), consider spike-in or background-region normalization instead — default normalization can mask true global change
• Multiple testing: thousands of peaks are tested simultaneously — always use th (FDR/BH), never raw p-values, to call significance
• Batch effects: check dba.plotPCA() for replicate clustering by condition (not by batch/day) before trusting dba.report() output
• Promoter window choice: ChIPseeker's tssRegion default (-3000, 3000) is looser than the common (-2000, 200); mismatched windows between analyses make "% promoter peaks" non-comparable across papers

See Also

• bio-chip-seq-peak-calling — generating the narrowPeak/broadPeak inputs DiffBind requires
• bio-chip-seq-peak-annotation — deeper ChIPseeker/annotation workflows
• bio-differential-expression-deseq2-basics — the DESeq2 statistics DiffBind reuses under the hood
• bio-pathway-analysis-go-enrichment — clusterProfiler enrichment on annotated peak gene lists