snakemake-workflow-engine

Snakemake — Python 工作流引擎

概述

Snakemake 是一个基于 Python 的工作流管理系统，可将分析从笔记本电脑扩展到 HPC 和云端。工作流由规则定义，并显式指定输入／输出文件依赖关系；Snakemake 会自动确定执行顺序，并并行运行相互独立的步骤。规则可以调用 shell 命令、Python／R／Julia 脚本，或内联 Python。每条规则都可使用独立的 conda 或 Singularity 环境，让工作流完全可复现。它广泛用于生物信息学领域的 NGS、基因组组装和变异检测流程。

适用场景

• 构建可复现的多步骤生物信息学流程（比对 → 排序 → 变异检测 → 注释）
• 无需修改代码，即可将同一工作流从本地开发扩展到 SLURM 集群
• 使用基于通配符的规则，以相同方式处理多个样本
• 自动管理依赖关系——只重新运行输入已更改的步骤
• 为每条规则部署独立的 conda 或 Singularity 环境，隔离工具
• 在投入计算资源前生成可视化 DAG 和试运行预览
• 如果需要 Groovy DSL + 数据流通道，或要使用 nf-core 社区流程库，请改用 Nextflow
• 对于简单的 shell 循环，请使用 bash 脚本；只有在存在 3 个以上带分支的连续步骤时，使用 Snakemake 才值得承担额外开销
• 对于包含动态任务图或基于时间调度的数据工程工作流，请改用 Prefect 或 Airflow

前置条件

• Python 包：snakemake、graphviz（用于 DAG 可视化）
• 环境：Python 3.11+；建议使用 conda/mamba 管理每条规则的环境
• 数据要求：输入文件、参考文件；在规则中定义输出路径

安装前请检查：当前环境中可能已提供该工具（例如在 pixi / conda 环境内）。先运行 command -v snakemake；如果返回路径，请跳过下面的安装命令。在 pixi 项目中运行时，请使用 pixi run snakemake 调用该工具，而不是直接运行 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：规则定义

每条规则都定义一个分析步骤，包括输入、输出和执行方法。

# 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：通配符与模式展开

通配符让一条规则可以处理任意数量的样本；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 质控流程

目标：对多个样本执行 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	命令行	1	1–N 或 all	最大并发作业/线程数
threads:	规则	1	1–N	每条规则的线程数（最多扩展至 --cores）
mem_mb:	resources:	无	整数	内存，单位为 MB（供 SLURM 配置文件使用）
runtime:	resources:	无	整数（分钟）	每个作业的最长运行时间
--profile	命令行	无	路径	用于执行器配置的 YAML 配置文件
--use-conda	命令行	False	标志	为每条规则启用 conda 环境
--use-apptainer	命令行	False	标志	启用 Singularity/Apptainer 容器
-n	命令行	False	标志	试运行（显示任务，不执行）
--forceall	命令行	False	标志	无论状态如何，都重新运行所有规则
--rerun-incomplete	命令行	False	标志	重新运行输出不完整的规则
configfile:	Snakefile	无	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 等人（2021）“使用 Snakemake 进行可持续数据分析” — F1000Research 10:33
• Snakemake 工作流目录 — 由社区维护的参考流水线