变异检测与 VCF 处理

引言

变异检测（variant calling）是基因组学分析的核心产出：把比对好的 BAM 转换成「样本携带哪些变异」的清单（VCF）。这是从「序列数据」到「生物学结论」的关键一跃，也是临床与研究应用的价值所在。相比上游的质控和比对，变异检测的「正确性」更难验证——因为真正的变异清单往往没有金标准。

难点在于「区分真实变异与噪声」。测序错误、比对错误、PCR 错误都会产生「看起来像变异」的信号。一个位点有 3 条读段支持变异、97 条支持参考，这到底是低丰度的真实变异（如肿瘤亚克隆），还是测序错误？判据必须综合碱基质量、比对质量、链偏倚、位置上下文——这正是变异检测算法要解决的核心问题。

工程上的挑战有三：一是工具的复杂度，GATK 的最佳实践有十几个步骤，每步都有参数；二是 VCF 的复杂性，一个 VCF 文件包含 INFO、FORMAT、FILTER 等多层信息，字段含义需要逐一理解；三是场景的多样性，胚系、体细胞、结构变异、拷贝数变异的检测方法完全不同，不能用一套流程。

本文按「变异类型 → GATK 流程 → HaplotypeCaller 原理 → GVCF 与联合分析 → VCF 解析 → 过滤 → 注释 → 体细胞 → 结构变异」的顺序展开。命令基于 GATK 4.5 与 bcftools 1.19。读完你应该能跑通一条变异检测流程，并理解每一步为什么必要。

目录

1. 变异的类型与检测目标
2. GATK 最佳实践流程
3. HaplotypeCaller 原理
4. GVCF 与联合基因分型
5. VCF 格式逐字段解析
6. 变异过滤：硬过滤与 VQSR
7. 变异注释：VEP 与 snpEff
8. 体细胞变异检测：Mutect2
9. 结构变异与 CNV 检测

1. 变异的类型与检测目标

「变异」是一个统称，不同尺度的变异需要完全不同的检测方法：

按「来源」又可分为：

• 胚系变异（germline）：来自父母、存在于所有细胞，杂合或纯合；
• 体细胞变异（somatic）：肿瘤中后天获得，只在部分细胞，丰度可变（VAF）；
• 新生变异（de novo）：父母没有、子代新出现的变异，罕见但重要。

这个分类决定了流程设计：胚系检测用「单样本或联合基因分型」，关注「是否存在」；体细胞检测必须「肿瘤 + 正常配对」，关注「是否仅在肿瘤中」，且要处理低丰度（VAF 5% 甚至更低）的挑战。

一个关键认知：变异检测的敏感性（sensitivity）与特异性（specificity）是权衡的。放宽阈值能检出更多真变异，但也引入更多假阳性；收紧则相反。选哪个取决于应用——临床诊断宁可漏检少（假阳性代价高），群体筛查则要保证不漏掉罕见变异。没有「绝对正确」的参数。

2. GATK 最佳实践流程

GATK（Genome Analysis Toolkit）是变异检测的事实标准，其「最佳实践」（Best Practices）流程是行业参考：

# 前置：BAM 已经过 MarkDuplicates + BQSR
# 1. 单样本变异检测，输出 GVCF（含所有位点，包括非变异）
gatk HaplotypeCaller \
  -R ref.fa -I sample.bam -O sample.g.vcf.gz \
  -ERC GVCF

# 2. 多样本联合（可选，提高准确性）
gatk GenomicsDBImport \
  -V s1.g.vcf.gz -V s2.g.vcf.gz --genomicsdb-workspace-path db \
  -L intervals.bed
gatk GenotypeGVCFs -R ref.fa -V gendb://db -O cohort.vcf.gz

# 3. 过滤（硬过滤或 VQSR）
gatk VariantFiltration -R ref.fa -V cohort.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \
  --filter-expression "FS > 60.0" --filter-name "HighFS" \
  -O filtered.vcf.gz

这套流程的关键设计：

• BQSR 前置：碱基质量重校准让质量值更准确，直接影响变异判定的可信度；
• GVCF 中间格式：HaplotypeCaller 先输出「基因组 VCF」，记录每个位点（含非变异），供多样本联合分析；
• 联合基因分型：多样本一起分析，共享变异位点信息，比单样本检测更准（尤其对低覆盖样本）。

