📚 课程目录第3章 Snakefile基础结构

Snakefile 基础结构：规则、输入输出与通配符

写 Snakemake 流水线，说白了就是在组装一套“自动投料机”。

你告诉机器：用什么原料，执行什么工序，得到什么成品。Snakefile 就是这台机器的说明书。

这章我会把 Snakefile 最核心的三个零件拆开讲清楚：规则（rule）、输入输出文件，还有通配符（wildcards）。理解它们，你就能写出可复用的分析流程，而不是重复粘贴代码。

规则：流水线上的工位

每个 rule 就是一个工位。它接受原材料（input），执行操作（shell 或 run），产出成品（output）。

最小可用规则示例：

rule trim_reads:
    input:
        "data/sample1.fastq"
    output:
        "results/sample1_trimmed.fastq"
    shell:
        "fastp -i {input} -o {output}"

你看，非常直白。input 指向原始数据，output 指向经过 fastp 处理后的文件，shell 里写了具体命令。Snakemake 会自动检查 output 是否已经存在，如果存在就跳过，这叫增量执行。

我在项目中遇到过有人把所有规则写成一个大 shell 脚本，改个参数就得全量重跑。用 rule 切分后，哪个文件更新就只重跑哪个工位，省下大把时间。

规则本体：rule 的构成要素

一个 rule 可以包含以下字段，我按使用频次排个序：

字段	必选/可选	作用
input	可选（但通常都有）	声明输入文件或通配符模式
output	必选	声明输出文件，Snakemake 据此决定是否执行
shell	可选（与 run 二选一）	shell 命令字符串，支持 {input} {output} 占位符
run	可选	Python 代码块，灵活度更高
params	可选	额外参数，可在 shell 中用 {params.name} 引用
threads	可选	指定线程数，控制并发资源
log	可选	日志文件路径，方便调试

小技巧：output 里写两个文件，Snakemake 会等到两个都生成才判定已完成。不信你试试 output: "a.txt", "b.txt"，中途杀掉进程只留一个，下次重跑会重新生成两个。

输入输出：连接工位的传送带

规则之间的依赖靠文件路径自动推导。你不需要手动写"等上一步完成才执行下一步"。Snakemake 通过 input/output 的路径匹配形成 DAG（有向无环图），然后按拓扑序执行。

串联示例：质控 + 比对

rule trim:
    input:
        "data/{sample}.fastq"
    output:
        "results/trimmed/{sample}.fq"
    shell:
        "fastp -i {input} -o {output}"

rule align:
    input:
        "results/trimmed/{sample}.fq"
    output:
        "results/bam/{sample}.bam"
    shell:
        "bwa mem ref.fa {input} | samtools sort -o {output}"

当你执行 snakemake results/bam/sampleA.bam 时，Snakemake 会自动找到 trim 规则先跑，因为 align 的 input 依赖于 trim 的 output。这就叫隐式依赖。

我曾经踩过一个坑：两个规则的 output 路径一模一样，结果 Snakemake 报错说输出文件冲突。记住，output 路径必须全局唯一，否则不知道谁生成了谁。

命名规则与路径模式

我习惯把规则按处理阶段命名：trim、align、call_variants。output 路径最好带上阶段名和样本名，比如 results/bam/{sample}.bam，这样方便后期批量管理和清理。

路径里可以用 / 建子目录，Snakemake 会自动创建不存在的目录。你想想看，这省了多少 mkdir 的麻烦。

通配符：一份配方打天下

通配符是 Snakefile 的魔法。它用 {wildcard_name} 的形式嵌入文件名，让一条规则可以处理任意数量的样本。

通配符匹配示例：

rule count_reads:
    input:
        "data/{sample}.fastq"
    output:
        "counts/{sample}.count.txt"
    shell:
        "echo $(wc -l < {input}) > {output}"

执行 snakemake counts/sample1.count.txt counts/sample2.count.txt，Snakemake 会执行两次规则，一次传入 sample=sample1，一次传入 sample=sample2。通配符的值从目标文件名里自动提取出来。

说白了，通配符让规则从“一次只能处理一个文件”升级成“一次处理一组文件”。

通配符约束：避免过度匹配

有时通配符会匹配到你不想要的路径。比如 data/{sample}.fastq 可能匹配 data/subdir/sample.fastq 吗？默认不会，因为 {sample} 只匹配单个目录层级。如果想匹配跨层级，用 {sample: \\w+} 是行不通的，得用 {sample: [^/]+} 这种正则。

