Peak Calling

Overview

This skill automatically performs core peak calling with MACS2 for ChIP-seq and ATAC-seq data, based on the BAM files in the current directory. It includes automatic experiment recognition and parameter selection.

Main steps include:

• Refer to the Inputs & Outputs section to check inputs and build the output architecture. All the output file should located in ${proj_dir} in Step 0.
• Always prompt user for genome_size to use (e.g. hs or mm). Never decide by yourself.
• Always prompt user if required control files are missing for ChIP-seq data.
• Always prompt user for the q value cutoff for peak calling.
• Detect experiment type (TF, histone mark, or ATAC-seq).
• Automatically decide whether to call narrow or broad peaks.
• Always use filtered BAM file (filtered.bam) if available.
• Detect sequencing type (single-end or paired-end) using SAM/BAM flags.
• Perform MACS3 peak calling accordingly.
• Generate a parameter log file (${sample}usedparameters.txt) with justification for each chosen option.

Inputs & Outputs

Inputs

${sample}.bam # filtered bam files

Outputs

all_peak_calling/
  peaks/
    ${sample}.narrowPeak # or ${sample}.broadPeak
  temp/
  logs/
    ${sample}_used_parameters.txt

Decision Tree

Step 0: Initialize Project

Call:

• mcpproject-init-toolsproject_init

with:

• sample: all
• task: peak_calling
• genome: provided by user

The tool will:

• Create ${sample}peakcalling directory.
• Return the full path of the ${sample}peakcalling directory, which will be used as ${proj_dir}.

Step 1. Identify and Classify BAM Files

Command Example

find . -name "*.bam" | sort

• Treatment BAMs: filenames contain TF names or histone marks (e.g., CTCF, H3K27me3, ATAC).
• Control BAMs: filenames contain “input”, “control”, or “IgG”.
• Prefer filtered BAMs (*.filtered.bam) if available.

Step 2. Detect Sequencing Type (Single-End or Paired-End)

Command Example

samtools flagstat sample.bam | egrep "properly paired|singletons"

• If properly paired > 0 → Paired-end (-f BAMPE)
• If singletons ≈ total → Single-end (-f BAM)

Step 3. Detect Experiment Type and Choose Peak Mode

Step 4. Execute MACS3 with Auto Parameters

Call:

• mcpmacs2-toolsrun_macs2

with:

• treatment_file: Path to treatment BAM file.
• control_file: Path to control/input BAM file. Required for ChIP-seq data. Prompt the user for the required file if not provided.
• genome_size: Always provided by user.
• name: Experiment name (prefix for output files).
• outdir: ${projdir}/peaks
• broad: If True, call broad peaks (for histone marks).
• broad_cutoff: Cutoff for broad region calling.
• qvalue: Q-value cutoff for peak detection. Prompt the user for the q value cutoff.
• format: use BAMPE for pair-end data, BAM for single-end data.
• nomodel: True for ATAC-seq, False for ChIP-seq.
• shift: Shift size in bp (e.g., -100 for ATAC-seq).
• extsize:"Extension size in bp (e.g., 200 for ATAC-seq).

Step 5. Generate Parameter Log File

After auto-selection, the skill writes a log file:

Example content:

Genome detected: 
Experiment type: H3K27me3 (broad histone)
Sequencing type: paired-end
Control used: input_control.bam
MACS3 mode: --broad --broad-cutoff 0.1 -q 0.05
Reasoning:
- Broad mark (H3K27me3) requires domain-level detection
- Control detected and applied
- Genome identified as <*>; using -g <*>
- Paired-end library; use -f BAMPE