snakemake-workflow-engine

Snakemake — Python ワークフローエンジン

概要

Snakemake は、ノート PC から HPC やクラウドまで分析を拡張できる、Python ベースのワークフロー管理システムです。ワークフローは、入力ファイルと出力ファイルの依存関係を明示したルールとして定義します。Snakemake が実行順序を自動的に解決し、互いに独立したステップを並列に実行します。ルールからシェルコマンド、Python／R／Julia スクリプト、またはインライン Python を呼び出せます。ルールごとに conda 環境または Singularity 環境を設定できるため、ワークフローを完全に再現できます。バイオインフォマティクス分野で、NGS、ゲノムアセンブリ、バリアントコールのパイプラインに広く利用されています。

使用する場面

• 再現可能な複数ステップのバイオインフォマティクスパイプラインを構築する場合（アライン → ソート → バリアントコール → アノテーション）
• コードを変更せずに、ローカルでの開発から SLURM クラスターまで同じワークフローを拡張する場合
• ワイルドカードベースのルールを使用して、複数のサンプルを同じ方法で処理する場合
• 依存関係を自動管理する場合。入力が変更されたステップだけが再実行されます
• ツールを分離するために、ルールごとに conda 環境または Singularity 環境をデプロイする場合
• 計算リソースを確保する前に、DAG の図やドライランのプレビューを生成する場合
• Groovy DSL とデータフローチャネルが必要な場合、または nf-core コミュニティのパイプラインライブラリを活用する場合は、代わりに Nextflow を使用してください
• 単純なシェルループには bash スクリプトを使用してください。分岐を含む 3 つ以上の連続したステップがある場合にのみ、Snakemake を使う価値があります
• 動的なタスクグラフや時間ベースのスケジューリングを含むデータエンジニアリングのワークフローには、代わりに Prefect または Airflow を使用してください

前提条件

• Python パッケージ：snakemake、graphviz（DAG の可視化用）
• 環境：Python 3.11 以降。ルールごとの環境には conda／mamba を推奨
• データ要件：入力ファイル、参照ファイル。出力パスはルールとして定義します

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

# Install via conda (includes optional dependencies)
conda install -c conda-forge -c bioconda snakemake

# Minimal pip install
pip install snakemake

# Verify
snakemake --version
# 8.x.x

クイックスタート

# Snakefile — minimal 2-rule pipeline
SAMPLES = ["sampleA", "sampleB"]

rule all:             # Target rule: request final outputs
    input:
        expand("results/{sample}.sorted.bam", sample=SAMPLES)

rule align:
    input:
        fastq="data/{sample}.fastq",
        ref="refs/genome.fa"
    output:
        bam="results/{sample}.sorted.bam"
    threads: 4
    shell:
        "bwa mem -t {threads} {input.ref} {input.fastq} "
        "| samtools sort -@ {threads} -o {output.bam}"

# Run: dry-run first, then execute
snakemake -n            # dry-run: show what would run
snakemake --cores 8     # execute with 8 cores

コア API

モジュール 1：ルールの定義

各ルールでは、入力、出力、実行方法を指定して、1 つの分析ステップを定義します。

# Shell rule: run a command with {input} and {output} placeholders
rule fastqc:
    input:
        fastq="data/{sample}.fastq"
    output:
        html="qc/{sample}_fastqc.html",
        zip="qc/{sample}_fastqc.zip"
    log:
        "logs/fastqc/{sample}.log"
    shell:
        "fastqc {input.fastq} -o qc/ 2> {log}"

# Run rule: inline Python for logic-heavy steps
rule parse_stats:
    input:
        txt="results/{sample}.flagstat.txt"
    output:
        csv="results/{sample}.stats.csv"
    run:
        import re, csv
        lines = open(input.txt).readlines()
        mapped = re.search(r"(\d+) mapped", "".join(lines)).group(1)
        with open(output.csv, "w") as f:
            csv.writer(f).writerow([wildcards.sample, mapped])

