📚 课程目录第4章 Snakemake核心机制

4. Snakemake核心机制

Snakemake 的工作流由若干条 rule（规则） 组成，每条规则描述“如何从一组输入文件生成一组输出文件”。理解规则的关键，是把一条规则看成一份可被 Snakemake 调度的“文件生成契约”：它声明输入、输出、执行命令，以及运行所需的参数、日志、线程和资源。Snakemake 根据规则之间的输入输出文件名匹配关系，自动推导出依赖图，并决定执行顺序。

规则的基本结构

下面是一条最简规则，用于对 FASTQ 做质量统计：

rule fastqc_raw:
    input:
        "data/sample1_R1.fastq.gz",
        "data/sample1_R2.fastq.gz"
    output:
        "qc/sample1_R1_fastqc.html",
        "qc/sample1_R2_fastqc.html"
    threads: 2
    log:
        "logs/fastqc/sample1.log"
    shell:
        "fastqc -t {threads} -o qc {input} > {log} 2>&1"

其中 input 与 output 是规则的核心。Snakemake 通过比较不同规则的输出与输入文件名，自动建立依赖关系。例如，若某条比对规则的输入是 qc/sample1_R1_fastqc.html，Snakemake 就会先执行 fastqc_raw。这种“按文件名连接”的机制，使工作流无需显式写执行顺序。

通配符：让规则适配多个样本

如果为每个样本写一条规则，工作流会迅速膨胀。通配符（wildcards）用占位符代替样本名，使一条规则可处理多个样本。例如：

rule fastqc_raw:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        "qc/{sample}_R1_fastqc.html",
        "qc/{sample}_R2_fastqc.html"
    threads: 2
    log:
        "logs/fastqc/{sample}.log"
    shell:
        "fastqc -t {threads} -o qc {input} > {log} 2>&1"

这里的 {sample} 是通配符。当 Snakemake 需要生成 qc/sample1_R1_fastqc.html 时，会把 sample 解析为 sample1，并据此推导输入文件为 data/sample1_R1.fastq.gz。通配符的匹配是“反向推导”：从目标输出文件名反推通配符取值，再代入输入、日志、参数等位置。

需要注意，通配符不能有歧义。若两条规则都能匹配同一输出文件，Snakemake 会报 AmbiguousRuleException。此时可用 ruleorder 指定优先级，或通过更具体的路径、wildcard_constraints 限制匹配范围。

输入输出与参数

input 和 output 支持命名列表，便于在命令中按名称引用：

rule trim_reads:
    input:
        r1 = "data/{sample}_R1.fastq.gz",
        r2 = "data/{sample}_R2.fastq.gz"
    output:
        r1 = "trimmed/{sample}_R1.trim.fastq.gz",
        r2 = "trimmed/{sample}_R2.trim.fastq.gz"
    params:
        adapter = "AGATCGGAAGAGC",
        minlen = 36
    threads: 4
    log:
        "logs/trim/{sample}.log"
    shell:
        "cutadapt -j {threads} -a {params.adapter} -m {params.minlen} "
        "-o {output.r1} -p {output.r2} {input.r1} {input.r2} > {log} 2>&1"

params 用于声明不参与依赖推导、但影响命令执行的参数，如接头序列、最小长度、阈值等。它与 input 的区别在于：input 是文件依赖，params 是运行参数。把参数集中写在 params 中，便于统一调整，也避免把参数硬编码在 shell 字符串里。

日志、线程与资源声明

log 指定日志文件路径，Snakemake 会在执行前自动创建日志目录。把标准输出和标准错误重定向到 {log}，可以在任务失败时快速定位问题。若省略 log，错误信息会混在 Snakemake 主输出中，不利于排查。

threads 声明任务使用的 CPU 线程数，命令中通过 {threads} 引用。它既影响实际运行，也影响调度：Snakemake 会确保同时运行的任务线程总数不超过 --cores 指定值。resources 用于声明内存、GPU、磁盘等非 CPU 资源，例如：

