single-cell-atac-seq-peak-calling-annotaion

ATAC-seq Peak Calling and Differential Accessibility

Call accessible chromatin peaks from ATAC-seq BAM files, annotate peaks to genomic features and genes, and identify differentially accessible regions between experimental conditions.

This is Step 2 of the bulk ATAC-seq pipeline.

---

What it does

1. Calls peaks with MACS2 (--nomodel --shift -100 --extsize 200) optimized for ATAC-seq
2. Filters peaks overlapping ENCODE blacklist regions
3. Builds a consensus peak set across all samples using bedtools merge
4. Counts reads per sample in consensus peaks (featureCounts or bedtools coverage)
5. Annotates peaks with genomic features (promoter, UTR, exon, intron, intergenic) and nearest gene
6. Runs differential accessibility analysis with DESeq2 on pseudobulk counts (recommended) or DiffBind
7. Applies FDR correction and generates ranked DAR table and volcano plot

---

Why this exists

If you ask a general AI to "call peaks from ATAC-seq data," it will:

• Use default MACS2 parameters (designed for ChIP-seq), not ATAC-seq-specific parameters
• Not build a consensus peak set across samples, making cross-sample comparison impossible
• Not remove ENCODE blacklist regions, leaving artifactual peaks from repetitive elements
• Use --format BAM instead of --format BAMPE (ATAC-seq is paired-end)
• Apply DESeq2 directly to individual reads per cell instead of pseudobulk aggregation for multi-sample data

This skill encodes the correct methodological decisions:

• Uses ATAC-seq-specific MACS2 flags: --nomodel --shift -100 --extsize 200 --format BAMPE
• Filters ENCODE blacklist regions (hg38 or mm10) that produce artifactual signal
• Builds a reproducible consensus peak set using IDR or bedtools merge across replicates
• Applies pseudobulk DESeq2 for multi-sample differential analysis (controls type I error)

---

Reference Methods

MACS2 ATAC-seq parameters:

• --nomodel: Skip the ChIP enrichment model building (ATAC-seq does not have a broad enrichment model)
• --shift -100 --extsize 200: Centers signal on Tn5 cut site by shifting reads 100 bp upstream and extending 200 bp
• --format BAMPE: Reads paired-end fragment coordinates directly from BAM
• --nolambda: Disables local lambda background estimation (optional; use with high-coverage data)