Identify chromatin regions with significantly different accessibility between conditions. Choose between consensus-peak workflows (DiffBind, DESeq2 directly) and peak-free sliding-window workflows (csaw). Pick the correct normalization (reads-in-peaks vs full-library vs spike-in) for the biological setting, correct batch effects with SVA/RUVseq, and report effect sizes with calibrated FDR.
BiocManager::install(c('DiffBind', 'DESeq2', 'edgeR', 'csaw', 'limma',
'ChIPseeker', 'sva', 'RUVSeq', 'GenomicRanges',
'TxDb.Hsapiens.UCSC.hg38.knownGene', 'org.Hs.eg.db'))
conda install -c bioconda subread bedtools homer
Inputs: per-sample BAM files (deduplicated, chrM-stripped), per-sample peak files (narrowPeak), and a sample sheet defining conditions and any covariates.
Tell your AI agent what you want to do: - "Run DiffBind with default DESeq2 backend on 3 control vs 3 treated ATAC samples" - "Use edgeR QL F-test because I only have 2 replicates per condition" - "Switch DiffBind normalization from reads-in-peaks to full-library because the treatment causes global compaction" - "Run csaw sliding-window analysis since differentiation collapses peak calls" - "Add SVA surrogate variables for hidden batch and re-fit" - "Filter to FDR < 0.05 AND abs(log2FC) >= 1, then annotate to nearest gene with ChIPseeker"
"Run DiffBind on three replicates per group with DESeq2 backend,
summits=250for fixed-width counting, native normalization (DBA_NORM_NATIVE, i.e. DESeq2 RLE), and report peaks with FDR < 0.05 and abs(log2FC) >= 1."
"I have 2 reps per condition. Use DiffBind with edgeR backend (better at low n) or skip DiffBind and run edgeR QL directly. Do not use DESeq2 + apeglm shrinkage because shrinkage at n=2 is over-aggressive."
"Treatment is an HDAC inhibitor that globally opens chromatin. Run DiffBind with
normalize=DBA_NORM_LIB(full library; also DiffBind's default) instead of reads-in-peaks (library=DBA_LIBSIZE_PEAKREADS), since RiP normalization will erase the global biology."
"Conditions are pre- and post-differentiation. Peak structure differs dramatically. Use csaw with
width=150windows,filterWindowsGlobalat 3-fold over background ($filter > log2(3)), edgeR QL F-test, then merge significant adjacent windows."
"Three batches of samples processed weeks apart. Run svaseq with n.sv=2 to estimate hidden surrogates, then re-fit DESeq2 with
~SV1 + SV2 + conditiondesign."
"Adjust for batch and donor as covariates. Use DiffBind
dba.contrast(..., design='~Batch + Donor + Condition')with explicit factor levels."
"After DA analysis, annotate peaks to genes using ChIPseeker with promoter region
-2000 to +500, then run GO enrichment on opened-vs-closed gene lists separately."
SampleID,Condition,Replicate,Batch,bamReads,Peaks,PeakCaller
ctrl_rep1,control,1,A,ctrl_rep1.bam,ctrl_rep1_peaks.narrowPeak,macs
ctrl_rep2,control,2,B,ctrl_rep2.bam,ctrl_rep2_peaks.narrowPeak,macs
treat_rep1,treated,1,A,treat_rep1.bam,treat_rep1_peaks.narrowPeak,macs
treat_rep2,treated,2,B,treat_rep2.bam,treat_rep2_peaks.narrowPeak,macs
Add Donor, Sex, Time, etc. columns as needed and include in the design formula.
| Replicate count | Method |
|---|---|
| 1 per group | Cannot do statistical DA; report log2FC only with caveat |
| 2 per group | edgeR QL (DiffBind edgeR backend); avoid DESeq2 + apeglm |
| 3 per group | DiffBind DESeq2 (default) is fine |
| 5+ per group | DESeq2 with apeglm shrinkage; can also use limma-voom |
| 10+ heterogeneous | DESeq2 LRT for primary effect + interactions |
| Treatment biology | Normalization |
|---|---|
| Localized differential (most experiments) | DiffBind DBA_NORM_NATIVE (RLE for DESeq2 / TMM for edgeR; DiffBind's own default is DBA_NORM_LIB) |
| Global accessibility shift (HDACi, DNMTi) | DBA_NORM_LIB OR exogenous spike-in |
| Outlier-prone counts | DBA_NORM_TMM (edgeR-style) |
| Cross-cell-type comparison | DBA_NORM_LIB (cell-type biology often global) |
| Goal | Threshold |
|---|---|
| Standard reporting | FDR < 0.05, abs(log2FC) >= 1 |
| Conservative high-confidence | FDR < 0.01, abs(log2FC) >= 1 |
| Primary cell types (smaller effects) | FDR < 0.05, abs(log2FC) >= 0.585 (1.5x) |
| Discovery for follow-up validation | FDR < 0.1, shrunken log2FC magnitude reported |
library=DBA_LIBSIZE_PEAKREADS) erases global accessibility shifts; keep the full-library default (DBA_NORM_LIB/DBA_LIBSIZE_FULL) whenever the biology is whole-genome.summits=250 to re-center peaks on summit; otherwise peak-width differences inflate counts.filterByExpr) before fitting; tail peaks inflate dispersion and FDR.tssRegion=c(-2000, 500) not the ChIPseeker default (-3000, 3000) -- the default over-counts promoter assignments by lumping in distal elements.