1. 变异准确性评估概述
做生物信息学这些年，我接触过不少变异检测的项目。不管是肿瘤体细胞突变检测，还是遗传病家系分析，大家最常问的一个问题就是：“这堆突变到底准不准？”
说白了，测序仪吐出来的VCF文件里，有的位点是真的，有的只是测序噪音或者比对错误。如果你不评估准确性，直接拿着这些结果去做临床决策，那风险可不小。
我个人习惯在任何一个变异检测流程里，先搭一套准确性评估的框架。为什么要做这件事？因为只有知道假阳性多少、假阴性多少，你才能对结果有信心，也才能去优化流程参数。
什么是变异准确性评估？
准确性评估，就是拿你检测出来的变异集合，跟某个“金标准”或“参考答案”去对比，然后算出各种评价指标。
- 真阳性（TP）：你检出了，且参考集里也有的变异。
- 假阳性（FP）：你检出了，但参考集里没有的变异（也就是你多报了）。
- 真阴性（TN）：你没检出，参考集里也确实没有的变异。
- 假阴性（FN）：你没检出，但参考集里实际存在的变异（你漏报了）。
把这四个数凑齐，就能算出一堆指标：召回率（Sensitivity/Recall）、精确率（Precision）、F1-score，还有我们后续要重点讲的ROC曲线和AUROC。
| 指标 | 公式 | 含义 |
|---|---|---|
| 召回率 (Recall) | TP / (TP + FN) | 真实变异中你找回来了多少 |
| 精确率 (Precision) | TP / (TP + FP) | 你报出的变异里有多少是真的 |
| F1-score | 2 × (P × R) / (P + R) | 召回与精确的调和均值 |
| 假阳性率 (FPR) | FP / (FP + TN) | 误报的比例（ROC曲线横轴） |
为什么需要准确性评估？
有一次我在做肿瘤ctDNA低频突变检测，用的是某个流行的caller。默认参数下跑出来假阳性极高，大概有一半的“突变”其实都是测序错误。如果不做评估，拿着这个结果给医生看，那简直是在害人。
评估的作用，我归纳了四点：
- 判断流程是否可用——比如召回率低于70%的流程要慎重。
- 指导参数调优——通过对比不同阈值的TPR/FPR，找到最佳平衡点。
- 比较不同工具——比如GATK vs. DeepVariant，哪个更靠谱？
- 建立置信区间——在临床报告里可以写“该位点检出概率为95%”。
应用场景一览
变异准确性评估不是实验室里关起门来做的研究，它出现在很多实际场景中。我列几个我碰过的：
| 场景 | 典型的参考标准 | 评估重点 |
|---|---|---|
| 肿瘤体细胞突变检测 | 混合标准品 (如 HorizonDx) | 低频变异 (VAF≤5%) 的召回率 |
| 遗传病家系WES分析 | 千人基因组高置信位点 | SNV精确率 > 99% |
| 病原体宏基因组检测 | 模拟数据 (in silico spike-in) | 物种级假阳性控制 |
| 新流程上线验证 | GIAB (Genome in a Bottle) | 全基因组级别的 Recall/Precision |
你想想看，如果你要上线一个自动化变异检测流水线，不做准确性评估就往外推数据，哪天出了问题追溯起来，你连自己算法水平如何都不清楚，那就尴尬了。
一个真实案例：GIAB基准集评估
GIAB（Genome in a Bottle）提供了NA12878等高置信度变异集。我去年帮一个团队评估他们自研的caller，流程大概是这样的：
# 下载GIAB v4.2.1 高置信度VCF
wget ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/NA12878_HG001/NISTv4.2.1/
# 用hap.py进行比对和评估
hap.py \
GIAB.HG001.vcf.gz \
my_caller_output.vcf.gz \
-f GIAB.HG001.bed \
-o eval_result \
--threads 8跑完之后会输出一个eval_result.summary.csv，里面直接就是TP、FP、FN的计数和指标。嗯，这里要注意：hap.py默认只比较同一基因型的变异，你要是想让它对杂合/纯合不敏感，可以加--no-roc参数。
核心逻辑：从评估到ROC
以我自己的经验，做变异准确性评估不是为了拿到一个静态的数字就完事。更重要的价值在于——通过改变阈值或质量分数，画出ROC曲线，找到最优阈值。比如你调整caller的--min-qual参数，每个阈值都会得到一组TPR和FPR，把这些点连起来，就是ROC曲线。曲线下面积AUROC越大，说明区分真假变异的能力越强。
下面这张图概括了本章的知识脉络：
从金标准参考集出发，比对得到四个核心计数（TP/FP/TN/FN），然后计算召回率、精确率、F1-score，最后通过改变质量阈值绘制出ROC曲线。这就是整个评估的骨架。
如何开始你的第一次评估？
如果你手头有一个VCF文件，还有对应的参考集（比如从GIAB下载的），建议按下面几步走：
- 确认格式：两个VCF的染色体命名、参考基因组版本必须一致。我遇到过hg19和hg38混用的坑。
- 归一化：用
bcftools norm -f ref.fa对两边都做左对齐和降维。 - 运行评估工具：推荐 hap.py 或 RTG Tools vcfeval。我个人偏爱rtg vcfeval，因为它对复杂变异（indel、MNV）处理更聪明。
- 解读输出：重点关注混淆矩阵和 ROC 数据（如果工具生成的话）。
# 使用rtg vcfeval的示例
rtg vcfeval \
-b GIAB_highconf.vcf.gz \
-c my_calls.vcf.gz \
-t SDF/GRCh38.sdf \
-o eval_out \
--sample NA12878,NA12878 \
--bed-regions GIAB_highconf.bed跑完后，eval_out/roc.tsv 里就是每个质量阈值下的TPR和FPR，可以直接喂给Python画ROC曲线。
讲到这里，你应该对变异准确性评估为什么重要、怎么用、有哪些核心指标有个整体印象了。后续我们会一步步深入ROC曲线的数学原理、Python实现、以及如何用评估结果优化你的caller参数。