bio-chipseq-peak-calling

## Version Compatibility

Reference examples tested with: MACS2 2.2+, MACS3 3.0+

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Peak Calling with MACS3

**"Call peaks from my ChIP-seq data"** → Identify significantly enriched regions (narrow peaks for TFs, broad peaks for histone marks) by comparing IP signal to input control.
- CLI: `macs3 callpeak -t chip.bam -c input.bam -f BAM -g hs -n sample`

MACS3 is the actively developed successor to MACS2. Commands are identical except the binary name. MACS2 is in maintenance mode.

## Basic Peak Calling

**Goal:** Call enriched regions from ChIP-seq alignments with input control normalization.

**Approach:** Compare treatment BAM signal against input control using MACS3 local Poisson model.

```bash
# Call peaks with input control (recommended)
macs3 callpeak -t chip.bam -c input.bam -f BAM -g hs -n sample --outdir peaks/

# For MACS2 (legacy), replace 'macs3' with 'macs2' - syntax is identical
```

## Without Input Control

**Goal:** Call peaks without a matched input/control sample.

**Approach:** Use MACS3 with genomic background estimation only (less accurate than with control).

```bash
# Not recommended, but possible
macs3 callpeak -t chip.bam -f BAM -g hs -n sample --outdir peaks/
```

## Narrow Peaks (TF, H3K4me3, H3K27ac)

**Goal:** Call sharp, well-defined peaks typical of transcription factors and active histone marks.

**Approach:** Use default narrow peak mode with q-value filtering and genome size correction.

```bash
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \                        # hs=human, mm=mouse, ce=worm, dm=fly
    -n sample_narrow \
    --outdir peaks/ \
    -q 0.05                        # q-value threshold
```

## Broad Peaks (H3K36me3, H3K27me3, H3K9me3)

**Goal:** Call diffuse, broad enrichment domains typical of repressive or elongation-associated histone marks.

**Approach:** Enable broad peak mode which links nearby enriched regions into broader domains.

```bash
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample_broad \
    --outdir peaks/ \
    --broad \                      # Broad peak mode
    --broad-cutoff 0.1             # Broad peak q-value
```

## Paired-End Data

**Goal:** Call peaks from paired-end sequencing using actual fragment sizes instead of modeled estimates.

**Approach:** Use BAMPE format so MACS3 calculates fragment size from mate pairs directly.

```bash
# MACS3 uses BAMPE format for paired-end
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAMPE \                     # Paired-end BAM
    -g hs \
    -n sample_pe \
    --outdir peaks/
```

## Multiple Replicates

**Goal:** Call peaks from multiple biological replicates pooled together for increased statistical power.

**Approach:** Provide all replicate BAMs to MACS3, which internally pools reads before peak calling.

```bash
# Pool replicates (MACS3 handles internally)
macs3 callpeak \
    -t rep1.bam rep2.bam rep3.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n pooled \
    --outdir peaks/
```

## Custom Genome Size

**Goal:** Call peaks for non-model organisms without a built-in genome size shortcut.

**Approach:** Provide the effective genome size as a numeric value instead of a species abbreviation.

```bash
# For non-model organisms or custom genomes
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g 2.7e9 \                     # Effective genome size in bp
    -n sample \
    --outdir peaks/
```

## Common Genome Sizes

| Genome | Flag | Effective Size |
|--------|------|----------------|
| Human | hs | 2.7e9 |
| Mouse | mm | 1.87e9 |
| C. elegans | ce | 9e7 |
| D. melanogaster | dm | 1.2e8 |

## Fixed Fragment Size

**Goal:** Call peaks when fragment size modeling fails or a specific extension size is needed.

**Approach:** Bypass model building and specify a fixed read extension size manually.

```bash
# If modeling fails or for ATAC-seq
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    --nomodel \                    # Skip model building
    --extsize 200 \                # Fixed extension size
    -n sample \
    --outdir peaks/
```

## Generate Signal Tracks

**Goal:** Produce normalized signal tracks for genome browser visualization alongside peak calls.

**Approach:** Enable bedGraph output with signal-per-million-reads normalization, then convert to bigWig.

```bash
# Generate bedGraph and bigWig files
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample \
    --outdir peaks/ \
    -B \                           # Generate bedGraph
    --SPMR                         # Signal per million reads

# Convert to bigWig (requires UCSC tools)
sort -k1,1 -k2,2n peaks/sample_treat_pileup.bdg > peaks/sample.sorted.bdg
bedGraphToBigWig peaks/sample.sorted.bdg chrom.sizes peaks/sample.bw
```

## Local Lambda for Broad Marks

**Goal:** Improve broad peak calling by disabling the genome-wide lambda estimate.

**Approach:** Use --nolambda to rely solely on local background estimation for very broad domains.

```bash
# Recommended for very broad marks
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    --broad \
    --nolambda \                   # Use local lambda only
    -n sample \
    --outdir peaks/
```

## Cutoff Analysis

**Goal:** Evaluate how different significance thresholds affect the number of called peaks.

**Approach:** Run MACS3 cutoff analysis mode to generate a table of peak counts at various q-value cutoffs.

```bash
# Test different q-value cutoffs
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    --cutoff-analysis \            # Generate cutoff analysis file
    -n sample \
    --outdir peaks/
```

## Output Files

| File | Description |
|------|-------------|
| *_peaks.narrowPeak | Peak coordinates (BED6+4) |
| *_peaks.broadPeak | Broad peak coordinates |
| *_summits.bed | Peak summit positions |
| *_model.r | R script for model visualization |
| *_treat_pileup.bdg | Treatment signal (with -B) |
| *_control_lambda.bdg | Control signal (with -B) |

## narrowPeak Format

```
chr1  100  200  peak_1  100  .  5.2  10.5  8.3  50
```
Columns: chr, start, end, name, score, strand, signalValue, pValue, qValue, peak

## Filter Peaks

**Goal:** Post-filter called peaks by statistical significance or signal strength.

**Approach:** Use awk on narrowPeak columns to apply q-value or signal-value cutoffs.

```bash
# Filter by q-value
awk '$9 > 2' peaks.narrowPeak > peaks.filtered.narrowPeak  # -log10(q) > 2 means q < 0.01

# Sort by signal strength
sort -k7,7nr peaks.narrowPeak > peaks.sorted.narrowPeak
```

## Key Parameters

| Parameter | Default | Description |
|-----------|---------|-------------|
| -t | required | Treatment BAM file(s) |
| -c | none | Control BAM file(s) |
| -f | AUTO | Format (BAM, BAMPE, BED) |
| -g | hs | Genome size |
| -n | NA | Output prefix |
| -q | 0.05 | Q-value cutoff |
| -p | none | P-value cutoff (overrides -q) |
| --broad | false | Broad peak calling |
| --nomodel | false | Skip model building |
| --extsize | 200 | Extension size (with --nomodel) |
| -B | false | Generate bedGraph |
| --SPMR | false | Signal per million reads |

## Related Skills

- peak-annotation - Annotate peaks to genes
- differential-binding - Compare peaks between conditions
- alignment-files - Prepare BAM files
- chipseq-visualization - Visualize peaks