macs3-peak-calling

MACS3 — ChIP-seq および ATAC-seq ピークコーラー

概要

MACS3（Model-based Analysis of ChIP-seq）は、ChIP-seq、ATAC-seq、CUT&RUN、CUT&TAG 実験から、リードが有意に濃縮された領域（ピーク）を特定します。ペアエンドデータからフラグメント長の分布をモデル化するか、シングルエンドデータでモノヌクレオソームのリードシフトからフラグメント長を推定し、ポアソンモデルを適用して、入力／IgG コントロールに対する濃縮倍率を特定します。MACS3 は、IGV や UCSC Genome Browser での可視化に使用できる、シグナルおよび q 値のトラック付きの BED 形式の narrowPeak（転写因子用）または broadPeak（ヒストンマーク用）ファイルを生成します。

使用する場面

• ChIP-seq 実験から転写因子結合ピークをコールする場合（--nomodel --extsize 200 を使用するか、MACS3 にフラグメント長を推定させます）
• ATAC-seq 実験からオープンクロマチン領域を特定する場合（--nomodel --shift -100 --extsize 200 -f BAMPE を使用します）
• --broad を使って幅広いヒストン修飾ピーク（H3K27me3、H3K9me3、H3K36me3）をコールする場合
• -B --SPMR を使って、ゲノムブラウザーで可視化するためのピークシグナルトラック（bedGraph/bigWig）を生成する場合
• DiffBind または DESeq2 への入力として MACS3 のピークを使用し、差次的結合解析を行う場合
• ヌクレオソーム解像度で ATAC-seq のピークをコールする場合は、HMMRATAC（MACS3 の一部）を使用します
• 代替手段として SPP または HOMER を使用します。MACS3 は ENCODE 推奨の標準ツールです

前提条件

• Python パッケージ：macs3（Python ≥ 3.8）
• 入力：ChIP-seq または ATAC-seq のアラインメントから得た、ソート済み BAM ファイル（インデックス付き）（例：STAR または Bowtie2 を使用）
• 任意：バックグラウンド正規化用の入力／IgG コントロール BAM

インストール前に確認してください：このツールは、現在の環境（例：pixi／conda 環境内）ですでに利用可能な場合があります。まず command -v macs3 を実行し、パスが返された場合は以下のインストールコマンドを省略してください。pixi プロジェクト内で実行する場合は、macs3 を直接実行するのではなく、pixi run macs3 を使用してください。

# Install with pip or conda
pip install macs3
# or
conda install -c bioconda macs3

# Verify
macs3 --version
# macs3 3.0.2

実行前の確認

解析コードを書く前に、ユーザーと以下の点を確認してください。

decisions:
  - id: D1
    param: assayType
    kind: required
    source: user
    ask: "What produced these reads - a transcription-factor ChIP, a histone mark, or open-chromatin (ATAC)?"
    default: null

  - id: D2
    param: controlSample
    kind: required
    source: data
    ask: "Is there a matched input or IgG control for this sample?"
    default: "none - background is estimated locally instead"

  - id: D3
    param: effectiveGenomeSize
    kind: derived
    source: upstream
    ask: "Which genome were the reads aligned to?"
    default: "carried from the alignment reference"

  - id: D4
    param: peakShape
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Are the enriched regions sharp summits or broad domains?"
    default: "narrow for TF and ATAC; broad for spreading histone marks"

  - id: D5
    param: significanceThreshold
    kind: required
    source: user
    ask: "How strong must enrichment be before a region is called a peak?"
    default: "q-value 0.05"

  - id: D6
    param: fragmentModel
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Should fragment length be modelled from the data, or fixed because the assay's signal is not paired summits?"
    default: "model for ChIP; fixed shift/extension for ATAC"

  - id: D7
    param: duplicateHandling
    kind: optional
    source: user
    depends_on: [D1]
    ask: "Are identical read positions PCR duplicates, or real independent insertions?"
    default: "one per position for ChIP; keep all for ATAC"

  - id: D8
    param: signalTracks
    kind: optional
    source: user
    ask: "Emit normalized coverage tracks alongside the peak calls?"
    default: "peaks only"

