第一个Snakemake工作流：编写简单的质控规则

终于到了动手环节。前面我们把环境搭好、概念理清，现在来写第一个真正的 Snakemake 工作流。

我个人习惯，学一个新工具时先让它跑通一个最简单的任务。不要一上来就搞复杂的全流程。先能输出一个文件，再考虑后续的批量、条件、报告。这叫「最小可行工作流」。

3.1 第一个 Snakefile 长什么样？

Snakemake 的核心文件叫 Snakefile。你想想看，它其实就是一个带规则的 Python 文件。所有规则（rule）都写在这个文件里。

我们直接来看一个最简单的 Snakefile：

rule hello:
    output:
        "hello.txt"
    shell:
        "echo 'Hello, Snakemake!' > hello.txt"

这个规则做了三件事：

• 声明规则名：rule hello:，名字叫 hello
• 定义输出：output: "hello.txt"，告诉 Snakemake 我要生成一个文件
• 执行命令：shell: 后面跟的就是实际要跑的 shell 命令

运行它试试看：

snakemake -j1

-j1 表示只用一个 CPU 核心。输出长这样：

Building DAG of jobs...
Using shell: /bin/bash
Provided cores: 1
Rules claiming more threads will be scaled down.
Job counts:
	count	jobs
	1	hello
	1

[Mon Jan  1 00:00:00 2024]
rule hello:
    output: hello.txt
    jobid: 0
    reason: Missing output files: hello.txt
    wildcards: --

Finished job 0.
1 of 1 steps (100%) done

看到「100% done」了吗？工作流跑完了。当前目录下多了一个 hello.txt，内容就是 Hello, Snakemake!。

3.2 编写第一条 QC 规则

光输出一个文本文件太简单了。我们用 FastQC 来做真正的质控。

假设你有一个 FASTQ 文件 sample1_R1.fastq，在项目的 data/ 目录下。我们的目标是让 Snakemake 自动对 sample1_R1.fastq 运行 FastQC，并把结果输出到 results/fastqc/。

先写好 Snakefile：

rule run_fastqc:
    input:
        "data/sample1_R1.fastq"
    output:
        "results/fastqc/sample1_R1_fastqc.html",
        "results/fastqc/sample1_R1_fastqc.zip"
    shell:
        "fastqc {input} -o results/fastqc/"

这里我加入了 input:。Snakemake 会检查输入文件是否存在。如果不存在，它会报错并停下，不会傻傻地跑出一个空结果。

运行命令不变：

mkdir -p results/fastqc
snakemake -j1

3.3 理解 Snakemake 的输入输出模式

刚才这个规则虽然能跑，但你仔细看：input 和 output 都是硬编码的。这样写有个大问题——如果换一个样本，就要改 Snakefile。

Snakemake 真正的威力在于「自动推断」。你只需要告诉它最终想要什么，它会反向推导出需要哪些输入、跑哪些规则。

看这个更通用的写法：

rule run_fastqc:
    input:
        "data/{sample}_R1.fastq"
    output:
        "results/fastqc/{sample}_R1_fastqc.html",
        "results/fastqc/{sample}_R1_fastqc.zip"
    shell:
        "fastqc {input} -o results/fastqc/"

注意到 {sample} 了吗？这叫通配符（wildcard）。当你执行：

snakemake results/fastqc/sample1_R1_fastqc.html -j1

Snakemake 会自动把 {sample} 替换成 sample1，然后推导出输入文件是 data/sample1_R1.fastq。

3.4 运行前检查：干运行模式

你写的规则到底对不对？Snakemake 有一个非常实用的功能——干运行（dry-run）。

snakemake -n

-n 参数表示只模拟执行，不真的跑。输出像这样：

Building DAG of jobs...
Job stats:
rule        count
run_fastqc      1

[Tue Jan  2 00:00:00 2024]
rule run_fastqc:
    input: data/sample1_R1.fastq
    output: results/fastqc/sample1_R1_fastqc.html,
            results/fastqc/sample1_R1_fastqc.zip
    jobid: 0
    reason: Missing output files: results/fastqc/sample1_R1_fastqc.html
    wildcards: sample=sample1
    resources: tmpdir=/tmp

Job counts:
	count	jobs
	1	run_fastqc
	1

你会看到它准备执行哪些规则、输入输出是什么、通配符的值是什么。这些信息在调试时非常有用。

3.5 查看结果

工作流跑完后，结果在 results/fastqc/ 目录下。

ls results/fastqc/
sample1_R1_fastqc.html  sample1_R1_fastqc.zip

你可以用浏览器打开 HTML 文件查看质控报告。如果是在远程服务器上，可以用 scp 下载到本地，或者用 python -m http.server 临时起一个 HTTP 服务再访问。

至此，你完成了第一个真正的 Snakemake 工作流。它很小，但已经包含了所有核心要素：规则、输入、输出、通配符、干运行。后面的内容就是在这个基础上不断扩展功能。