📚 课程目录第4章第一个质控规则

第一个质控规则：编写 FastQC 规则

好，咱们直接开工。第一个要写的规则就是跑 FastQC，对单个样本的 fastq 文件做质量检查。这个步骤在真实项目中几乎是必做的，拿到测序数据第一件事，就是看看它到底干不干净。我习惯把它叫作「开门见山的第一道坎」，你想想看，数据质量不行，后面的分析再花哨也是白费力气。

规则目标：给每个样本打一份质量报告

说白了，我们要让 Snakemake 自动对每个传入的 fastq 文件执行 FastQC，生成报告文件。这里有一个关键点：FastQC 的输出通常有两个文件 – .zip 和 .html。我们的规则可以只关注其中一个，也可以统统捕获，后面到底用哪个取决于你。我个人倾向于把 .zip 也保留，因为它方便后续的 MultiQC 汇总。

先看一眼最基本的 Snakefile 结构：

rule fastqc:
    input:
        "data/samples/{sample}.fastq"
    output:
        "results/fastqc/{sample}_fastqc.html",
        "results/fastqc/{sample}_fastqc.zip"
    run:
        shell("fastqc -o results/fastqc {input}")

注意这里的通配符 {sample}，Snakemake 会根据实际需要解析出来。比如你有一个 fastq 叫 NA12878.fastq，它会把 {sample} 自动设为 NA12878。这个机制很灵活，我刚开始用的时候还总想着手动写死路径，后来才发现通配符才是亲爹。

避坑：output 目录必须提前存在

我曾经犯过一个低级错误：FastQC 已经有 -o 参数了，但如果你指定的 results/fastqc/ 目录还没创建，FastQC 会直接报错退出。解决办法有两个：

在 Snakefile 里加一条 shell("mkdir -p results/fastqc") 前置命令。
或者让 Snakemake 管理目录：把 output 中的文件层级目录写清楚，Snakemake 会自动创建父目录。但注意，FastQC 的 -o 参数本身不会自动创建目录，你依然需要手动保证目录存在。

我推荐第二种做法，更干净：在 run 块里先用 shell("mkdir -p " + os.path.dirname(output[0]))。看起来多一行，但能救命。

提醒：有些系统版本的 FastQC 不认相对路径，建议在 shell 命令里用绝对路径，或者提前用 os.path.abspath 转换一下。

给规则加点“料”：日志和资源控制

实际生产环境里，我们肯定不能光写一个 rule 就跑。你还需要考虑日志输出、线程数和内存限制。不然哪天集群上跑崩了，你连错在哪都不知道。我一般会这样写：

rule fastqc:
    input:
        "data/samples/{sample}.fastq"
    output:
        html  = "results/fastqc/{sample}_fastqc.html",
        zip   = "results/fastqc/{sample}_fastqc.zip"
    log:
        "logs/fastqc/{sample}.log"
    threads: 2
    resources:
        mem_mb = 2048
    run:
        # 确保输出目录存在
        shell("mkdir -p results/fastqc logs/fastqc")
        shell("fastqc -o results/fastqc --threads {threads} {input} 2>&1 | tee {log}")

看到 --threads 了吗？FastQC 在多核下能稍微快一点，但其实它核心瓶颈是 I/O。调成 2 就够，别浪费资源。内存给 2GB 是因为大 fastq 文件解压后很容易吃掉 1.5GB+。我遇到过有的老师说给 512MB 就行，结果跑俩小时被 OOM kill 了……嗯，宁可多给一点。

传递参数：用 params 让配置更清爽

如果你想把输出目录或线程数变成可配置的，别硬编码在 run 块里。Snakemake 提供了 params 字段：

FASTQC_DIR = config.get("fastqc_dir", "results/fastqc")

rule fastqc:
    input:
        "data/samples/{sample}.fastq"
    output:
        html  = f"{FASTQC_DIR}/{{sample}}_fastqc.html",
        zip   = f"{FASTQC_DIR}/{{sample}}_fastqc.zip"
    params:
        outdir = FASTQC_DIR,
        threads = 2
    log:
        "logs/fastqc/{sample}.log"
    run:
        shell("mkdir -p {params.outdir} logs/fastqc")
        shell("fastqc -o {params.outdir} --threads {params.threads} {input} 2>&1 | tee {log}")

这样你只需要改配置文件，不用改规则本身。我曾在 30 多个样本的项目里，因为换集群改了输出目录，只动了一行 config 就全解决了，爽快。

验证输出：Snakemake 如何判断这条规则跑好了？

Snakemake 靠的是 output 文件的存在和修改时间。如果一个 rule 的输出文件全部存在且比输入新，它就不会重跑。这其实是个大坑——有时候 FastQC 报错了，但依然生成了一个残缺的 .html 文件，Snakemake 会误以为成功了。怎么防？

你可以在 run 块末尾加一个检查命令：比如 shell("grep -c 'Per base sequence quality' {output.html} > /dev/null")，如果关键字段不存在，抛异常。
或者直接让 FastQC 用 --noextract 和 --format fastq 强制严格模式。

我个人更倾向于后者，简单粗暴：

shell("fastqc -o {params.outdir} --threads {params.threads} --noextract {input} 2>&1 | tee {log}")

小技巧：如果你想把 zip 文件也保留用于 MultiQC，但不想在 Snakemake 里显式声明其为 output（减少文件管理），可以使用 touch() 作为伪目标，但那样容易造成重跑不一致。我一般老老实实把 zip 放在 output 里，Snakemake 自己会处理好。

本章知识体系图

下面这张图概括了 FastQC 规则在 Snakemake 工作流中的位置和关键环节：

实战小贴士：多条规则之间如何串起来

FastQC 规则写完后，你通常还会跟着一个 MultiQC 汇总规则。如何把它们串起来？最简单的方式就是让 MultiQC 规则的输入依赖于 FastQC 的输出。比如：

rule multiqc:
    input:
        expand("results/fastqc/{sample}_fastqc.zip", sample=SAMPLES)
    output:
        "results/multiqc/multiqc_report.html"
    run:
        shell("multiqc -o results/multiqc results/fastqc")

这里 SAMPLES 是你定义好的样本列表。Snakemake 会自动识别依赖关系，先跑所有 FastQC，再跑 MultiQC。你只管 snakemake -c1 触发就好了。

  核心要点：
  用通配符 {sample} 动态匹配输入输出，别写死路径。
输出目录要提前创建好，或者在 shell 里用 mkdir -p。
加 threads、mem_mb、log 让规则更健壮。
输出文件是否完好需要额外检查（或靠 FastQC 的严格模式）。

好，FastQC 规则就先聊到这里。你把它跑起来，看到每个样本的 _fastqc.html 都生成好，就可以放心地往下走了。下一步我们可以看看怎么把这些质量报告汇总成一个漂亮的 MultiQC 报告——嗯，这个我们后面再细说。