📚 课程目录第3章单个样本GVCF生成

单个样本 GVCF 生成：HaplotypeCaller 实战详解

终于到了真正动手生成 GVCF 的环节。这节课我们只聊一个工具：HaplotypeCaller，而且是它的 GVCF 模式。说白了，就是把每个样本的变异信息，以“概率快照”的形式存下来，为后面的联合基因分型铺路。我最早做大规模 WGS 时，对这一步的理解很肤浅，结果联合分析阶段卡了好几天。嗯，后来才明白——GVCF 质量直接决定最终结果的可靠性。

HaplotypeCaller 的 GVCF 模式

GATK4 的 HaplotypeCaller 默认输出 VCF，但加上 --emit-ref-confidence GVCF 后，就会进入 GVCF 模式。它会针对每个位点输出经过局部重比对的单倍型信息，并在未变异区域生成 block 格式的参考置信度记录，大大减小文件体积。

提示：如果你用 GATK4.0 以上版本，建议始终使用 HaplotypeCaller --emit-ref-confidence GVCF，不要再用老版的 -T HaplotypeCaller --emitRefConfidence GVCF。命令行写法变了，我见过有人混用新旧参数导致报错。

核心参数一览

我把它常用的参数整理成一张表，方便你对照。记住，这里没有“魔改”参数，只有适合你数据场景的合理搭配。

参数	作用	我常用的值 / 说明
`--emit-ref-confidence`	输出模式	`GVCF` 或 `BP_RESOLUTION`。99% 场景选 `GVCF`
`--gvcf-gq-bands`	GQ 合并区间	默认值 `5,20,60`；低深度时可改为 `5,10,20,40,60`
`--max-alternate-alleles`	最大候选等位基因数	默认 6；WES 样本可降到 4，WGS 保持在 6
`--pair-hmm`	HMM 算法实现	推荐 `AVX_LOGLESS_CACHING`（自动检测硬件）
`-XL / -L`	区间限制	WES 用 `-L bed文件`；WGS 可跳过
`--native-pair-hmm-threads`	HMM 线程数	设置为 `4-8`，不要超过物理核数
`-RF`	读取过滤	常用 `--disable-read-filter NotDuplicateReadFilter` 保留重复 (配合 MarkDuplicates 已标记的情况)

一句话总结：GVCF 模式下的核心参数就这几个。我把 80% 的地方都用默认值，只调 -ERC 和 --gvcf-gq-bands。其余参数除非有明确瓶颈，否则别乱动。

输出模式：GVCF vs BP_RESOLUTION

GATK 的 HaplotypeCaller 有两种 GVCF 输出模式：

GVCF（推荐）：产生 block 格式，相同等位基因频率的连续区间会尽量合并。文件小，联合分析时可以高效解压缩。
BP_RESOLUTION：每个碱基都单独输出。体积巨大，几乎只在调试或需要极精细控制时使用。我从来不用。

举个例子，一条 1Mb 的区域，GVCF 模式可能只输出 500 条记录，BP_RESOLUTION 会输出 100 万条。你想想看在 10 万样本级别的项目中，这是个什么量级？

所以记住：非特殊情况，一律用 -ERC GVCF。

常用自定义参数详解

聊几个我经常被问到、也在项目中踩过坑的参数。

1. --gvcf-gq-bands

这个参数定义了多少种 GQ 区间。GQ 就是基因型质量分数，0-99。默认分成三档：5、20、60。比如 GQ = 18 的位点会归到 5-20 区间，用 END 标记合并。为什么要合并？
因为对于非变异位点，精确到每个碱基的 GQ 价值不大，合并后文件更小，又不会丢失关键信息。

我在做深度很低的混池测序（<5X）时，把 bands 改成了 5,10,20,40,60，这样低置信度区间更细，联合分析时能更准确区分真纯合和缺失。不过常规 30X WGS 用默认值就好。

2. --max-alternate-alleles