# Script rule: delegate to external R/Python/Julia script
rule plot_coverage:
    input:
        depth="results/{sample}.depth.txt"
    output:
        pdf="results/{sample}.coverage.pdf"
    script:
        "scripts/plot_coverage.R"
    # In the R script, access via snakemake object:
    # depth_file <- snakemake@input[["depth"]]
    # pdf_path <- snakemake@output[["pdf"]]

モジュール 2：ワイルドカードとパターン展開

ワイルドカードを使うと、1 つのルールで任意の数のサンプルを処理できます。expand() は必要なファイルパスをすべて生成します。

# Define sample list (from config or glob)
SAMPLES = ["ctrl_rep1", "ctrl_rep2", "treat_rep1", "treat_rep2"]

rule all:
    input:
        # expand() generates: qc/ctrl_rep1_fastqc.html, qc/ctrl_rep2_fastqc.html, ...
        expand("qc/{sample}_fastqc.html", sample=SAMPLES),
        expand("results/{sample}.bam", sample=SAMPLES)

# Access wildcard values inside shell/run
rule align:
    input:
        "data/{sample}.fastq"
    output:
        "results/{sample}.bam"
    shell:
        "echo Processing {wildcards.sample}; "
        "bwa mem refs/genome.fa {input} | samtools view -b > {output}"

# Wildcard constraints prevent ambiguous matches
rule process:
    input:
        "data/{sample}_{rep}.fastq"
    output:
        "results/{sample}_{rep}.txt"
    wildcard_constraints:
        sample="[A-Za-z]+",   # letters only
        rep="\d+"             # digits only

# multiext: multiple outputs sharing a common path base
rule bwa_index:
    input:
        "refs/genome.fa"
    output:
        multiext("refs/genome.fa", ".amb", ".ann", ".bwt", ".pac", ".sa")
    shell:
        "bwa index {input}"

モジュール 3：設定とパラメーター

設定ファイルを使って設定を外部化し、params を使ってファイル依存関係を設けずにルールごとの値を渡します。

# Snakefile: declare config file
configfile: "config/config.yaml"

# config/config.yaml:
# samples: [ctrl, treat]
# threads:
#   align: 8
#   sort: 4
# min_mapq: 20

SAMPLES = config["samples"]

rule filter_reads:
    input:
        "results/{sample}.bam"
    output:
        "results/{sample}.filtered.bam"
    params:
        mapq=config["min_mapq"]    # from config, not a file
    threads:
        config["threads"]["sort"]
    shell:
        "samtools view -q {params.mapq} -b {input} > {output}"

# Dynamic params via lambda functions
rule trim:
    input:
        fastq="data/{sample}.fastq"
    output:
        trimmed="trimmed/{sample}.fastq"
    params:
        # Adapt quality threshold based on sample name
        quality=lambda wildcards: 25 if "ctrl" in wildcards.sample else 20
    shell:
        "fastp -q {params.quality} -i {input.fastq} -o {output.trimmed}"

モジュール 4：リソースと環境

スケジューラーとの連携のために計算リソースを宣言し、ツールの分離には conda／Singularity を使用します。

# Resource declaration (used by SLURM/LSF profiles)
rule variant_calling:
    input:
        bam="results/{sample}.deduped.bam",
        ref="refs/genome.fa"
    output:
        vcf="variants/{sample}.vcf.gz"
    resources:
        mem_mb=16000,       # memory in MB
        runtime=240,        # max walltime in minutes
        disk_mb=20000       # scratch disk space
    threads: 8
    shell:
        "bcftools mpileup -f {input.ref} {input.bam} "
        "| bcftools call -m -Oz -o {output.vcf}"

# Conda environment per rule (for reproducibility)
rule star_align:
    input:
        reads="data/{sample}.fastq",
        genome_dir="refs/star_index/"
    output:
        bam="star_out/{sample}/Aligned.sortedByCoord.out.bam"
    conda:
        "envs/star.yaml"
    # envs/star.yaml:
    # channels:
    #   - bioconda
    # dependencies:
    #   - star=2.7.10b
    #   - samtools=1.17
    threads: 8
    shell:
        "STAR --runThreadN {threads} --genomeDir {input.genome_dir} "
        "--readFilesIn {input.reads} --outSAMtype BAM SortedByCoordinate"