3. HaplotypeCaller 原理

HaplotypeCaller 的核心创新是「局部从头组装」（local de novo assembly），这与传统的「逐位点比较」完全不同：

传统方法（如 bcftools mpileup）：
  逐位点看读段碱基 → 与参考比较 → 统计差异
  缺点：无法区分「两个相邻 indel」和「一个复杂替换」

HaplotypeCaller：
  1. 找「活跃区域」（有变异信号的区域）
  2. 在该区域丢弃现有比对，重新组装读段 → 得到候选单倍型
  3. 用 Smith-Waterman 把单倍型比回参考 → 得到变异
  4. 计算样本携带哪几种单倍型的似然（pair-HMM）

这个「丢弃比对再重组装」的策略对 indel 尤其有效。例如参考是 ACGT，读段显示 A--T（两个碱基缺失）和 A-GT（插入），传统方法会纠结是缺失还是插入，而组装法能直接构造出正确单倍型。

关键参数：

--pcr-indel-model CONSERVATIVE 值得注意：如果建库经过 PCR，PCR 会在 indel 附近引入系统性错误，这个参数能降低假阳性 indel。但对于 PCR-free 建库，设成默认（NONE）即可。

4. GVCF 与联合基因分型

GVCF（Genomic VCF）是 GATK 的一个巧妙设计：它不只记录「检出的变异」，而是记录所有位点的信息，包括「这里没有变异」的证据。

普通 VCF：只列出变异位点
  chr1  100  .  A  G  ...   ← 有变异

GVCF：列出变异位点 + 非变异区块
  chr1  100  .  A  G  ...         ← 变异
  chr1  200  .  A  <NON_REF>  ...  ← 非变异块（"这里和参考一样"）

这个设计的价值在于联合基因分型（joint genotyping）：把多个样本的 GVCF 合并，GATK 能看到「样本 A 在某个位点是变异，样本 B 在同一个位点是参考」，从而在同一位点上对全部样本统一基因分型。这比「每个样本单独检测再合并」准确得多，尤其对低覆盖样本。

# 多样本联合分析
ls *.g.vcf.gz > samples.list
gatk GenomicsDBImport --sample-name-map sample_map.txt \
  --genomicsdb-workspace-path cohort_db -L intervals.bed
gatk GenotypeGVCFs -R ref.fa -V gendb://cohort_db -O cohort.vcf.gz

GVCF 的代价是文件大（记录所有位点，约是普通 VCF 的 3-5 倍）和 GenomicsDBImport 较慢。但对于群体研究（数十到数千样本），联合分析带来的准确性提升是值得的。单样本临床检测可以跳过联合，直接对 GVCF 做 GenotypeGVCFs。

5. VCF 格式逐字段解析

VCF（Variant Call Format）是变异的标准格式。结构分三部分：头部、列头、记录行。

##fileformat=VCFv4.2
##contig=<ID=chr1,length=248956422>
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS  ID  REF  ALT  QUAL  FILTER  INFO  FORMAT  sample1
chr1    100  .   A    G    50    PASS    DP=30;AF=0.5  GT:DP:GQ  0/1:30:99

基因型（GT） 是最关键的字段：0/1 表示杂合（一条参考一条变异），1/1 表示纯合变异，0/0 表示纯合参考，./. 表示缺失。多等位位点如 1/2 表示携带两种不同的变异等位基因。

INFO 里的关键字段：

常用 bcftools 操作：

bcftools view -v snps cohort.vcf.gz -Oz -o snps.vcf.gz    # 只看 SNP
bcftools view -f PASS cohort.vcf.gz                        # 只看 PASS
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\n' cohort.vcf.gz   # 提取字段
bcftools stats cohort.vcf.gz > stats.txt                   # 统计
bcftools isec -p dir a.vcf.gz b.vcf.gz                     # 比较两个 VCF

6. 变异过滤：硬过滤与 VQSR

原始 VCF 里混杂着真变异与假阳性，必须过滤。两种策略：

硬过滤（hard filtering）：对每个位点独立应用阈值。简单、快速，适合小样本量（< 30 样本）。