这个参数控制 HaplotypeCaller 在每个位点最多考虑几种候选等位基因。默认 6。如果你做全外显子组（WES），很多位点深度高但变异少，设成 6 会浪费计算。我一般 WES 设成 4，WGS 保持 6。超过这个数的超多等位基因位点极少，没必要让 HMM 去算。

3. --pair-hmm

这是计算 Pr(D|H) 的核心算法。GATK4 默认用 AVX_LOGLESS_CACHING，比老版的 ORIGINAL 快 2-3 倍。如果你跑在 ARM 架构的机器上，可能需要显式指定 LOGLESS_CACHING。我记得有一次帮同事远程看问题，他的 Hail 集群是 ARM 的，用默认参数跑了一周没跑完，改成这个后三天搞定。

警告：--pair-hmm AVX_LOGLESS_CACHING 需要 CPU 支持 AVX 指令集。老机器（比如 Intel E5-26xx v2 之前）可能无法使用。如果报错 Illegal instruction，回退到 LOGLESS_CACHING。

4. 线程控制

HaplotypeCaller 支持多线程：

-XX:ParallelGCThreads：Java GC 线程，一般设为 2-4 即可
--native-pair-hmm-threads：HMM 计算的 native 线程数，我建议设为物理核数的一半左右。比如 40 核的机器设成 20，因为还有 IO 和其他 Java 进程竞争。

我曾经犯过一个错：把 --native-pair-hmm-threads 设成 64，结果操作系统直接 OOM，因为每个线程占用的内存叠加了。注意：这个线程数 不等于 并行样本数，是单样本内部的并行级别。

一条完整的命令模板

下面是我在 30X WGS 项目中最常用的命令结构：

gatk --java-options "-Xmx16g -XX:ParallelGCThreads=2" \
    HaplotypeCaller \
    -R /ref/hg38.fa \
    -I sample.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF \
    --gvcf-gq-bands 5 20 60 \
    --max-alternate-alleles 6 \
    --pair-hmm AVX_LOGLESS_CACHING \
    --native-pair-hmm-threads 8 \
    --tmp-dir /tmp/gatk

注意 --tmp-dir 一定要指定一个足够大的临时空间，默认可能会写到 /tmp 导致磁盘爆满。我习惯在每个项目目录下建一个 tmp_gatk 文件夹。

GVCF 生成流程一览

下面这张图展示了一个样本从输入到 GVCF 输出的核心步骤。你可以在脑子里过一遍：

验证输出文件

跑完之后，我习惯先做一个快速检查：

tabix -f sample.g.vcf.gz
bcftools stats sample.g.vcf.gz | grep "number of records"

看看记录条数是否合理：30X 人类 WGS 大约 5000 万～1 亿条（含 block）。如果只有几百万，可能参数搞错了，比如误用了 -L 限定了很小的区间。

我踩过的坑

坑1：Java 堆内存不足。 我一开始图省事，用 -Xmx8g 跑 60X 的外显子组，结果频繁 GC 崩溃。后来改成 -Xmx16g 再也没出过问题。对于 30X WGS，16GB 是底线，推荐 32GB。

坑2：临时目录磁盘爆满。 有一次跑 200 个样本的批量任务，把 --tmp-dir 写到了 /tmp，结果系统直接卡死。从那以后我每次都会在命令里显式指定 --tmp-dir，并且监控磁盘使用率。

坑3：GATK 版本不一致。 我遇到过团队里有人用 GATK3.8 的 HaplotypeCaller，有人用 GATK4.2，中间换了参数名字，联合分析时才发现 GVCF 格式有细微差异。建议整个项目统一一个 GATK 版本，最好用 4.2 以上的稳定版。

好了，单样本 GVCF 生成的核心内容就是这些。记住，参数调优永远是服务于你的数据特性：测序深度、reads 长度、物种参考质量。没有万能配方，但有合理的默认起点——我上面给出的那个命令模板就是一个不错的起点。