# Singularity/Apptainer container
rule gatk_haplotypecaller:
    input:
        bam="results/{sample}.bam",
        ref="refs/genome.fa"
    output:
        gvcf="gvcfs/{sample}.g.vcf.gz"
    container:
        "docker://broadinstitute/gatk:4.4.0.0"
    shell:
        "gatk HaplotypeCaller -I {input.bam} -R {input.ref} "
        "-O {output.gvcf} -ERC GVCF"

モジュール 5：実行とクラスター用プロファイル

ローカル、クラスター、またはクラウドで実行できます。プロファイルを使うと、Snakefile を変更せずにエグゼキューターを設定できます。

# Local execution
snakemake --cores 8                   # use 8 CPU cores
snakemake --cores all                 # use all available cores

# Dry run: show tasks without executing
snakemake -n --cores 8
# Output: 12 of 24 steps are complete. 12 jobs to run.

# Force rerun (ignore existing outputs)
snakemake --forceall --cores 8

# Visualize DAG as PDF
snakemake --dag | dot -Tpdf > workflow_dag.pdf

# SLURM cluster profile (profiles/slurm/config.yaml)
# executor: slurm
# jobs: 50
# default-resources:
#   mem_mb: 2000
#   runtime: 60
# use-conda: true

# Run with profile (cluster submit + monitor)
snakemake --profile profiles/slurm --cores 128

# Override resources at runtime
snakemake --profile profiles/slurm \
    --set-resources variant_calling:mem_mb=32000 --cores 128

# Override threads
snakemake --set-threads align=16 --cores 64

モジュール 6：特殊な出力タイプとユーティリティ

一時ファイル、保護された出力、チェックポイント、出力検証を扱います。

# temp: auto-delete after downstream rules consume it
rule sort_bam:
    input:
        "results/{sample}.raw.bam"
    output:
        temp("results/{sample}.sorted_temp.bam")  # deleted after indexing
    shell:
        "samtools sort {input} -o {output}"

# protected: write-protect final outputs (prevent overwrite)
rule final_report:
    input:
        "results/{sample}.vcf.gz"
    output:
        protected("reports/{sample}.final.vcf.gz")
    shell:
        "cp {input} {output}"

# directory: rule that outputs a directory
rule denovo_assembly:
    input:
        fastq="data/{sample}.fastq"
    output:
        directory("assemblies/{sample}/")
    shell:
        "spades.py -s {input.fastq} -o {output}"

# touch: create empty flag file (for ordering-only dependencies)
rule validate_bam:
    input:
        "results/{sample}.bam"
    output:
        touch("checkpoints/{sample}.validated")
    shell:
        "samtools quickcheck {input} && echo OK"

# ensure: validate output properties before considering rule complete
rule download_reference:
    output:
        ensure("refs/genome.fa", min_size=1_000_000)
    shell:
        "wget -O {output} https://example.com/genome.fa"

主要な概念

ルールの解決と DAG

Snakemake はターゲットから逆方向に処理します。必要な出力ファイルの一覧を指定すると、それらの生成に必要なルールの DAG を構築します。現在のターゲットに必要のないルールは無視されます。

# rule all: declare all final outputs here
# Without this, snakemake runs only the first rule
rule all:
    input:
        expand("results/{sample}.vcf.gz", sample=SAMPLES),
        expand("qc/{sample}_fastqc.html", sample=SAMPLES)

ワイルドカードと expand の違い

• ルールの入力／出力にある {sample} = ワイルドカード：実行時に Snakemake が値を設定します
• expand("results/{sample}.bam", sample=SAMPLES) = Python：文字列のリストを今すぐ生成します（rule all で使用）

よくあるワークフロー

