Snakefile 基本结构：规则、输入、输出、shell 命令

Snakefile 是 Snakemake 的核心配置文件。说白了，它就是告诉 Snakemake 你要做什么、用到哪些数据、产出哪些东西。很多初学者一上来就被 Python 语法和 rule 嵌套吓住，其实根本不用慌。

我个人习惯把 Snakefile 拆成三部分来看：规则（rule）、输入输出（input/output）、以及真正干活的 shell 命令。今天咱们就把这三块讲透。

规则（rule）—— 管线的“零件”

每条 rule 就是流水线上的一道工序。比如在 GATK 管线上，你会有 bwa_alignment、mark_duplicates、base_recalibrator 这样的规则。它们各自只干一件事，组合起来就是完整的变异检测流程。

我曾经见过有人把整个基因组比对、去重、重校正塞进同一条 rule 里……说真的，那调试起来简直要命。所以我的建议是：一条 rule 只做一个原子操作。

rule bwa_alignment:
    input:
        "data/sample_1.fastq.gz"
    output:
        "results/sample_1.sorted.bam"
    shell:
        """
        bwa mem -t 4 ref/hg38.fa {input} | samtools sort -o {output}
        """

注意看，这里的 input 和 output 都是字符串路径，shell 里用花括号 {input}、{output} 引用它们。这是最基础的写法。你想想看，如果把 input 写成 ["data/sample_1.fastq.gz"] 列表也行，但那样在 shell 中引用时需要多一层 {input[0]}，我个人觉得非必要不折腾。

要点：规则名是唯一的，Snakemake 通过规则名和输出文件名来追踪依赖关系。如果两条规则试图生成同一个输出文件，Snakemake 会报错。

输入（input）—— 上游的“半成品”

输入可以是本地文件，也可以是其他规则的输出。Snakemake 会自动根据文件名匹配规则依赖。举个例子：假如 mark_duplicates 这个规则的输入是 "results/sample_1.sorted.bam"，而恰好 bwa_alignment 规则的输出也是这个文件，那 Snakemake 就会自动先跑 bwa_alignment，再跑 mark_duplicates。

rule mark_duplicates:
    input:
        "results/sample_1.sorted.bam"
    output:
        "results/sample_1.dedup.bam"
    shell:
        """
        gatk MarkDuplicates \
            -I {input} \
            -O {output} \
            -M results/sample_1.metrics.txt
        """

这里有个小坑：文件路径必须完全一致，包括目录和扩展名。我之前在项目里因为路径里多了一个 ./，结果 Snakemake 把它当成两个不同的文件，死活不匹配依赖关系。排查了半小时……嗯，从那以后我统一用绝对路径或固定相对路径。

输出（output）—— 你想保留的结果

输出文件是每条 rule 的“产品”。Snakemake 会把输出文件作为下一级规则的输入候选。但要注意：输出文件不需要显式声明所有中间产物。比如你只需要最终的 VCF，那么中间 BAM 文件只要被某个规则作为输入引用，Snakemake 就会自动保留。如果你不引用它，Snakemake 可能会在管线结束后删除它（这由 --delete-all-output 或资源限制决定）。

我习惯把所有重要输出放在 results/ 下，按样本名组织。这样既清晰，又便于并行。对于临时文件，比如 sample_1.sam（BWA 的原始输出），如果后续不需要，我会在 shell 命令里顺手删除，省得占空间。

技巧：可以用通配符 {sample} 让规则适配多个样本。比如 output: "results/{sample}.sorted.bam"，这样一条规则就能处理所有样本，无需重复写。

shell 命令—— 真正干活的代码

shell 块里写的就是你平时在终端敲的命令。Snakemake 会用 Python 的 subprocess 执行它。花括号里可以引用 input、output、params（自定义参数）、wildcards（通配符值）等。

我遇到过一个常见问题：shell 命令太长，比如 GATK 的命令行动辄几十个参数。这时候 不要硬拼成一行，用三对引号包起来，每条参数写一行，反斜杠换行，可读性会好很多。

rule base_recalibrator:
    input:
        bam = "results/{sample}.dedup.bam",
        known_sites = "known_sites/hg38.dbsnp.vcf"
    output:
        "results/{sample}.recal.table"
    params:
        ref = "ref/hg38.fa"
    shell:
        """
        gatk BaseRecalibrator \
            -R {params.ref} \
            -I {input.bam} \
            --known-sites {input.known_sites} \
            -O {output}
        """

你看，我用 input.bam 和 input.known_sites 分别引用两个输入文件。这样语义清晰，以后改输入路径只动 rule 头部，不用改 shell 正文。

