📚 课程目录第1章 HaplotypeCaller概述

1. HaplotypeCaller概述

说实话，我刚接触 GATK 那会儿，对 HaplotypeCaller 这名字有点怵。
听上去像是个很高深的算法。用了一段时间才发现，它的核心思路其实很直白：
对着测序数据，把可能存在的变异位点一个一个揪出来，然后算概率。

我个人习惯把 HaplotypeCaller 看作整个 GATK 流程里的“发动机”。
前面 BWA 比对、MarkDuplicates 去重，都是为了给它喂干净数据。
它输出的 GVCF 文件，质量直接决定了后续联合分型的结果靠不靠谱。

  一句话总结： HaplotypeCaller 是 GATK Best Practices 中用于种系变异检测的核心工具。它不直接硬比对，而是在局部区域重新组装单倍型，再用 Bayesian 模型计算每个位点的基因型概率。

为什么叫“单倍型”？

你想想看，传统的 SNP calling 工具，比如 Samtools 的 mpileup，是逐碱基去看的。
但 HaplotypeCaller 不是。它先把一个区域里的 reads 重新组装成短序列片段。
这些片段，就是局部的单倍型。然后它把这几种单倍型跟参考基因组比对，找差异。

我在项目里遇到过一位同事，他坚持认为 GATK 的 HaplotypeCaller 是个“黑盒子”。
其实不是。它做的就是两件事：

组装：在含有变异的活跃区域，基于 De Bruijn 图构建局部单倍型
评估：把每个样本的 reads 映射到这些单倍型上，用概率模型算出最可能的基因型

说白了，它不是在“找”变异，而是在“算”变异。这个区别很关键。

GATK 流程里的位置

下面这张图，是我项目里一直贴在墙上的。它把 HaplotypeCaller 在整个流程中的角色画清楚了。
（注意：这是单样本阶段的最后一步，也是联合分型的入口。）

原始 FASTQ
    │
    ▼
BWA MEM 比对  →  BAM/SAM
    │
    ▼
MarkDuplicates 去重  →  clean BAM
    │
    ▼
BaseRecalibrator 碱基质量校正  →  recalibrated BAM
    │
    ▼
╔══════════════════════════════════════╗
║    HaplotypeCaller（单样本 GVCF）    ║
╚══════════════════════════════════════╝
    │
    ▼
单个样本 GVCF (每个样本一个)
    │
    ▼
CombineGVCFs / GenomicsDBImport 合并
    │
    ▼
GenotypeGVCFs 联合分型  →  群体 VCF

经验之谈： 我建议你在跑 HaplotypeCaller 之前，先确认 BAM 文件的染色体顺序是 lexicographical（比如 chr1, chr10, chr11...）还是 natural（chr1, chr2, chr3...）。不一致会导致合并时报错。

核心原理：活跃区域与 Bayesian 模型

HaplotypeCaller 是怎么工作的？我拆成几个关键步骤来讲。

定义活跃区域
扫描 BAM，找那些有足够多不一致碱基或 Insertion/Deletion 证据的区域。
这一步会忽略那些所有 reads 都跟参考一致的区域，节省大量计算。
组装局部单倍型
在这个区域里，用 De Bruijn 图算法，把 reads 重新连接成候选单倍型。
为什么要重新组装？
因为原始的 reads 可能跨一个变异长度不够，组起来才能看到完整的变异背景。
确定基因型似然值
把每个样本的 reads 比对到刚才组装出的单倍型上。
用 PairHMM 算法算一个矩阵：每个 read 对应每个单倍型的可能性。
计算 PL 值与 GQ 值
基于 Bayesian 公式，给定观测的 reads 数据，计算三种基因型（HomRef、Het、HomVar）的后验概率。
输出 PL（Phred-scaled likelihood）和 GQ（Genotype Quality）。
输出 GVCF 记录
无论是否发现变异，每个位点都会有一条记录。
纯合参考位点被压缩成 REF 块（非变异块），变异位点则单独列出。

注意： 活跃区域长度默认是 100 bp 左右，但如果你用了 --max-assembly-region-size 调整，要小心内存消耗。我曾经在一个高覆盖度的外显子数据上设了 500 bp，结果一个样本就跑崩了。

步骤	关键参数	我的建议
活跃区域检测	`--active-probability-threshold`	默认 0.002，一般不动
组装	`--max-assembly-region-size`	外显子用 200，全基因组用 100
PairHMM 计算	`--pair-hmm-implementation`	建议用 `VECTOR_LOGLESS_CACHING`，更快
输出质量过滤	`--standard-min-confidence-threshold-for-calling`	默认 30，低深度项目可以降到 20

应用场景：什么数据适合用它

HaplotypeCaller 不是万能的。它主要针对 种系变异，也就是从受精卵就有的变异。
体细胞突变、肿瘤正常配对、单细胞数据，都不推荐直接用这个工具。

我把它常见的应用场景归纳成三类：

全基因组测序（WGS）
覆盖度 30x 以上，HaplotypeCaller 能很好地处理基因组上的重复区域和 Indel。
全外显子组测序（WES）
靶向区域的深度高，组装出来的单倍型更准，但注意染色体间易位区域容易假阳性。
群体遗传学项目
多份样本同时用 HaplotypeCaller + GVCF 模式，再联合分型。
这是目前千人基因组项目级别的主流做法。

应用场景对比表：

场景	推荐模式	注意事项
单人全基因组	单样本 GVCF	内存要够，至少 8G 每个线程
家系分析	GVCF 模式 + 合并	用 GenomicsDBImport 合并更快
大群体分型（>100 人）	GVCF 模式 + GenomicsDB	按染色体拆分，并行跑
特殊区域（MHC、着丝粒）	单样本 GVCF + --max-assembly-region-size 调小	调小防止组装发散

与其他工具的横向对比

除了 GATK，市面上还有 FreeBayes、DeepVariant、Platypus 这些工具。
它们各有侧重。我在选型时一般这么判断：

FreeBayes：支持多等位基因，但 Indel 的处理不如 HaplotypeCaller 精细
DeepVariant：基于深度学习，准确率高，但需要 GPU，不容易做群体联合分型
Platypus：快，但已经多年不维护，社区逐渐放弃

HaplotypeCaller 的优势在于：它的 GVCF 模式可以保持每个样本的原始信号，联合分型时才做最终决策。这个“先算概率、再定基因型”的架构，让它在多样本项目中特别稳。

避坑提醒： 我曾经在 FreeBayes 和 HaplotypeCaller 之间做过对比，发现 FreeBayes 在某些 Indel 富集区域会多报 20% 的假阳性。而 HaplotypeCaller 因为有重新组装这一步，对这些区域的过滤效果明显更好。

本章小节：你学到了什么

HaplotypeCaller 不是逐碱基比对，而是局部组装单倍型，再用 Bayesian 模型定基因型
它在 GATK Best Practices 中扮演从 BAM 到 GVCF 的关键角色
GVCF 模式的设计，天然适合多样本联合分型
应用场景以种系变异为主，WGS 和 WES 都可以，但要注意参数调整