ワークフロー 1：標準的な NGS QC パイプライン

目標：複数のサンプルに対して FastQC → トリミング → アライン → ソート → 重複除去 → flagstat を実行します。

configfile: "config/config.yaml"
SAMPLES = config["samples"]

rule all:
    input:
        expand("qc/{sample}_fastqc.html", sample=SAMPLES),
        expand("results/{sample}.flagstat.txt", sample=SAMPLES)

rule fastqc:
    input:  "data/{sample}.fastq"
    output: "qc/{sample}_fastqc.html", "qc/{sample}_fastqc.zip"
    shell:  "fastqc {input} -o qc/"

rule trim:
    input:  "data/{sample}.fastq"
    output: "trimmed/{sample}.fastq"
    shell:  "fastp -q 20 -i {input} -o {output}"

rule align:
    input:
        fastq="trimmed/{sample}.fastq",
        ref="refs/genome.fa"
    output: temp("results/{sample}.raw.bam")
    threads: 8
    shell:
        "bwa mem -t {threads} {input.ref} {input.fastq} | samtools view -b > {output}"

rule sort_dedup:
    input:  "results/{sample}.raw.bam"
    output:
        bam="results/{sample}.bam",
        bai="results/{sample}.bam.bai"
    threads: 4
    shell:
        "samtools sort -@ {threads} {input} | samtools markdup -r - {output.bam} "
        "&& samtools index {output.bam}"

rule flagstat:
    input:  "results/{sample}.bam"
    output: "results/{sample}.flagstat.txt"
    shell:  "samtools flagstat {input} > {output}"

ワークフロー 2：SLURM クラスターでの実行

目標：ジョブごとにリソースを割り当て、同じ Snakefile を HPC にデプロイします。

# 1. Create profiles/slurm/config.yaml
mkdir -p profiles/slurm
cat > profiles/slurm/config.yaml << 'EOF'
executor: slurm
jobs: 100
default-resources:
  mem_mb: 4000
  runtime: 60
use-conda: true
latency-wait: 30
rerun-incomplete: true
EOF

# 2. Add resources to compute-heavy rules in Snakefile
# resources:
#   mem_mb=16000, runtime=120

# 3. Submit
snakemake --profile profiles/slurm --cores 256 -n  # dry-run
snakemake --profile profiles/slurm --cores 256      # submit

# 4. Monitor
snakemake --profile profiles/slurm --report report.html  # after completion

主要なパラメーター

パラメーター	コンテキスト	デフォルト	範囲／オプション	効果
--cores	CLI	1	1–N または all	同時実行ジョブ／スレッドの最大数
threads:	Rule	1	1–N	ルールあたりのスレッド数（--cores に応じて拡張）
mem_mb:	resources:	None	整数	メモリ（MB、SLURM プロファイルで使用）
runtime:	resources:	None	整数（分）	ジョブあたりの最大実行時間
--profile	CLI	None	パス	エグゼキューター設定用の YAML プロファイル
--use-conda	CLI	False	フラグ	ルールごとの conda 環境を有効化
--use-apptainer	CLI	False	フラグ	Singularity／Apptainer コンテナーを有効化
-n	CLI	False	フラグ	ドライラン（タスクを表示し、実行しない）
--forceall	CLI	False	フラグ	ステータスに関係なく、すべてのルールを再実行
--rerun-incomplete	CLI	False	フラグ	出力が不完全なルールを再実行
configfile:	Snakefile	None	YAML パス	YAML から設定ディクショナリーを読み込む

ベストプラクティス

1. 必ず rule all を定義する：定義しないと、Snakefile の最初のルールだけが実行されます。rule all はすべての最終出力をまとめ、Snakemake はそれらの生成に必要な処理をすべて実行します。

2. 大きな中間ファイルには temp() を使う：重複除去前の BAM ファイル、ソートされていない BAM ファイル、中間アセンブリには temp() を指定すると、使用後に自動削除されます。これにより、ディスク容量を大幅に節約できます。

