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}_used_parameters.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:

• mcp__project-init-tools__project_init

with:

• sample

  : all

• task

  : peak_calling

• genome

  : provided by user

The tool will:

• Create

  ${sample}_peak_calling

  directory.
• Return the full path of the

  ${sample}_peak_calling

  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

--call-summits -q 0.01

--call-summits -q 0.05

--broad --broad-cutoff 0.1 -q 0.05

--call-summits -q 0.05

--broad

--nomodel --shift -100 --extsize 200 -q 0.05

Step 4. Execute MACS3 with Auto Parameters

Call:

• mcp__macs2-tools__run_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).

• out_dir

  : ${proj_dir}/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