D1 を最初に確認するのは、後続の4つの設定が D1 に基づいて決まるためです。特に ATAC では、ChIP のデフォルト設定3つ（フラグメントモデルなし、負のシフト、重複リードを保持）を一度に反転させます。そのため、ATAC を ChIP のデフォルト設定で処理すると、一見もっともらしく見えるものの、誤ったピークが生成されます。

クイックスタート

# Call peaks for TF ChIP-seq (narrow peaks, with input control)
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample_tf \
    --outdir peaks/ \
    -q 0.05

# Output: peaks/sample_tf_peaks.narrowPeak
wc -l peaks/sample_tf_peaks.narrowPeak

ワークフロー

ステップ 1：入力 BAM ファイルを準備する

MACS3 には、ゲノムアラインメントから得たソート済み・インデックス付き BAM ファイルが必要です。

# Sort and index ChIP and control BAMs (if not already done)
samtools sort -@ 8 chip_raw.bam -o chip.bam
samtools sort -@ 8 input_raw.bam -o input.bam
samtools index chip.bam
samtools index input.bam

# Check read counts
echo "ChIP reads: $(samtools view -c -F 4 chip.bam)"
echo "Input reads: $(samtools view -c -F 4 input.bam)"

ステップ 2：ナローピークをコールする（TF ChIP-seq）

転写因子結合部位の特定には、デフォルトモードを使用します。

# TF ChIP-seq with input control
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n tf_chip \
    --outdir peaks/ \
    -q 0.05 \
    --keep-dup auto

echo "Peaks called: $(wc -l < peaks/tf_chip_peaks.narrowPeak)"
echo "Summit file: peaks/tf_chip_summits.bed"

# Without input control (less recommended)
macs3 callpeak \
    -t chip.bam \
    -f BAM \
    -g hs \
    -n tf_noinput \
    --outdir peaks/ \
    --nolambda

ステップ 3：ブロードピークをコールする（ヒストンマーク）

H3K27me3 や H3K36me3 のような広がりのあるヒストン修飾には、--broad を使用します。

# H3K27me3 broad histone mark
macs3 callpeak \
    -t h3k27me3.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n h3k27me3 \
    --outdir peaks/ \
    --broad \
    --broad-cutoff 0.1 \
    -q 0.05

echo "Broad peaks: $(wc -l < peaks/h3k27me3_peaks.broadPeak)"

# H3K4me3 (sharp mark — use narrow peaks)
macs3 callpeak \
    -t h3k4me3.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n h3k4me3 \
    --outdir peaks/ \
    -q 0.05

ステップ 4：ATAC-seq ピークをコールする

ATAC-seq では、Tn5 挿入部位に対する特別な処理が必要です。

# ATAC-seq with paired-end BAM (recommended)
macs3 callpeak \
    -t atac.bam \
    -f BAMPE \
    -g hs \
    -n atac_sample \
    --outdir peaks/ \
    --nomodel \
    --nolambda \
    -q 0.05 \
    --keep-dup all

echo "ATAC peaks: $(wc -l < peaks/atac_sample_peaks.narrowPeak)"

# Single-end ATAC-seq: shift reads to center on Tn5 cut site
macs3 callpeak \
    -t atac_se.bam \
    -f BAM \
    -g hs \
    -n atac_se \
    --outdir peaks/ \
    --nomodel \
    --shift -100 \
    --extsize 200 \
    --keep-dup all

ステップ 5：可視化用のシグナルトラックを生成する

ゲノムブラウザーで可視化するため、bedGraph および bigWig ファイルを生成します。

# Generate bedGraph normalized to million reads (SPMR)
macs3 callpeak \
    -t chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n chip_track \
    --outdir tracks/ \
    -B \
    --SPMR \
    --keep-dup auto