注意：通配符默认匹配 非斜杠字符（除非你用 {sample: .*} 匹配全路径）。我建议尽量限制匹配范围，避免不小心吞掉相邻目录。比如想匹配 data/raw 下的文件，就写 data/raw/{sample}.fastq，别写 data/{sample}.fastq 然后指望只匹配 raw 目录。

多通配符与组合

一条规则可以有多个通配符，比如同时区分样本和染色体：

rule split_chrom:
    input:
        "alignments/{sample}.bam"
    output:
        "split/{sample}/{chrom}.bam"
    shell:
        "samtools view -b {input} {chrom} > {output}"

运行时 snakemake split/sampleA/chr1.bam split/sampleB/chr22.bam，Snakemake 会自主排列组合，每个工位只执行自己需要的组合。

知识体系总览

下面这张图概括了 Snakefile 的核心逻辑。从目标文件开始，逆向推导依赖链，通配符帮你实现批量处理。

实战：将质控流程写成 Snakefile

下面我写一个完整的质控流程 Snakefile，包含通配符、多规则串联，以及参数传递。

完整示例（用于双端测序质控）：

rule fastqc_raw:
    input:
        r1 = "data/{sample}_R1.fastq.gz",
        r2 = "data/{sample}_R2.fastq.gz"
    output:
        "qc_raw/{sample}_R1_fastqc.html",
        "qc_raw/{sample}_R2_fastqc.html"
    params:
        threads = 4
    shell:
        "fastqc -t {params.threads} -o qc_raw {input.r1} {input.r2}"

rule trim_galore:
    input:
        r1 = "data/{sample}_R1.fastq.gz",
        r2 = "data/{sample}_R2.fastq.gz"
    output:
        r1_trim = "trimmed/{sample}_R1_val_1.fq.gz",
        r2_trim = "trimmed/{sample}_R2_val_2.fq.gz"
    params:
        quality = 20
    shell:
        "trim_galore --paired --quality {params.quality} -o trimmed {input.r1} {input.r2}"

rule fastqc_trimmed:
    input:
        r1 = "trimmed/{sample}_R1_val_1.fq.gz",
        r2 = "trimmed/{sample}_R2_val_2.fq.gz"
    output:
        "qc_trimmed/{sample}_R1_fastqc.html",
        "qc_trimmed/{sample}_R2_fastqc.html"
    shell:
        "fastqc -o qc_trimmed {input.r1} {input.r2}"

执行 snakemake qc_trimmed/NA12878_R2_fastqc.html -j4，Snakemake 会自动先跑 fastqc_raw，再跑 trim_galore，最后跑 fastqc_trimmed。注意看，trim_galore 的输出文件命名依赖于 trim_galore 自己的输出格式，你必须提前了解工具的命名规则，否则通配符匹配不上。

踩坑经验：我曾经因为 trim_galore 的新版本改了输出后缀，导致下游 fastqc_trimmed 的 input 永远找不到文件，报错说“MissingInput”。解决办法是在 output 里写死后缀，或者在 shell 里用 mv 重命名为统一的模式。我倾向后者：在 shell 命令后加 mv trimmed/*_{sample}*.fq.gz trimmed/{sample}.fq.gz 这样输出可控。

执行命令与观察日志

跑 Snakemake 时加 -p 打印每条命令，加 --dry-run 预览但不执行。我每次写新规则都会先干跑一次，看看依赖链是否正确。

snakemake qc_trimmed/NA12878_R2_fastqc.html -j4 -p --dry-run

输出会显示类似这样的计划：

Building DAG of jobs...
Job counts:
    count   jobs
    1       fastqc_raw
    1       trim_galore
    1       fastqc_trimmed
    3

[Mon Day 10:00:00 2024]
rule fastqc_raw:
    input: data/NA12878_R1.fastq.gz, data/NA12878_R2.fastq.gz
    output: qc_raw/NA12878_R1_fastqc.html, qc_raw/NA12878_R2_fastqc.html
    jobid: 2
    reason: Missing output files: qc_raw/NA12878_R2_fastqc.html
    wildcards: sample=NA12878
    ...

看 wildcards 那一行，Snakemake 自动提取了样本名。如果通配符没有按预期匹配，--dry-run 能最快暴露问题。

小结：你该记住什么

规则是一个工位，input/output 定义了前后依赖。
通配符让规则复用，用 {sample} 表示一个可变值。
Snakemake 自动建 DAG，从目标文件逆向推导需要执行哪些规则。
先 dry-run 再跑，可以省掉很多重跑时间。

这些就是 Snakefile 的骨架。接下来你可以往骨架上挂更多功能，比如配置文件、条件分支、报告生成。但无论多复杂，底层的规则、输入输出、通配符这三个零件永远不会变。