gatk VariantFiltration -R ref.fa -V raw.vcf.gz \
  --filter-expression "QD < 2.0"       --filter-name "LowQD" \
  --filter-expression "FS > 60.0"      --filter-name "HighFS" \
  --filter-expression "MQ < 40.0"      --filter-name "LowMQ" \
  --filter-expression "SOR > 3.0"      --filter-name "HighSOR" \
  --filter-expression "MQRankSum < -12.5" --filter-name "LowMQRankSum" \
  --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowReadPosRankSum" \
  -O filtered.vcf.gz

VQSR（Variant Quality Score Recalibration）：用机器学习建模，需要大量已知变异位点（dbSNP、HapMap）作为「真集」和「假集」训练。适合大样本量（> 30 样本），准确性更高但需要足够的训练数据。

gatk VariantRecalibrator -R ref.fa -V cohort.vcf.gz \
  --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
  --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
  --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G.vcf.gz \
  --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.vcf.gz \
  -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
  -mode SNP -O snp.recal --tranches-file snp.tranches
gatk ApplyVQSR -R ref.fa -V cohort.vcf.gz --recal-file snp.recal \
  --tranches-file snp.tranches --truth-sensitivity-filter-level 99.7 -O snp.vqsr.vcf.gz

关键取舍：样本少用硬过滤，样本多用 VQSR。VQSR 需要至少几十个样本才能训练出可靠的模型，样本太少会「过拟合」到训练集，反而更差。此外 VQSR 依赖已知位点资源（dbSNP 等），对非人类物种可能不适用。

7. 变异注释：VEP 与 snpEff

变异检测给出「哪里变了」，注释（annotation）给出「变了有什么影响」。主流工具是 VEP 和 snpEff：

# VEP 注释
vep -i cohort.vcf.gz -o annotated.vcf \
  --cache --assembly GRCh38 --offline \
  --vcf --symbol --terms SO --plugin LoF \
  --fields "SYMBOL,Consequence,HGVSc,HGVSp,gnomAD_AF,CLIN_SIG"

# snpEff 注释
snpEff -v GRCh38.105 cohort.vcf.gz > annotated.vcf

注释的核心输出是「后果（consequence）」，按严重程度分级：

除了功能后果，注释还整合了人群频率（gnomAD）、致病性（ClinVar）、保守性（phyloP）等数据库。过滤罕见变异是临床解读的常用策略：一个等位基因频率 > 1% 的变异，通常不是罕见病的致病变异（因为太常见）。

注释的工程要点：VEP 的缓存（cache）需要与参考版本严格匹配（GRCh38 用 GRCh38 缓存）；--offline 模式用本地缓存，避免网络依赖；大批量注释可以并行（VEP 支持 --fork）。

8. 体细胞变异检测：Mutect2

肿瘤体细胞变异检测与胚系检测有本质区别：必须区分「肿瘤特有的变异」与「个体固有的胚系变异」。GATK 的方案是 Mutect2，需要肿瘤 + 正常配对样本：

gatk Mutect2 -R ref.fa \
  -I tumor.bam -tumor tumor_sample \
  -I normal.bam -normal normal_sample \
  --germline-resource af-only-gnomad.vcf.gz \
  --panel-of-normals pon.vcf.gz \
  -O somatic.vcf.gz

# 过滤假阳性
gatk FilterMutectCalls -R ref.fa -V somatic.vcf.gz -O filtered.vcf.gz

Mutect2 的三个关键资源：

• 正常样本：用于区分胚系变异（正常样本也有的变异）；
• gnomAD（germline resource）：人群频率库，用于识别常见胚系多态；
• Panel of Normals（PoN）：一组正常样本的「技术噪声库」，用于识别反复出现的假阳性位点（如比对错误、测序artifact）。

体细胞检测的核心难点是低丰度变异（low VAF）。一个 VAF 5% 的变异，在 100x 覆盖下只有 5 条读段支持，很容易被当成测序错误。Mutect2 用「局部组装 + 贝叶斯模型」来评估，但仍需足够的覆盖深度。临床肿瘤检测常用 500x 甚至 1000x 的深度来捕捉低丰度克隆。

