📚 课程目录第3章 VQSR输入文件详解

第3章 VQSR输入文件详解：VCF、训练集、真值集、资源集

做VQSR之前，得先把输入文件搞明白。很多新手一上来就怼参数，结果模型跑出来乱七八糟——其实八成是输入文件没整对。我做了几年变异过滤，见过最坑的坑，基本都出在这几个文件上。

3.1 核心输入：VCF文件

VQSR的起点，就是你那个合并后的VCF。所有的变异位点、基因型、质量分数，都从这儿来。个人习惯，在跑VQSR之前，先确认VCF里至少包含以下几个关键字段：

QUAL — 位点质量分数，模型用它做主要特征
INFO/DP — 总深度
INFO/QD — 质量与深度比值
INFO/FS — Fisher精确检验的链偏
INFO/SOR — 链偏的比值
INFO/MQ — RMS映射质量
INFO/MQRankSum — 映射质量的秩和检验
INFO/ReadPosRankSum — 读取位置的秩和检验

注意： 我曾经遇到过一个案例，VCF里没有QD和FS，VQSR直接报错。后来发现是GATK版本不对，旧版本不生成某些注解。升级到4.x就解决了。建议跑VQSR之前先用 bcftools view -h 检查一下VCF头部的INFO行。

3.2 训练集 (Training Set)

训练集是VQSR用来「学习」变异好坏的参考。说白了就是一堆已知靠谱的变异位点。模型拿这些位点的特征向量去拟合一个高斯混合模型，然后给每个位点打分。

常见的训练集来源：

数据集	描述	我常用的选择
HapMap 3.3	广泛的人类变异参考，约1.6M个常见SNP	必选，尤其是SNP模型
1000G 高置信度	千人基因组高可信位点	SNP和INDEL都可用
dbSNP	已知变异库，但质量参差不齐	我一般不单独用，混合用更稳
Omni 2.5M 芯片数据	高密度芯片验证过的位点	对SNP模型特别好使

一个实战经验： 如果你处理的是非人物种（比如细菌、植物），没有现成的高质量训练集。我做过一个玉米项目，只能用模拟数据 + 一小部分Sanger测序验证的位点当训练集。结论是：训练集质量比数量重要。宁缺毋滥。

3.3 真值集 (Truth Set)

真值集和训练集容易搞混。其实区别很简单：训练集参与模型训练，真值集不参与。真值集只用来在模型训练完成后，计算敏感性和特异性，从而寻找最优的VQSLOD阈值。

GATK官方推荐的真值集：

SNP真值集：HapMap 3.3（和训练集可以是同一份数据，但用法不同）
INDEL真值集：Mills and Devine数据

嗯，这里有个小陷阱：真值集必须有明确的真/假标签。GATK通过VCF的FILTER字段来判断——PASS代表真阳性，FAIL代表假阳性。我曾经见过一个团队把整份dbSNP当真值集用，结果模型评估出来的阈值完全离谱，因为dbSNP里有很多低质量位点没标FAIL。

3.4 资源集 (Resource Set)

资源集不参与训练，也不参与阈值计算。它的作用是告诉VQSR：「这些位点是被外部数据库收录过的，优先保留」。说白了就是给模型一个先验信息。

最常用的资源集：

dbSNP — 收录过的位点模型会给更高权重
1000G — 尤其适合低频变异

一句话总结： 训练集教模型怎么打分，真值集教模型怎么切分，资源集教模型不要乱杀已知位点。

3.5 输入文件的对齐与准备

这几个文件必须和你的VCF使用相同的参考基因组版本。我吃过这个亏，有一次用了hg19的参考，训练集混进了hg38的，结果模型收敛到奇怪的方向去了。检查版本很简单：

# 查看参考基因组版本
bcftools view -h your_input.vcf | grep "reference"
# 或者看contig行
bcftools view -h your_input.vcf | grep "##contig"

另外，VQSR要求所有输入文件必须经过与VCF相同的预处理流程：

排序 — 染色体坐标一致
索引 — 必须要有.tbi或.csi索引
部分重叠 — 训练集/真值集/资源集的位点不必全覆盖VCF，但重叠部分越大模型越稳

避坑提醒： 我曾经遇到一个案例，VCF中有多等位位点，但训练集里全是双等位位点。VQSR处理时会跳过不匹配的位点，导致有效训练集大小骤减。建议先用 bcftools norm -m -any 分解多等位位点。

3.6 四类输入的关系框架

3.7 实战配置文件示例

下面是一个我常用的GATK VQSR命令模板，帮你理解这四个文件到底怎么用。注意看 --resource 的 known=false/true 和 training=true/false 组合：

gatk VariantRecalibrator \
   -R reference.fasta \
   -V merged_snps.vcf.gz \
   --resource:hapmap,known=false,training=true,truth=true,prior=15.0 \
       hapmap_3.3.vcf.gz \
   --resource:omni,known=false,training=true,truth=false,prior=12.0 \
       1000G_omni2.5.vcf.gz \
   --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 \
       dbsnp_151.vcf.gz \
   --resource:mills,known=false,training=true,truth=true,prior=12.0 \
       Mills_and_1000G_gold_standard.indels.vcf.gz \
   -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
   -mode SNP \
   -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 \
   -O recal_snps.recal \
   --tranches-file recal_snps.tranches \
   --rscript-file recal_snps.R

解读参数：

hapmap 被设为训练+真值，优先值15，意味着模型会拿它当黄金标准。
omni 只训练不真值，因为它虽然是高质量芯片数据，但覆盖位点不如HapMap全。
dbsnp 只当资源集，不训练不真值，优先值2，纯粹用来保位点。
mills 是INDEL的真值集，这里在SNP模式下其实用不到，但放在 -mode SNP 中会被自动忽略。

看到这里你应该明白了：VQSR不是随便丢一堆文件就能跑。每个输入文件扮演什么角色，你得心里有数。我记得有一次帮同事排查问题，他死活调不出好的敏感度，结果发现训练集和真值集用了同一个文件的同一列，导致模型过拟合——真值集不该参与训练的，他却把 training=true 也加上了。

嗯，这章的内容差不多这样。下一件事就是把这些文件准备好，然后开始调参了。记录好自己的每个文件来源和版本号，后期复盘会省很多时间。