4. 生物信息学典型工作流剖析
变异检测（Variant Calling）是重测序分析中最具代表性的工作流之一。它把测序仪产出的原始碱基序列，逐步转换为“样本在参考基因组哪些位置存在何种变异”的结构化结论。这条链路之所以值得单独剖析，是因为它同时具备三类典型特征：数据形态多次转换（FASTQ→BAM→VCF）、任务间存在强依赖（比对必须先于排序，去重必须先于变异检测）、可并行粒度不均匀（按样本并行容易，按区间并行需要额外合并）。理解这条链，等于掌握了大多数生信工作流编排的通用骨架。
一、任务链的六个阶段与数据形态
以人类全基因组重测序（WGS，约 30×）为例，从原始数据到变异结果通常经过以下阶段：
- 质控与修剪：输入
sample_R1.fastq.gz/sample_R2.fastq.gz，用 FastQC 评估碱基质量与接头污染，用 fastp 或 Trimmomatic 切除接头和低质量尾端，输出修剪后的 FASTQ。 - 比对：用 BWA-MEM 将 reads 比对到参考基因组（如 GRCh38），输出 SAM；经 samtools 转为 BAM。此阶段决定后续所有坐标的正确性。
- 排序与去重：samtools sort 按坐标排序，Picard MarkDuplicates 或 samtools markdup 标记 PCR 重复，避免重复 reads 造成假阳性变异。
- 碱基质量重校准（BQSR）：GATK BaseRecalibrator 基于已知变异位点学习系统误差，ApplyBQSR 输出重校准后的 BAM。这一步对 Illumina 数据尤其重要。
- 变异检测：GATK HaplotypeCaller 在活跃区域重建单倍型并输出 GVCF；多样本场景下用 GenomicsDBImport 合并，再 GenotypeGVCFs 联合分型。
- 过滤与注释：VariantFiltration 按 QD、FS、MQ 等指标硬过滤，或用 VQSR 做软过滤；随后用 SnpEff/VEP 注释功能影响。
关键点在于：每一步的输入都是上一步的输出，且中间文件体积巨大（30× WGS 的 BAM 常达 80–100 GB）。这意味着编排系统必须处理大文件传递、断点续跑和磁盘清理，而不能只关心命令是否执行成功。
二、依赖关系的三种类型
把上述阶段抽象出来，任务依赖可分为三类，它们对调度策略的要求完全不同：
- 线性依赖：比对→排序→去重→BQSR，严格串行，无法并行。这类依赖决定了工作流的关键路径长度。
- 分叉-汇聚依赖：HaplotypeCaller 可按染色体或区间拆分为多个任务并行执行，各自产出局部 GVCF，再汇聚到 GenomicsDBImport。分叉数由区间划分决定，汇聚点必须等待全部分支完成。
- 跨样本依赖：联合分型要求所有样本的 GVCF 就绪后才能开始。这是典型的“多对一”汇聚，也是队列最容易积压的位置。
下面这张图把六个阶段、数据形态转换与依赖类型放在同一条链路上，便于对照理解。
三、编排视角下的三个工程要点
第一，资源画像要按阶段区分。比对阶段是 CPU 与内存密集型，BWA-MEM 单样本常需 8–16 核、20–30 GB 内存；排序阶段是 I/O 密集型，瓶颈在磁盘读写；HaplotypeCaller 内存需求随区间增大而上升，单区间建议不超过 50 Mb。若用统一资源模板调度所有任务，要么浪费，要么频繁 OOM。合理做法是为每个任务类型定义独立的 CPU、内存、磁盘与预期运行时长，并据此设置重试上限。
第二，中间文件的命名与生命周期要显式管理。推荐采用 {sample}.{stage}.{version}.bam 这类可追溯命名，并在任务定义中声明哪些中间文件在下一阶段完成后可删除。例如 BQSR 完成后，去重 BAM 若不再被其他分支引用即可清理，否则磁盘会在多样本并行时迅速耗尽。
第三，失败重试要区分“可重试”与“不可重试”。因节点抢占、临时网络抖动导致的失败可自动重试；而参考基因组索引缺失、FASTQ 校验和不匹配、VCF 样本名冲突等属于确定性错误，重试只会浪费资源。编排系统应通过退出码与日志关键字识别这两类失败，并给出明确的人工介入提示。
四、一个可直接参考的任务依赖清单
下表把每个任务的输入、输出、依赖与并行策略列清楚，可作为工作流定义文件的填写依据。字段只保留编排真正需要的部分。
| 任务 | 输入 | 输出 | 前置依赖 | 并行策略 |
|---|---|---|---|---|
| fastp | R1/R2 FASTQ.gz | trimmed FASTQ | 无 | 按样本并行 |
| BWA-MEM | trimmed FASTQ + 参考索引 | SAM | fastp | 按样本并行 |
| samtools sort | SAM | sorted BAM | BWA-MEM | 按样本并行 |
| markdup | sorted BAM | dedup BAM | sort | 按样本并行 |
| BQSR | dedup BAM + 已知位点 | recal BAM | markdup | 按样本并行 |
| HaplotypeCaller | recal BAM + 区间 | 局部 GVCF | BQSR | 按区间并行 |
| GenomicsDBImport | 全部样本 GVCF | GenomicsDB | 所有 HaplotypeCaller | 汇聚，不可并行 |
| GenotypeGVCFs | GenomicsDB | 原始 VCF | GenomicsDBImport | 可按区间并行 |
| VQSR / 硬过滤 | 原始 VCF | 过滤 VCF | GenotypeGVCFs | 单任务 |
| VEP 注释 | 过滤 VCF | 注释 VCF | 过滤 | 可按染色体并行 |
把这张表与前面的链路图对照，可以看到编排的核心工作并不是“把命令串起来”，而是为每个任务声明正确的依赖、资源与并行边界。一旦依赖声明错误，轻则重复计算，重则产出坐标错位或样本混淆的 VCF，而这类错误往往在最终注释阶段才暴露，排查成本极高。因此，在把工作流投入生产前，建议先用一个小样本（如 chr20 或 chr22）跑通全链路，确认每个中间文件的样本名、参考版本与坐标一致，再扩展到全基因组与多样本规模。