3. 設定とコードを分離する：サンプルリスト、スレッド数、ファイルパス、しきい値は config.yaml に記述します。Snakefile に値をハードコードすると、パイプラインの柔軟性と再利用性が損なわれます。

4. まず snakemake -n でテストする：ドライランでは、実行されるルールとその順序が正確に表示されます。DAG が正しいことを確認するため、本番実行の前に必ず実行してください。

5. すべての shell ルールに log: を使う：ツールの stderr／stdout をルールごとのログファイルにリダイレクトします（2> {log}）。ログがないと、クラスターのジョブ失敗時のデバッグはほぼ不可能です。

6. 本番環境でルールのベンチマークを行う：benchmark: "benchmarks/{rule}/{sample}.txt" を追加して、実際の実行時間とメモリ使用量を測定します。これは SLURM のリソース要求を調整するために不可欠なデータです。

よく使うレシピ

レシピ：ファイルからサンプルリストを生成する

# Auto-discover samples from input directory (no hardcoded list)
from pathlib import Path

SAMPLES = [p.stem.replace(".fastq", "") for p in Path("data/").glob("*.fastq")]
print(f"Found {len(SAMPLES)} samples: {SAMPLES[:3]}...")

rule all:
    input:
        expand("results/{sample}.bam", sample=SAMPLES)

レシピ：設定に基づいて条件付きで実行する

configfile: "config/config.yaml"

# Only run deduplication for WGS (not amplicon) data
rule dedup:
    input:
        "results/{sample}.sorted.bam"
    output:
        "results/{sample}.deduped.bam"
    run:
        if config.get("assay_type") == "WGS":
            shell("samtools markdup -r {input} {output}")
        else:
            shell("cp {input} {output}")

レシピ：複数のサンプルを集約する

# Collect all per-sample stats into one summary table
rule multiqc:
    input:
        expand("qc/{sample}_fastqc.zip", sample=SAMPLES),
        expand("results/{sample}.flagstat.txt", sample=SAMPLES)
    output:
        "multiqc/multiqc_report.html"
    shell:
        "multiqc qc/ results/ -o multiqc/"

トラブルシューティング

問題	原因	解決策
AmbiguousRuleException	複数のルールが同じ出力に一致する	wildcard_constraints: を追加する、ruleorder rule_a > rule_b を使う、または出力名を変更する
MissingOutputException	ルールは完了したが、出力ファイルが存在しない	shell の作業ディレクトリーを確認する。出力パスとディスク容量を確認する
TargetFileException	rule all が、どのルールでも生成できないファイルを要求している	expand() の引数がワイルドカード名と一致することを確認する。解決の流れを追跡するには -n を使う
クラスターのジョブがすべて失敗する	ツールに対してリソースが少なすぎる	mem_mb または runtime を増やす。squeue でクラスターのキューを確認する
Conda 環境のビルドに失敗する	パッケージの競合、またはチャンネルの指定ミス	conda-forge を bioconda より前に追加する。パッケージのバージョンを固定する
ルールが予期せず再実行される	出力ファイルのタイムスタンプが入力ファイルより古い	snakemake --touch で出力ファイルのタイムスタンプを更新する。または削除して再実行する
保護された出力で PermissionError が発生する	protected() ラッパーが適用されている	--force で保護を解除する。または削除し、protected() を使わずに再生成する

関連スキル

• samtools-bam-processing — Snakemake パイプラインでよく使われる BAM のソートおよびインデックス作成ルール
• bedtools-genomic-intervals — 下流のアノテーションルールでの区間操作
• neuropixels-analysis — Snakemake が役立つ、複雑な複数ステップのパイプラインの例

参考資料

• Snakemake ドキュメント — ルール、ワイルドカード、プロファイル、API リファレンス
• Snakemake GitHub — ソースコード、リリース、課題トラッカー
• Mölder et al.（2021）「Snakemake による持続可能なデータ解析」— F1000Research 10:33
• Snakemake ワークフローカタログ — コミュニティが管理するリファレンスパイプライン