resources:
    mem_mb = 8000,
    gpu = 1

运行时可用 --resources mem_mb=32000 gpu=1 限制总资源。这样，即使机器有 64 GB 内存，也可避免多个高内存任务同时挤爆内存。

构建第一个可运行流程

下面给出一个完整的最小工作流，包含原始数据质控、接头修剪和汇总报告三个环节。假设样本名为 sample1 和 sample2，原始数据位于 data/。

SAMPLES = ["sample1", "sample2"]

rule all:
    input:
        "qc/multiqc_report.html"

rule fastqc_raw:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        "qc/{sample}_R1_fastqc.html",
        "qc/{sample}_R2_fastqc.html"
    threads: 2
    log:
        "logs/fastqc/{sample}.log"
    shell:
        "fastqc -t {threads} -o qc {input} > {log} 2>&1"

rule trim_reads:
    input:
        r1 = "data/{sample}_R1.fastq.gz",
        r2 = "data/{sample}_R2.fastq.gz"
    output:
        r1 = "trimmed/{sample}_R1.trim.fastq.gz",
        r2 = "trimmed/{sample}_R2.trim.fastq.gz"
    params:
        adapter = "AGATCGGAAGAGC",
        minlen = 36
    threads: 4
    log:
        "logs/trim/{sample}.log"
    shell:
        "cutadapt -j {threads} -a {params.adapter} -m {params.minlen} "
        "-o {output.r1} -p {output.r2} {input.r1} {input.r2} > {log} 2>&1"

rule fastqc_trimmed:
    input:
        "trimmed/{sample}_R1.trim.fastq.gz",
        "trimmed/{sample}_R2.trim.fastq.gz"
    output:
        "qc/{sample}_R1.trim_fastqc.html",
        "qc/{sample}_R2.trim_fastqc.html"
    threads: 2
    log:
        "logs/fastqc_trimmed/{sample}.log"
    shell:
        "fastqc -t {threads} -o qc {input} > {log} 2>&1"

rule multiqc:
    input:
        expand("qc/{sample}_R1_fastqc.html", sample=SAMPLES),
        expand("qc/{sample}_R2_fastqc.html", sample=SAMPLES),
        expand("qc/{sample}_R1.trim_fastqc.html", sample=SAMPLES),
        expand("qc/{sample}_R2.trim_fastqc.html", sample=SAMPLES)
    output:
        "qc/multiqc_report.html"
    log:
        "logs/multiqc.log"
    shell:
        "multiqc qc -o qc > {log} 2>&1"

运行命令为：

snakemake --cores 8 --resources mem_mb=16000 --use-conda

执行时，Snakemake 先解析 rule all 的目标文件 qc/multiqc_report.html，发现它由 multiqc 规则生成；multiqc 的输入又依赖各样本的 FastQC 报告，于是继续向上追溯，最终形成完整依赖图。若某一步失败，日志会写入对应 logs/ 文件，便于定位。

运行前的检查与常见问题

正式执行前，建议先用 snakemake -n -r 做干跑（dry-run），查看将要执行的任务及原因。-n 表示只打印不执行，-r 会显示每个任务被触发的原因。若发现依赖图不符合预期，通常是通配符匹配错误或输出路径写错。

常见问题包括：输出文件未声明导致重复执行；log 目录不存在导致重定向失败；threads 与命令中实际线程数不一致；以及 params 中引用了未定义的变量。排查时优先看 logs/ 下对应日志，再检查规则声明是否完整。

建议把样本列表、参考基因组路径、公共参数写在 config.yaml 中，通过 config["key"] 引用。这样同一套规则可在不同项目间复用，只需替换配置文件，而不必修改规则本身。

至此，一个包含质控、修剪、汇总的最小可运行流程已经建立。它的价值不在于步骤多少，而在于每条规则都清晰声明了输入、输出、参数、日志、线程与资源，使 Snakemake 能够自动推导依赖、并行调度，并在失败时留下可追溯的日志。