规则深入：局部变量、动态文件名与扩展函数

这一章，我们扎进 Snakemake 的「规则内部」。说白了，就是搞清楚每个规则里那几个特殊变量——params、log、benchmark，还有怎么让文件名「活起来」。我自己刚学的时候，差点被expand和multiext绕晕，后来发现其实就是个模板引擎。

局部变量三件套：params、log、benchmark

每条规则运行时，Snakemake 会给你三个「专属小仓库」。它们不是全局的，只属于当前这条规则。

params：给 shell 命令传额外参数

写对接流程时，shell 里经常需要传递阈值、路径前缀、算法选项。直接硬编码？太丑了。用 params 把这些参数抽出来，统一管理。

我举个例子：你的对接脚本需要指定一个打分函数，还有输出目录。

rule docking:
    input:
        receptor = "protein.pdbqt",
        ligand = "ligand.pdbqt"
    output:
        "results/{ligand_name}_docked.pdbqt"
    params:
        score_func = "vina",
        out_dir = "docking_outputs"
    shell:
        "vina --receptor {input.receptor} --ligand {input.ligand} "
        "--score_func {params.score_func} "
        "--out {params.out_dir}/{output[0]}"

params 里可以放字符串、数字、列表，甚至动态计算的值。我在做高通量虚拟筛选时，经常把 CPU 核心数也放进去：

params:
    cpu = 4,
    exhaustiveness = 8

然后在 shell 里用 {params.cpu} 引用。嗯，这样改参数只需改规则头部，不用在长命令里翻找。

log：让日志文件不再迷路

日志管理是很多人忽略的环节。Snakemake 的 log 变量会自动记录每条规则的运行日志。你不需要手动重定向 stdout/stderr。

rule docking:
    input:
        "protein.pdbqt",
        "ligands/{ligand}.pdbqt"
    output:
        "results/{ligand}_docked.pdbqt"
    log:
        "logs/docking_{ligand}.log"
    shell:
        "vina --receptor {input[0]} --ligand {input[1]} "
        "--out {output} &> {log}"

注意！我踩过坑：log 只是定义了一个路径，Snakemake 并不会自动帮你重定向输出。你必须在 shell 里手动写 > {log} 或 2>&1。否则日志文件是空的。

benchmark：性能分析的利器

想看看每条规则跑了多久？内存占用多少？用 benchmark 变量。它会自动生成一个 TSV 文件，记录运行时间、最大内存、CPU 时间等。

rule docking:
    input:
        "protein.pdbqt",
        "ligands/{ligand}.pdbqt"
    output:
        "results/{ligand}_docked.pdbqt"
    benchmark:
        "benchmarks/docking_{ligand}.txt"
    shell:
        "vina --receptor {input[0]} --ligand {input[1]} --out {output}"

运行后，Snakemake 会自动填充基准数据。你什么都不用写。我习惯在流程跑完后，用 pandas 读取所有 benchmark 文件，快速定位慢的步骤。

动态文件名：用花括号让名字活起来

Snakemake 最迷人的地方：文件名可以用 {wildcard} 匹配。注意，花括号在 input/output 里是通配符，在 shell 里却是 Python 格式化语法。别搞混。

比如你想处理 1000 个配体，每个配体一个文件：

rule docking:
    input:
        "receptor.pdbqt",
        "ligands/{ligand}.pdbqt"
    output:
        "results/{ligand}_docked.pdbqt"
    shell:
        "vina --receptor {input[0]} --ligand {input[1]} --out {output}"

{ligand} 就是通配符。当你执行 snakemake results/ligand001_docked.pdbqt 时，Snakemake 会自动把 {ligand} 替换为 ligand001。我最初看到这个设计时觉得：「这不就是 Makefile 的 % 吗？」但 Snakemake 更强：支持多个通配符，还能用正则约束。

通配符约束举例：

rulecomplexity:
    wildcard_constraints:
        ligand = "[a-zA-Z0-9]+"
    input:
        "ligands/{ligand}.pdbqt"
    output:
        "results/{ligand}_docked.pdbqt"
    shell:
        "..."

加约束可以防止奇怪的字符导致匹配失败。我遇到过配体文件名包含下划线被错误分割的情况，加上正则就稳了。

扩展函数：expand 和 multiext

这两个函数是 Snakemake 的「批量生成器」。很多教程把它们放在最后讲，但我觉得应该尽早掌握。

expand：一次生成多个文件名

当你需要在 input 或 output 里列出多个文件时，expand 可以替代手写循环。它会做笛卡尔积。

# 生成所有配体的对接结果文件列表
LIGANDS = ["lig001", "lig002", "lig003"]
RULES = ["vina", "glide"]

rule all:
    input:
        expand("results/{ligand}_{rule}_docked.pdbqt",
               ligand=LIGANDS, rule=RULES)

这会产生 6 个文件名：lig001_vina_docked.pdbqt、lig002_vina_docked.pdbqt……lig003_glide_docked.pdbqt。默认是笛卡尔积，你也可以用 zip 模式做一一对应：

expand("results/{ligand}_docked.pdbqt",
       zip, ligand=LIGANDS)
# 只输出 lig001_docked.pdbqt, lig002_docked.pdbqt, ... 一一对应

multiext：一个文件名变成多个后缀

multiext 用来处理同一基础名、不同后缀的情况。比如你有一个配体文件，想同时产生 .pdbqt、.log、.txt 三个输出。

rule prepare_ligand:
    input:
        "raw/{ligand}.sdf"
    output:
        multiext("prepared/{ligand}", ".pdbqt", ".log", ".txt")
    shell:
        "obabel {input} -O {output.pdbqt} -h 2>&1 | tee {output.log}"

注意：multiext 返回的是一个命名元组，可以通过 {output.pdbqt}、{output.log} 分别访问。比手写三个 output 方便得多。

我自己的习惯：处理文件格式转换时，multiext 简直是天赐之物。以前要写三条规则，现在一条搞定。

SVG 框架图：局部变量与动态文件名关系

下面这张图总结了本章的核心结构。你一看就能明白每个变量在规则中的位置，以及 expand 和 multiext 如何生成文件名。

常见陷阱与避坑指南

这里列几个我实际踩过的坑，你记一下，能省不少调试时间。

实战小贴士：组合使用

最后，我把本章所有知识串起来，写一个完整的对接流程片段：

LIGANDS = ["mol1", "mol2", "mol3", "mol4"]

rule all:
    input:
        expand("results/{ligand}_docked.pdbqt", ligand=LIGANDS),
        expand("logs/docking_{ligand}.log", ligand=LIGANDS),
        expand("benchmarks/docking_{ligand}.txt", ligand=LIGANDS)

rule docking:
    input:
        receptor = "receptor.pdbqt",
        ligand = "ligands/{ligand}.pdbqt"
    output:
        "results/{ligand}_docked.pdbqt"
    params:
        exhaustiveness = 16,
        cpu = 4,
        extra = ""
    log:
        "logs/docking_{ligand}.log"
    benchmark:
        "benchmarks/docking_{ligand}.txt"
    shell:
        "vina --receptor {input.receptor} "
        "--ligand {input.ligand} "
        "--out {output} "
        "--exhaustiveness {params.exhaustiveness} "
        "--cpu {params.cpu} "
        "{params.extra} "
        "&> {log}"