📚 课程目录第2章工作流基础

规则、输入输出与通配符——Snakemake 的“灵魂三件套”

讲批量对接流程，第一步就是把“要对什么、怎么对、对完生成什么”给说清楚。Snakemake 里，这个“说清楚”的工具就是 规则（rule）。我个人习惯把规则看成 Snakemake 的最小执行单元——每个规则只做一件事，但做得明白。

规则定义：写好你的“配方”

一个规则，说白了就是一个“配方”。它告诉 Snakemake：你需要什么文件进来，你想用什么脚本或命令处理，以及你希望生成什么文件。

举个例子。最简单的规则长这样：

rule prepare_protein:
    input:
        "raw/uniprot.fasta"
    output:
        "prepared/targets.fasta"
    shell:
        "grep '>' {input} > {output}"

你看，这个规则做了什么？它只做了一件事：从原始 FASTA 里提取序列头信息。输入、输出、shell 命令，三个部分清清楚楚。

规则的核心元素其实就这几个：

input——你需要的依赖文件
output——你要生成的文件
shell 或 run——具体执行什么操作

别小看这个结构。我在项目中遇到过有人把十几个处理步骤写在一个规则里，结果出了问题根本不知道在哪一步。你想想看，拆成多个规则，每个规则只做一件事，调试起来舒服多了。

规则元素	作用	常用写法
input	声明依赖文件	字符串、列表、字典
output	声明生成文件	字符串（推荐带通配符）
shell	执行 Shell 命令	用 {} 引用 input/output
run	执行 Python 代码	直接写 Python 逻辑

提示：shell 里引用文件时用 {input} 或 {output}，Snakemake 会自动替换成你定义好的路径。别硬编码路径，否则失去可复用性。

输入输出声明：让 Snakemake 知道谁依赖谁

Snakemake 的工作方式，是“由终及始”的。什么意思？你告诉它“我要生成这个最终文件”，它会自己倒推：生成这个文件需要哪些中间文件？生成中间文件又需要哪些原始数据？

这依赖的就是规则之间的输入输出连接。

看个例子：

rule dock_one:
    input:
        protein="prepared/targets.pdbqt",
        ligand="ligands/compound_1.pdbqt"
    output:
        "results/compound_1_out.pdbqt"
    shell:
        "vina --receptor {input.protein} "
        "--ligand {input.ligand} "
        "--out {output}"

这里的 input.protein 和 input.ligand 都是命名输入。命名输入的好处是什么？当你有多个输入文件时，用名字引用比用索引更清晰。嗯，这里要注意：如果输入只有一个文件，直接用 {input} 就够了，没必要硬加个名字。

输出文件名的命名方式也值得聊两句。我个人建议输出路径里带上配体 ID 或蛋白质 ID，这样你一眼就能看出这个结果是谁的。

警告：输入文件必须真实存在，Snakemake 会先检查它们。如果文件缺失，规则不会执行，而是直接报错。我曾经因为这个特性排查过一整天“为什么流程卡住了”——其实就是配体文件路径写错了。

通配符：让一条规则适配成千上万个小分子

好了，上面那个 dock_one 规则只能处理一个化合物。但我们要对接 1000 个化合物呢？写 1000 条规则？

不可能。这时候就轮到通配符登场了。

通配符，说白了就是在文件名里用花括号 {} 定义一个“变量”。Snakemake 会从你指定的最终目标文件里，自动提取这个变量的值，然后替换到规则中。

比如：

rule dock_all:
    output:
        "results/{ligand_id}_out.pdbqt"
    input:
        protein="prepared/targets.pdbqt",
        ligand="ligands/{ligand_id}.pdbqt"
    shell:
        "vina --receptor {input.protein} "
        "--ligand {input.ligand} "
        "--out {output}"

你猜怎么着？只要你告诉 Snakemake “我要生成 results/compound_42_out.pdbqt”，它会自动识别出 {ligand_id} 是 compound_42，然后去找 ligands/compound_42.pdbqt。

为什么会这样？因为 Snakemake 把输出文件名上的通配符，映射到了输入文件名的相同通配符上。一个通配符，贯穿整个规则。

通配符占位	实际值举例	匹配文件
`{ligand_id}`	`compound_42`	`ligands/compound_42.pdbqt`
`{protein_name}`	`1a2b`	`proteins/1a2b.pdbqt`
`{system}`	`docking_set_1`	`configs/docking_set_1.yaml`

提示：通配符的名字最好有意义。用 {id} 不如用 {ligand_id}，因为几个月后你回来看代码，能立刻明白这个 id 代表什么。

通配符约束：别让 Snakemake 猜错

通配符好用，但也可能“过度匹配”。你想想看，如果你有两个文件：ligands/compound_42.pdbqt 和 ligands/compound_42_extra.pdbqt，一个规则里的 {ligand_id} 到底会匹配成 compound_42 还是 compound_42_extra？

Snakemake 的默认行为是尽量匹配最短的可能。但如果你不希望它猜错，就需要 通配符约束（wildcard_constraints） 上场了。

约束可以是正则表达式，也可以是简单的字符串匹配规则。比如：

rule dock_with_constraint:
    output:
        "results/{ligand_id}_out.pdbqt"
    input:
        ligand="ligands/{ligand_id}.pdbqt"
    wildcard_constraints:
        ligand_id="compound_\d+"
    shell:
        "vina --ligand {input.ligand} --out {output}"

这个约束意味着：{ligand_id} 只能匹配“compound_”后面跟一串数字的形式。像 compound_42_extra 这种带有额外字符的，就匹配不上了。

约束不只在规则里写，也可以在 Snakefile 全局定义：

wildcard_constraints:
    ligand_id = "compound_\d+"
    protein_id = "[a-z0-9]{4}"

全局定义的约束，会应用到所有规则的同名通配符上。我个人喜欢在规则里单独写约束，因为这样更清晰，不会让全局文件变得臃肿。

注意：我曾经因为没有加通配符约束，导致流程意外匹配了错误的文件——一个叫 compound_42_backup.pdbqt 的文件被当成了配体。结果对接后的结果全不对，排查了整整一个下午。从那以后，只要文件名不唯一，我一定加约束。

知识结构小结

让我用一张图把这一章的核心逻辑串起来。不用死记硬背，看完这张图你大概就明白了：

这张图的核心逻辑是：规则定义输出 → 输出带通配符 → 输入自动映射同一通配符 → 通配符约束限制匹配范围 → 最终执行。记住这个链条，后面写复杂流程时就不会乱。

核心要点：

每个规则只做一件事，输入输出清晰
通配符让一个规则能处理所有配体
通配符约束防止文件名误匹配
Snakemake 由终及始自动推导依赖关系

如何验证你的理解？

动手试试吧。假设你要对接 5 个化合物，试试用通配符规则来执行。我建议你先从一个规则开始，确认能跑通，再扩展到多个配体。批量流程不是一步到位的——就像搭积木，一块一块搭，稳当。