📚 课程目录第4章 Snakefile的编写

Snakefile 的骨架：规则结构、全局参数与配置

说实话，很多人刚接触 Snakemake 就被它那个「类 Python」的语法唬住了。其实你把它拆开看，Snakefile 就干三件事：定规则、设全局变量、加载配置。把这三块理清楚，你的流程就跑顺了。

一、规则的结构——流程的最小单元

Snakefile 里最重要的概念就是 rule。我刚开始写 Snakefile 的时候，总想一个规则干完所有事，结果调试起来痛苦得要命。后来我学乖了：一个规则只干一件事，干完就传给下一个。

一个典型的规则包含这么几个部分：

rule fastq_trim:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        "results/{sample}_trimmed_R1.fastq.gz",
        "results/{sample}_trimmed_R2.fastq.gz"
    params:
        adapters="adapters.fa",
        min_length=36
    log:
        "logs/{sample}_trim.log"
    threads: 4
    resources:
        mem_mb=2000
    run:
        shell("cutadapt -a {params.adapters} ... -o {output[0]} -p {output[1]} {input[0]} {input[1]} 2> {log}")

看着长是吧？其实就分三块：声明区（input/output）、配置区（params/threads/resources）、执行区（run/shell）。

个人习惯：我永远把 log 放在 output 旁边。为什么？因为一旦跑崩了，你能第一时间找到日志看是哪里炸的。这是我在一个 50 个样本的大项目里摔出来的教训。

二、全局参数——别在规则里写死

你在每个规则里写 threads=4，等哪天想改成 8 线程，就得一个个改。我见过有人改漏了，结果流程跑一半报错，排查了半天才发现是线程数不统一。何必呢？

好在 Snakemake 提供了全局参数：

# 在 Snakefile 开头定义
cores: 8
jobs: 10
restart_times: 3
latency_wait: 60
use_conda: True

# 或者直接在命令行指定
snakemake -c 8 -j 10 --restart-times 3 --latency-wait 60

我个人倾向于把全局参数写在 Snakefile 里。为什么？因为这样别人拿到你的 Snakefile，一眼就能看到你这个流程的「硬性指标」——跑这个项目需要多少核心、允许重试几次、等待文件多久。这些信息对团队协作太重要了。

避坑：我曾经在集群上跑流程，忘了设 latency_wait。结果 NFS 文件系统写入慢，Snakemake 以为文件没生成，直接报错重跑。浪费了 5 个多小时。后来我就学乖了，至少设 30-60 秒。

三、配置文件——把参数从代码里抽出来

真正的项目中，参数是经常变的东西。参考基因组版本、质控阈值、比对软件参数……你把它们硬编码在 Snakefile 里，每换一个数据集就得改代码。这违背了「开发运维一体化」的精神。

Snakemake 支持 YAML 和 JSON 格式的配置文件。我推荐 YAML，可读性更好：

# config.yaml
reference:
  genome: "/data/references/hg38.fa"
  annotation: "/data/references/gencode.v38.annotation.gtf"

trimming:
  adapter: "/data/adapters/truseq.fa"
  min_length: 36
  quality: 20

qc:
  min_reads: 10000
  min_median_quality: 30

resources:
  threads: 8
  memory_gb: 16

然后在 Snakefile 里用 configfile 加载：

configfile: "config.yaml"

rule fastq_qc:
    input:
        "data/{sample}_R1.fastq.gz"
    output:
        "results/{sample}_fastqc.html"
    params:
        threads=config["resources"]["threads"]
    run:
        shell("fastqc -t {params.threads} {input} -o results/")

你想想看，这样改参数多方便？换个数据集，只需要改 config.yaml 里的几行，什么都不用动。我在团队里推的就是这个模式——分析师只改配置，不用碰流程。

四、规则依赖——让流程自己串联起来

Snakefile 最强大的地方，在于它不需要你手动指定规则顺序。它通过 文件名匹配 来自动识别依赖关系。

比如你有两个规则：