另一个概念是 VAF（Variant Allele Frequency）：变异读段数 / 总读段数。它近似反映肿瘤纯度 × 克隆比例。但要注意 VAF 受肿瘤纯度影响——纯度 30% 的样本里，一个克隆性杂合变异（理论上 VAF 50%）实际 VAF 可能只有 15%。解读 VAF 必须结合肿瘤纯度。

9. 结构变异与 CNV 检测

大于 50 bp 的变异（结构变异，SV）和拷贝数变异（CNV）用完全不同的方法：

结构变异（SV）：检测缺失、插入、倒位、易位。短读方法依赖「不一致的比对信号」：

# Manta：基于 split-read 与 discordant pair
configManta.py --normalBam normal.bam --tumorBam tumor.bam \
  --referenceFasta ref.fa --runDir manta_run
python manta_run/runWorkflow.py -j 16

# Sniffles：长读 SV 检测（更适合）
sniffles -i hifi.bam -v sv.vcf --reference ref.fa

SV 检测的信号类型：

拷贝数变异（CNV）：检测大片段的拷贝数变化（缺失、重复）。基于「覆盖度」而非「序列差异」：

# CNVkit：基于覆盖度与 off-target 归一化
cnvkit.py batch tumor.bam --normal normal.bam \
  --targets targets.bed --fasta ref.fa --output-dir cnv_out

# GATK gCNV：贝叶斯模型
gatk CollectReadCounts -I sample.bam -L intervals.bed -O counts.hdf5

CNV 检测的关键是「归一化」：测序覆盖度受 GC 含量、重复序列、批次影响，必须先校正这些系统性偏差，才能看出真实的拷贝数变化。CNVkit 用「off-target 区域的覆盖度」作为基线来校正，是较稳健的做法。

一个工程建议：SV 与 CNV 检测对读长高度敏感。短读（150 bp）难以精确确定 SV 断点，且对小于 1 kb 的 SV 不敏感。有条件时应补充长读测序，或用「短读 + 长读」联合策略。参见 测序原理 了解长读的优势。

权衡取舍

常见坑清单

1. BAM 缺读组（RG）：GATK 报错或样本混淆；比对时必设 -R，SM 字段正确。
2. 样本少却用 VQSR：训练不足导致过拟合，结果更差；< 30 样本用硬过滤。
3. 参考版本与注释资源不匹配：GRCh37 的 dbSNP 用到 GRCh38 数据上；全流程统一版本。
4. 忽略 QD/FS 过滤：假阳性直接进入结果；用标准阈值过滤（QD < 2、FS > 60）。
5. 把胚系当体细胞：没做配对对照，把所有变异都当体细胞；体细胞检测必须配对。
6. VAF 解读忽略肿瘤纯度：把低 VAF 当低丰度克隆，实际是纯度低；结合纯度解读。
7. GVCF 未做联合分析：多样本各自检测再合并，丢失联合分型的准确性；用 GenomicsDBImport。
8. bcftools 未做 norm：多等位位点未拆分，下游注释出错；用 bcftools norm -m -both。
9. 注释数据库过时：ClinVar/gnomAD 版本陈旧，漏掉最新致病性证据；定期更新。
10. SV 用短读硬扛：小 SV 漏检、断点不准；关键 SV 用长读验证。

小结

变异检测是基因组学的价值出口，它的准确性直接决定下游解读的可靠性。HaplotypeCaller 的局部组装、GVCF 的联合分型、VQSR 的机器学习校准，都是「在噪声中找信号」的工程智慧。理解这些方法的假设与适用边界，比记住命令更重要。

工程上，变异检测最该建立的认知是「场景决定方法」：胚系、体细胞、结构变异、拷贝数变异，各自有专属的流程与参数。用错方法不会报错，但结果不可信——这是最危险的失败模式。另一个要点是「版本一致性」：参考、dbSNP、注释数据库必须同源，任何一处不匹配都会导致坐标错乱。

下一步可以看 RNA-seq 转录组分析 了解表达层面的分析，或 生信流程编排 了解如何把这条复杂流程工程化。如果你的变异结果异常多或异常少，先用 bcftools stats 看看分布，再回头检查 SAM/BAM 与 samtools 那一环的过滤是否合适。