Snakemake 基本语法与规则

说实话，刚开始接触 Snakemake 的时候，我脑子里全是“这不就是 Makefile 加了个 Python 壳吗？”做了几个项目之后我才发现，这套规则系统远比想象中灵活。尤其是跟 Seurat 这类 R 包搭在一起的时候，Snakemake 的规则结构能帮我们把单细胞分析拆得干干净净。

规则结构：一个规则就是一个分析步骤

在 Snakemake 里，每个分析步骤都封装在一条 rule 里。说白了，就是告诉管道：你需要什么文件进来，要生成什么文件出去，中间跑什么命令。

一条最基础的长这样：

rule seurat_qc:
    input:
        "data/counts/cellranger/SCRNA_01/outs/filtered_feature_bc_matrix.h5"
    output:
        "results/qc/SCRNA_01_seurat_qc.rds"
    shell:
        "Rscript scripts/run_qc.R {input} {output}"

我个人习惯把规则名写成 项目_步骤 的格式，比如 seurat_qc、seurat_normalize。你想想看，十个规则放在一个 Snakefile 里，名字不乱才怪。我用这种命名方式，扫一眼就知道这一步在干什么。

规则的核心字段我列一下：

• name：规则名，唯一标识
• input：输入文件，可以是一个或多个
• output：输出文件，Snakemake 靠它判断要不要重跑
• shell / run / script：执行什么命令或脚本
• params：额外参数，比如阈值、版本号
• log：日志文件路径
• threads：分配多少线程
• resources：内存、硬盘等资源声明

我在项目中遇到过一个问题：有人把输入文件写成了通配符，结果 Snakemake 把所有样本都当成同一个文件来解析，跑出来全是错的。记住，输入文件要么是具体路径，要么通配符跟输出保持一致。

输入输出：管道的骨架

Snakemake 的依赖关系是靠输入输出推断的。你不需要手动写“先跑 A 再跑 B”，只需要告诉它：B 需要 A 的输出作为输入，它自己就会排顺序。

看个典型的多步传递：

rule normalize:
    input:
        "results/qc/{sample}_seurat_qc.rds"
    output:
        "results/norm/{sample}_seurat_norm.rds"
    shell:
        "Rscript scripts/run_normalize.R {input} {output}"

rule cluster:
    input:
        "results/norm/{sample}_seurat_norm.rds"
    output:
        "results/cluster/{sample}_seurat_cluster.rds"
    shell:
        "Rscript scripts/run_cluster.R {input} {output}"

你看，连中间文件名都不用记。Snakemake 会自动把 normalize 的产物喂给 cluster。嗯，这里要注意：输入输出的文件路径一旦定下来就别轻易改，否则所有下游规则都要动。我吃过这个亏，改了一次路径名，调试了整整一个下午。

关于多输入输出的写法，我用表格总结一下：

有一类特殊输入叫 expand，它能帮我们把多个样本组合成一条规则。我经常用它来合并 Seurat 对象：

rule merge_samples:
    input:
        expand("results/norm/{sample}_seurat_norm.rds", sample=SAMPLES)
    output:
        "results/merged/merged_seurat.rds"
    shell:
        "Rscript scripts/merge.R {input} {output}"

这里 SAMPLES 是一个全局变量，存所有样本名。这样写的好处是，新增一个样本只需要改列表，不用改规则。

通配符：让规则动起来

通配符是 Snakemake 最让我觉得“爽”的特性。你可以把一条规则写成模板，Snakemake 根据实际需要自动填充不同的样本名或参数。

通配符用花括号包裹，比如 {sample}、{step}、{date}。注意，通配符的名字必须跟输出文件里的占位符一致，否则 Snakemake 不知道替换什么。

rule qc_by_sample:
    output:
        "results/{sample}/qc_report.html"
    input:
        "data/{sample}/raw.h5"
    shell:
        "Rscript scripts/qc.R {input} {output}"

为什么这个设计很聪明？你想象一下，如果你有 10 个样本，传统做法要么写 10 条规则，要么写一个循环脚本。Snakemake 只需要这一条规则，它会在执行时根据用户指定的最终目标反推出所有通配符的值。

不过，通配符也有坑。我曾经踩过的一个坑是：通配符在输入和输出中不一致。比如输出里用 {sample}，输入里却写成了 {sample_id}。Snakemake 会报错说通配符不匹配。检查这个问题最好的办法就是多看日志，日志里会清楚告诉你每条规则用哪些通配符值。

通配符还能组合使用：

rule multi_dimension:
    output:
        "results/{tissue}/{treatment}/{sample}.rds"
    input:
        "data/{tissue}/{treatment}/{sample}.h5"
    shell:
        "Rscript scripts/process.R {input} {output}"

这里用了三个通配符，Snakemake 会为每个组合生成一条任务。如果你有 3 种组织、2 种处理、5 个样本，那就是 30 条任务。全部自动生成，不需要你手写。

日志：你的第一调试工具

很多新手觉得日志是可有可无的东西。我刚开始也是这么想的，直到有一次管道跑了 6 个小时，最后报了一个奇怪的错——没有日志，我根本不知道错在哪一步。

Snakemake 对日志的支持非常原生。你只需要在规则里加一个 log: 字段，然后重定向标准输出和错误流进去：

rule seurat_pca:
    input:
        "results/norm/{sample}_seurat_norm.rds"
    output:
        "results/pca/{sample}_seurat_pca.rds"
    log:
        "logs/{sample}_pca.log"
    shell:
        "Rscript scripts/run_pca.R {input} {output} &> {log}"

注意这里的 &> {log}，它把标准输出和标准错误都写进日志文件。我曾经试过只写 >，结果错误信息全丢了，排查了三个小时才发现是 R 包版本冲突。

日志文件本身也支持通配符。你想想看，每个样本有自己的日志，诊断起来多方便。出错了直接看对应样本的日志，不用在几百行控制台输出里翻来翻去。

我个人习惯在日志的第一行打印一个时间戳，方便跟系统日志对比：

rule seurat_pca:
    log:
        "logs/{sample}_pca.log"
    shell:
        """
        echo "[$(date)] Start processing {sample}" &> {log}
        Rscript scripts/run_pca.R results/norm/{sample}_seurat_norm.rds results/pca/{sample}_seurat_pca.rds &>> {log}
        echo "[$(date)] Done" &>> {log}
        """

这里用了 &>> 追加写，确保时间戳和 R 输出都在同一个日志文件里。这个小习惯帮我快速定位过很多问题。

说到日志，我还想提一句：日志文件最好跟输出文件放在同一个项目目录结构下，而不是扔到 /tmp 里。我一般会在项目根目录建一个 logs/ 文件夹，按样本名分。这样 ls logs/ 一眼就能看出哪些样本跑完了、哪些报错了。

你如果用一个跑批任务跑 50 个样本，日志管理不好简直是一场灾难。Snakemake 的日志机制加上通配符，天然就帮你把日志分好了——每个样本、每个步骤都有独立的日志文件。

好了，这章的内容就是这些。规则结构、输入输出、通配符、日志，这四个概念是后续所有复杂管道的基础。把这章吃透了，后面讲 Seurat 集成、并行调度、集群提交的时候，你才会觉得顺手。