注意：shell 命令中的 {input} 如果对应的是列表，会被自动展开成空格分隔的字符串。但如果你需要更复杂的处理，比如遍历列表，那就得用 Python 助手的 expand 函数或者在 shell 里写循环。

一个完整的迷你 GATK 管线示例

下面我把上面几条规则串起来，展示一个从 FASTQ 到重校正 BAM 的 Snakefile 片段。你可以看到规则之间是如何通过文件名自动链接的。

rule all:
    input:
        "results/sample_1.recal.bam"

rule bwa_alignment:
    input:
        "data/sample_1.fastq.gz"
    output:
        "results/sample_1.sorted.bam"
    shell:
        """
        bwa mem -t 4 ref/hg38.fa {input} | samtools sort -o {output}
        """

rule mark_duplicates:
    input:
        "results/sample_1.sorted.bam"
    output:
        "results/sample_1.dedup.bam"
    shell:
        """
        gatk MarkDuplicates -I {input} -O {output} -M results/sample_1.metrics.txt
        """

rule base_recalibrator:
    input:
        bam = "results/sample_1.dedup.bam",
        known_sites = "known_sites/hg38.dbsnp.vcf"
    output:
        table = "results/sample_1.recal.table"
    shell:
        """
        gatk BaseRecalibrator -R ref/hg38.fa -I {input.bam} --known-sites {input.known_sites} -O {output.table}
        """

rule apply_bqsr:
    input:
        bam = "results/sample_1.dedup.bam",
        table = "results/sample_1.recal.table"
    output:
        "results/sample_1.recal.bam"
    shell:
        """
        gatk ApplyBQSR -R ref/hg38.fa -I {input.bam} --bqsr-recal-file {input.table} -O {output}
        """

这个文件简单到只有 50 行不到，但它已经是一个可运行的变异检测前端流程了。你只需要在终端执行 snakemake --cores 4，Snakemake 就会自动判断先跑哪个规则、并行执行哪些步骤。

Snakefile 规则与数据流示例 data/sample.fastq.gz rule bwa_alignment bwa mem → samtools sort sample.sorted.bam rule mark_duplicates MarkDuplicates sample.dedup.bam rule base_recalibrator BaseRecalibrator sample.recal.bam 每条规则接收来自上游的输出文件，并产生新的输出。 Snakemake 通过文件名自动解析依赖关系。

规则组合：input functions 与动态输入

上面展现的静态输入适合小规模项目。但实际 GATK 管线经常会遇到动态输入的情况，比如你要同时用多个已知位点文件。这时候可以用 input functions——它是一个返回路径或路径列表的 Python 函数。

我举个例子：在 ApplyBQSR 步骤中，可能需要把 dbsnp 和 gnomad 两个 VCF 都作为 --known-sites 传进去。如果用静态写法，你得写两遍 --known-sites 参数。聪明的做法是定义一个函数来动态拼接。

def known_sites_input(wildcards):
    # wildcards 包含当前样本名等通配符
    return [
        "known_sites/hg38.dbsnp.vcf",
        "known_sites/hg38.gnomad.vcf"
    ]

rule apply_bqsr:
    input:
        bam = "results/{sample}.dedup.bam",
        table = "results/{sample}.recal.table",
        known_sites = known_sites_input
    output:
        "results/{sample}.recal.bam"
    shell:
        """
        gatk ApplyBQSR \
            -R ref/hg38.fa \
            -I {input.bam} \
            --bqsr-recal-file {input.table} \
            {input.known_sites} \
            -O {output}
        """

这样当 Snakemake 执行 apply_bqsr 时，会自动调用 known_sites_input 函数，传入当前的 wildcards 对象。函数返回的列表会被展开成两个文件路径，然后依次被放入 shell 命令中。

小建议：input function 是 Snakemake 最强大的特性之一。如果你发现某个规则的输入随着样本或参数变化，就应该写一个函数，而不是把所有可能文件都列在静态输入里。

一种常见的误解：shell 与 run 的区别

有些教程会告诉你还可以用 run: 或 script: 来执行 Python 代码或外部的 R/Python 脚本。我个人在 GATK 管线上几乎只用 shell:，因为 GATK 命令本身就是成熟的命令行工具，用 shell 最直接、最易调试。

如果你非得用 run: 来写 Python 逻辑，比如需要根据输入文件名动态生成参数，那也可以。但记得加上 shell 模式的回退——万一 Python 环境里没装 GATK 的 Python 绑定呢？对吧。

嗯，关于 Snakefile 的基本结构，今天就聊这么多。你只要记住：一条规则一个任务，输入输出清晰，shell 命令写完整，搭建可复用的 GATK 管线就不难。

← 上一章 📖 返回目录 下一章 →