# Convert bedGraph to bigWig for IGV/UCSC
# Requires bedGraphToBigWig and chrom.sizes
sort -k1,1 -k2,2n tracks/chip_track_treat_pileup.bdg > tracks/chip_sorted.bdg
bedGraphToBigWig tracks/chip_sorted.bdg genome/hg38.chrom.sizes tracks/chip.bw

echo "BigWig track: tracks/chip.bw"

ステップ 6：ピークに注釈を付けて解析する

narrowPeak 出力を解析し、ピークにゲノム上の特徴を注釈付けします。

import pandas as pd

# Load narrowPeak file
# Columns: chrom, start, end, name, score, strand, signalValue, pValue, qValue, peak
cols = ["chrom", "start", "end", "name", "score", "strand",
        "signalValue", "pValue", "qValue", "peak"]
peaks = pd.read_csv("peaks/tf_chip_peaks.narrowPeak", sep="\t",
                    header=None, names=cols)

print(f"Total peaks: {len(peaks)}")
print(f"Peaks on chr1: {(peaks['chrom'] == 'chr1').sum()}")
print(f"Median peak width: {(peaks['end'] - peaks['start']).median():.0f} bp")
print(f"Peaks with q-value < 0.01: {(peaks['qValue'] > 2).sum()}")  # -log10(q) > 2

# Filter high-confidence peaks
high_conf = peaks[peaks["qValue"] > 2].copy()  # q < 0.01
high_conf["width"] = high_conf["end"] - high_conf["start"]
print(f"\nHigh-confidence peaks: {len(high_conf)}")
high_conf.to_csv("high_confidence_peaks.bed", sep="\t", index=False, header=False,
                 columns=["chrom", "start", "end", "name", "score", "strand"])

主要なパラメーター

パラメーター	デフォルト	範囲／選択肢	効果
-t / --treatment	必須	BAM/BED/SAM	ChIP または ATAC の処理対象ファイル
-c / --control	—	BAM/BED/SAM	Input／IgG コントロール。指定しない場合は --nolambda を省略します
-g / --gsize	必須	hs、mm、ce、dm、または整数	有効ゲノムサイズ。hs=2.7e9（ヒト）、mm=1.87e9（マウス）
-q / --qvalue	0.05	0–1	ピーク検出の FDR 閾値
-p / --pvalue	—	0–1	P 値のカットオフ（厳密に制御する場合は q 値の代わりに使用）
--broad	オフ	フラグ	拡散性ヒストンマークのブロードピークを検出します
--broad-cutoff	0.1	0–1	ブロード領域の統合に使う Q 値のカットオフ
--nomodel	オフ	フラグ	フラグメント長のモデリングをスキップします。ATAC-seq では必須です
--extsize	200	50–1000	--nomodel 設定時のフラグメント拡張サイズ
--shift	0	-500–500	リードのシフト量（bp）。ATAC-seq では --extsize 200 とともに -100 を使用します
--keep-dup	1	auto、all、整数	重複の処理。auto はポアソンモデルを使用し、all はすべて保持します（ATAC-seq）
-B / --bdg	オフ	フラグ	bedGraph シグナルトラックを書き出します
--SPMR	オフ	フラグ	bedGraph を 100 万リードあたりのシグナルに正規化します

よく使うレシピ

レシピ 1：複数サンプルの一括ピーク検出

#!/bin/bash
# Call peaks for multiple ChIP-seq samples with the same input
INPUT="input.bam"
GENOME="hs"
OUTDIR="peaks"
mkdir -p "$OUTDIR"

SAMPLES=(H3K4me3 H3K27ac H3K27me3 CTCF)
MODES=(narrow narrow broad narrow)

