3. 基于真实FASTQ的质控流程

测序下机后的第一件事不是比对，而是判断数据本身是否可用。FastQC 负责对单个 FASTQ 文件做逐碱基质量、GC 分布、接头污染、重复序列等检查，MultiQC 则把多个样本的 FastQC 报告聚合成一张总览表。本节围绕真实 FASTQ 文件，把这两步组织成一条可复现的 Snakemake 流程，重点解决三个工程问题：样本如何驱动、输出如何组织、日志如何留痕。

样本列表驱动：让规则与样本解耦

质控流程最容易写死的地方是样本名。如果规则里逐个列出 sample1_R1.fastq.gz、sample2_R1.fastq.gz，每加一个样本就要改 Snakefile，流程很快会失控。正确做法是用一个样本列表文件驱动整个流程。常见格式是两列或三列 TSV：

sample  fq1                     fq2
S1      raw/S1_R1.fastq.gz      raw/S1_R2.fastq.gz
S2      raw/S2_R1.fastq.gz      raw/S2_R2.fastq.gz
S3      raw/S3_R1.fastq.gz      raw/S3_R2.fastq.gz

在 Snakefile 顶部读取该表，构造样本名到文件路径的映射：

import pandas as pd

samples = pd.read_csv("samples.tsv", sep="\t")
SAMPLE = samples["sample"].tolist()
FQ1 = dict(zip(samples["sample"], samples["fq1"]))
FQ2 = dict(zip(samples["sample"], samples["fq2"]))

这样规则里只出现 {sample} 通配符，新增样本只需在 TSV 中追加一行，Snakemake 会自动推导出需要生成的新目标文件。样本列表是流程的“单一事实来源”，也是后续所有规则共享的输入契约。

输出目录组织：按工具与样本分层

质控产物数量多且类型杂，目录结构必须提前规划。推荐按“工具/样本”两级组织，把 FastQC 的 HTML 与 zip 放在 qc/fastqc/{sample}/，MultiQC 的聚合报告放在 qc/multiqc/，日志统一放在 logs/ 下并镜像规则名与样本名。这样做的直接好处是：清理中间结果时可按工具目录整体删除，排查问题时能凭路径反推是哪个规则、哪个样本产生的文件。

FastQC 规则可以写成：

rule fastqc:
    input:
        r1 = lambda w: FQ1[w.sample],
        r2 = lambda w: FQ2[w.sample]
    output:
        html1 = "qc/fastqc/{sample}/{sample}_R1_fastqc.html",
        zip1  = "qc/fastqc/{sample}/{sample}_R1_fastqc.zip",
        html2 = "qc/fastqc/{sample}/{sample}_R2_fastqc.html",
        zip2  = "qc/fastqc/{sample}/{sample}_R2_fastqc.zip"
    log:
        "logs/fastqc/{sample}.log"
    threads: 2
    shell:
        """
        fastqc -t {threads} -o qc/fastqc/{wildcards.sample} \
            {input.r1} {input.r2} > {log} 2>&1
        """

这里有两个细节值得注意。第一，FastQC 的 -o 只指定输出目录，文件名由工具按输入文件名自动生成，因此 output 中声明的文件名必须与 FastQC 实际产物一致，否则 Snakemake 会认为规则未完成。第二，log: 指令把标准输出和标准错误重定向到日志文件，而不是打印到终端，这是后续排障的关键。

MultiQC 聚合与依赖关系

MultiQC 的输入是上一步所有样本的 FastQC zip 文件。它不需要逐个声明，而是扫描目录：

rule multiqc:
    input:
        "qc/fastqc"
    output:
        "qc/multiqc/multiqc_report.html"
    log:
        "logs/multiqc.log"
    shell:
        """
        multiqc qc/fastqc -o qc/multiqc -n multiqc_report.html \
            > {log} 2>&1
        """

把目录作为 input 会让 Snakemake 等待该目录下所有由其他规则生成的文件就绪，从而自动建立“先 FastQC、后 MultiQC”的依赖顺序，无需手工写 ruleorder 或显式列出每个 zip。最终用一条 rule all 汇总目标：

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

执行 snakemake -j 8 --use-conda 时，Snakemake 会先并行跑完所有样本的 FastQC，再触发 MultiQC。若某个样本的 FASTQ 缺失，流程会在该样本处报错并停止，而不会静默跳过。

日志记录：把“跑过了”变成“可追溯”

Snakemake 的 log: 指令不只是把输出存盘，它让每个规则的运行痕迹与样本一一对应。FastQC 的日志里会记录处理的读段数、模块执行情况；MultiQC 的日志会记录扫描到多少个 FastQC 报告、哪些样本被纳入。当 MultiQC 报告里某个样本缺失时，第一反应应该是查 logs/fastqc/{sample}.log，而不是重新跑一遍。

一个常见的坑是：FastQC 对损坏的 gzip 文件可能返回非零退出码，但 Snakemake 默认会因非零退出而终止。如果希望流程继续并保留问题样本的日志，可以在规则中加 shell 的容错处理，或在 log: 之外用 benchmark: 记录资源消耗。对于教学与生产都推荐的做法是：让流程在异常处失败，保留日志，人工确认后再决定是否跳过该样本，而不是让错误被静默吞掉。

至此，一条由样本列表驱动、输出分层、日志留痕的质控流程就成型了。它的价值不在于命令本身，而在于把“哪些样本、经过什么处理、产物在哪里、出错看哪里”这四件事固定成了可重复执行的结构。后续的比对、定量规则都可以沿用同一套样本列表与目录约定，流程的扩展成本会显著降低。