Skip to main content

STEP 6: Peak Calling

Peak calling identifies regions of the genome where reads are significantly enriched, representing protein-DNA binding sites (ChIP-seq) or open chromatin regions (ATAC-seq).

MACS2 Peak Calling

MACS2 (Model-based Analysis of ChIP-Seq) is the standard tool for calling peaks. There is also lanceotron, a machine learning based peak caller, which can be used as an alternative however we will focus on MACS2 here.

Peak calling essentially determines where the signal (read coverage) is significantly higher than background noise, using a statistical model to identify these enriched regions.

source: https://link.springer.com/protocol/10.1007/978-1-0716-2481-4_13

Loading MACS2

On the CCB:

module load macs2

MACS2 Parameters Explained

  • -t - Treatment BAM file (your sample)
  • -c - Control BAM file (input DNA or IgG control)
  • -f - Format: BAM (single-end) or BAMPE (paired-end)
  • -g - Genome size: mm (mouse), hs (human)
  • -n - Name prefix for output files
  • --outdir - Output directory
  • -q - Q-value (FDR) cutoff (default 0.05)
  • --broad - Call broad peaks instead of narrow peaks
  • --nomodel - Skip model building (useful for ATAC-seq)
  • --shift - Shift reads (useful for ATAC-seq to center on nucleosomes)
  • --extsize - Extend reads to this size
  • --keep-dup - How to handle duplicate reads (auto or all)

Transcription Factor ChIP-seq (Narrow Peaks)

For transcription factors and proteins with specific binding sites:

macs2 callpeak -t treatment.bam -c control.bam \
  -f BAM -g mm -n sample_name \
  --outdir peaks/ \
  -q 0.01

Histone Modification ChIP-seq (Broad Peaks)

For histone marks that cover broader regions:

macs2 callpeak -t treatment.bam -c control.bam \
  -f BAM -g mm -n sample_name \
  --broad --broad-cutoff 0.1 \
  --outdir peaks/

ATAC-seq Peak Calling

ATAC-seq typically uses narrow peak calling with fragment size adjustment:

macs2 callpeak -t sample.bam \
  -f BAMPE -g mm -n sample_name \
  --nomodel --shift -100 --extsize 200 \
  --keep-dup all \
  --outdir peaks/ \
  -q 0.01

Output Files

MACS2 generates several important files:

  1. narrowPeak or broadPeak1 BED format with peak locations and scores:
chromosome  start  end  name  score  strand  fold_enrichment  -log10(p-value)  -log10(q-value)  summit
chr1  1000  2000  peak_1  100  .  5.2  10.5  8.3  500

  • summits.bed Precise peak summit locations (useful for motif analysis)

  • peaks.xls Tab-delimited file with detailed peak information

Quality Control

After peak calling, assess quality:

# Count number of peaks
wc -l sample_name_peaks.narrowPeak

# Check FRiP (Fraction of Reads in Peaks)
# Higher FRiP = better enrichment (>5% is good, >20% is excellent)

Filtering Peaks

Filter peaks by score or fold enrichment:

# Keep peaks with fold enrichment > 5
awk '$7 > 5' sample_name_peaks.narrowPeak > filtered_peaks.bed

# Keep peaks with q-value > 10 (-log10 scale)
awk '$9 > 10' sample_name_peaks.narrowPeak > high_confidence_peaks.bed

Visualization

Convert peaks to bigBed for genome browser viewing:

# Sort peaks
sort -k1,1 -k2,2n peaks.narrowPeak > peaks_sorted.bed

# Convert to bigBed
bedToBigBed peaks_sorted.bed chrom.sizes peaks.bb

Next Steps

After peak calling:

  1. Annotate peaks (assign to nearest genes)
  2. Motif analysis (find transcription factor binding motifs)
  3. Differential peak analysis between conditions
  4. Functional enrichment analysis

DiffBind and ChIPseeker are popular R packages for downstream peak analysis.