for i in "${!SAMPLES[@]}"; do
    sample="${SAMPLES[$i]}"
    mode="${MODES[$i]}"
    echo "Calling peaks: $sample ($mode)"
    
    if [ "$mode" == "broad" ]; then
        BROAD_FLAG="--broad --broad-cutoff 0.1"
    else
        BROAD_FLAG=""
    fi
    
    macs3 callpeak \
        -t "${sample}.bam" \
        -c "$INPUT" \
        -f BAM \
        -g "$GENOME" \
        -n "$sample" \
        --outdir "$OUTDIR" \
        $BROAD_FLAG \
        -q 0.05 \
        --keep-dup auto
    
    echo "$sample: $(wc -l < $OUTDIR/${sample}_peaks.*Peak) peaks"
done

レシピ 2：IDR（再現不能な検出率）を用いた再現性のあるピーク

# Call peaks on individual replicates (lenient thresholds for IDR)
for rep in rep1 rep2; do
    macs3 callpeak \
        -t "chip_${rep}.bam" \
        -c input.bam \
        -f BAM \
        -g hs \
        -n "tf_${rep}" \
        --outdir peaks/ \
        -p 0.1 \
        --keep-dup auto
done

# Run IDR to find reproducible peaks
# pip install idr
idr --samples peaks/tf_rep1_peaks.narrowPeak peaks/tf_rep2_peaks.narrowPeak \
    --input-file-type narrowPeak \
    --output-file peaks/tf_idr_peaks.txt \
    --idr-threshold 0.05 \
    --plot

echo "IDR peaks: $(wc -l < peaks/tf_idr_peaks.txt)"

想定される出力

出力	形式	説明
*_peaks.narrowPeak	BED6+4	シグナル、p 値、q 値、ピーク頂点のオフセットを含むナローピーク
*_peaks.broadPeak	BED6+3	ブロードピーク（--broad 指定時）：chrom、start、end、signal、p-val、q-val
*_summits.bed	BED3+2	スコア付きのピーク頂点位置（1 bp）。モチーフ解析に使用します
*_treat_pileup.bdg	bedGraph	処理対象のシグナルトラック（-B 指定時）
*_control_lambda.bdg	bedGraph	コントロール／ローカル lambda トラック（-B 指定時）
*_model.r	R スクリプト	フラグメントサイズモデル。Rscript *_model.r を実行するとプロットを作成します

トラブルシューティング

問題	原因	解決策
検出されるピークが非常に少ない	q 値が厳しすぎる、またはリード深度が低い	-p 1e-3 に緩和します。シーケンス深度を確認します（アライメント済みリード 10M 以上を推奨）
ピークが多すぎる（>100k）	閾値が緩すぎる、または Input コントロールがない	--control input.bam を追加します。-q 0.01 を使用し、signalValue でフィルタリングします
「no reads」でピーク検出が失敗する	BAM ファイルがソートまたはインデックス化されていない	MACS3 の前に samtools sort と samtools index を実行します
ATAC-seq のピークがミトコンドリアに検出される	mtDNA 含有量が高い	次のようにフィルタリングします：`samtools view -h chip.bam
フラグメントモデルの作成に失敗する	リードが少なすぎる、またはリード長が標準的でない	モデリングをスキップするため、--nomodel --extsize 200 を追加します
bedGraph 出力が非常に大きい	正規化されていない高カバレッジデータ	100 万リードあたりのシグナルに正規化するため、--SPMR を追加します
--broad でナローピークが検出されない	シグナルが実際には鋭い	ChIP の標的を確認します。TF と H3K4me3 にはナローモードが必要です
gsize が一致しない	アセンブリに対して誤ったゲノムサイズを使用している	hg19／hg38 には hs、mm9／mm10 には mm を使用するか、正確な整数を指定します

参考資料

• MACS3 GitHub: macs3-project/MACS — ソースコード、ドキュメント、変更履歴
• Zhang Y et al. (2008) “Model-based Analysis of ChIP-Seq (MACS)” — Genome Biology 9:R137. DOI:10.1186/gb-2008-9-9-r137
• ENCODE ATAC-seq パイプライン — MACS3 を使用する ENCODE 標準の ATAC-seq ワークフロー
• IDR フレームワーク — 再現性のあるピーク検出のための再現不能な検出率