rule trim:
    input:
        "data/{sample}_R1.fastq.gz"
    output:
        "results/{sample}_trimmed_R1.fastq.gz"
    shell:
        "cutadapt -o {output} {input}"

rule align:
    input:
        "results/{sample}_trimmed_R1.fastq.gz"
    output:
        "results/{sample}.bam"
    shell:
        "bwa mem reference.fa {input} | samtools sort -o {output}"

align 的 input 跟 trim 的 output 格式一致，Snakemake 就自动知道要先跑 trim，再跑 align。这就是它的核心逻辑——用文件名做契约。

  关键点：output 的文件名模板决定了流程能不能串联起来。名字对不上，流程就断开了。我建议你在写规则前，先把整个流程的「文件名流」画出来（见下面的结构图），再动手写代码。

五、本章知识结构图

下面这张图帮你理清 Snakefile 的完整骨架。从配置入口到规则执行，再到自动串联，一目了然：

六、从零开始的流程三连问

每次我开始一个新的质控项目，都会问自己这三个问题：

输入长什么样？ —— 是双端测序还是单端？文件名有规律吗？比如 sample1_R1.fastq.gz，这就是通配符 {sample} 的用武之地。
中间产物怎么命名？ —— 从 trimmed 到 aligned 到 sorted 到 dedupped，每个阶段用什么后缀？这个得想清楚。我见过有人用 .tmp 当中间文件，结果后面全乱套了。
最终输出是什么？ —— BAM 文件？质控报告？统计表格？Snakemake 只关心你最终要什么，它会倒推需要哪些规则。

把这三问回答了，你的 Snakefile 结构基本就定型了。剩下的就是填充规则的具体内容。

组件	作用	常见陷阱
configfile	加载可修改的参数	路径写绝对还是相对？我建议用相对路径 + 项目根目录占位符
全局参数	定义线程、重试、等待等运行行为	忘了设 latency_wait，在 NFS 上必炸
rule	定义单步处理逻辑	output 通配符必须唯一，否则会冲突
input/output	通过文件名建立依赖关系	文件夹没创建？用 `touch` 或 `directory()`
params	传递参数到 shell 命令	小心引号转义，我每次都用 `lambda wildcards: ...`

我的一个小习惯：每次写完 Snakefile，我都会运行 snakemake --lint。这工具会帮你检查命名冲突、未使用的通配符、格式问题。别等跑流程了才发现语法错，那太浪费时间了。

七、一个完整的迷你样例

上面讲了一堆理论，咱们来点儿实际的。下面是一个「原始数据 → 质控报告 → 统计汇总」的 Snakefile：

configfile: "config.yaml"

rule all:
    input:
        expand("results/qc/{sample}_fastqc.html", sample=config["samples"]),
        "results/summary.tsv"

rule fastqc:
    input:
        "data/{sample}_R1.fastq.gz"
    output:
        html="results/qc/{sample}_fastqc.html",
        zip="results/qc/{sample}_fastqc.zip"
    params:
        threads=config["resources"]["threads"]
    log:
        "logs/{sample}_fastqc.log"
    threads: config["resources"]["threads"]
    shell:
        "fastqc -t {threads} -o results/qc/ {input} 2> {log}"

rule summarize_qc:
    input:
        "results/qc/{sample}_fastqc.html"
    output:
        "results/summary.tsv"
    run:
        # 这里用 Python 写解析逻辑
        with open(output[0], "w") as out:
            out.write("sample\tstatus\n")
            out.write(f"{wildcards.sample}\tPASS\n")

你看，rule all 是入口，它说「我要这两个东西」。Snakemake 一看，_fastqc.html 由 rule fastqc 负责，summary.tsv 由 rule summarize_qc 负责。fastqc 的 input 需要原始数据，Snakemake 就去 data 目录下找。整个链条就这么自动跑起来了。

嗯，这就是 Snakefile 的精华——你只告诉它「要什么」和「怎么造」，它自己会搞定「先做什么」和「后做什么」。