BQSR概述与重要性

做 NGS 数据分析的同行应该都有同感——碱基质量值这玩意儿，看似只是一个数字，可它直接决定了你后续找变异的可信度。我记得刚开始接触 GATK 那会儿，老板丢给我一批全外显子数据，我跳过 BQSR 直接跑了 HaplotypeCaller，结果出来的候选位点哗啦啦一堆假阳性。后来一查，原来是测序仪给的原始质量分跟实际错误率有偏差。嗯，BQSR 要解决的就是这个偏差。

什么是碱基质量重校准

说白了，BQSR（Base Quality Score Recalibration）是一种数据驱动的校正方法。它拿你现有的 BAM 文件，结合一个已知的、经过验证的高质量变异集合（比如 dbSNP、1000 Genomes 的可靠位点），去统计测序仪报告的原始碱基质量值与真实错误率之间的差异，然后训练一个模型，最后用这个模型重新计算每个碱基的质量得分。

你看下面这个流程图，就清楚它在整个 GATK 最佳实践里站什么位置了。

你看，BQSR 处在标记重复之后、HaplotypeCaller 变异检测之前。这一步做不好，后面所谓的"高可信位点"可能就是空中楼阁了。

为什么要做 BQSR

测序仪在吐出碱基时，会附带一个 Phred 质量得分（通常是简称 Q 值）。公式是 Q = -10 × log₁₀(P_error)，比如 Q30 意味着该碱基出错的概率是 1/1000。但问题是——

测序仪给的质量分，往往不是真实错误率的完美体现。

原因有几方面：

• 测序化学本身的偏好 — 比如某些碱基组成（GC 含量极端区域）会导致荧光信号串扰，读出的质量分虚高。
• 循环位置效应 — 一条 read 的开头和末尾错误率通常不同，仪器可能对 read 中段的质量分估计不准。
• 测序仪器与试剂批间差异 — 不同 flow cell、不同 lane、不同 run 的误差模式不一样，但软件给的默认质量转换参数是固定的。
• 原始碱基上下文 — 前面一个碱基是什么、后面一个碱基是什么，都会影响当前碱基的错误率。举个例子，我在一个肿瘤 panel 项目里遇到过：在 GGC 三联体区域，测序仪总是把 G 错读成 A，但错误率比 Q 值指示的高了将近 10 倍。

BQSR 干的事情就是：用一个"已知可信"的变异集作为金标准，把所有这些系统性偏差找出来，然后修正每个碱基的质量值，让它更贴合你手上这批数据的真实错误率。

BQSR 在 GATK 流程中的位置

GATK 的最佳实践流程（我用的版本是 4.x 系列）大致分以下几个阶段：

1. 数据准备阶段：原始 FASTQ → BWA 比对 → 排序 → 标记重复（MarkDuplicates）
2. 质量校正确认阶段：使用 BQSR 工具对 BAM 进行基水平质量重校准
   • 先是 BaseRecalibrator 生成重校准表（.recal 文件）
   • 然后用 ApplyBQSR 将重校准表应用到原始 BAM 上，输出校准后的 BAM
3. 变异检测阶段：HaplotypeCaller、Mutect2 等工具基于校准后的质量值进行 calling
4. 过滤与注释阶段：VQSR、硬过滤、Funcotator 注释

你给我留意，某些教程会偷懒不跑 BQSR。我的建议是——除非你百分之百确定你的数据来自一个经过严格质控的、且运行参数完全标准的平台（比如某些医院内部的 Illumina 自动化管道），否则一定不要跳过去。

跳过 BQSR 前后效果对比（一个真实 WES 样本）

指标	跳过 BQSR	启用 BQSR
原始 SNP 数（硬过滤前）	145,203	138,470
dbSNP 已知位点占比	82.7%	89.2%
Ti/Tv 比值	2.01	2.19
假阳性候选位点（经验证）	~5.2%	~1.8%

这个表里的 Ti/Tv 比值，如果是人全外显子数据，正常范围应该在 2.0~2.2 之间。跳过 BQSR 的那一列明显偏低，说明很多真正的变异被假阳性污染，把比值拉下来了。我自己做项目时，只要看到 Ti/Tv 小于 2.0，第一反应就是回去检查 BQSR 有没有跑对。

BQSR 的前提条件

跑 BQSR 前，你需要准备这么几个东西：

• 已标记重复的 BAM 文件 — 这一步是必须的，因为在标记重复之前做重校准意义不大（重复片段会干扰质量值统计）。
• 参考基因组 FASTA — 必须与比对时用的是同一个版本，比如 GRCh38/hg38。
• 已知变异点资源集 — 至少需要一份可靠的、经过严格过滤的 SNP/Indel 集合。常用的有：
  • dbSNP（经过 VQSR 过滤的 build 155 或更晚版本）
  • 1000 Genomes 高可信位点（Mills & 1000G gold standard）
  • 如果做肿瘤样本，最好配合 gnomAD 的常见变异。

BQSR 是怎么工作的（简要原理）

我不打算在这里堆数学公式，你只需知道核心是两个步骤：

1. 统计偏差：遍历 BAM 文件中所有碱基，针对每个碱基记录其特征（原始质量分、read 的位置、上下游碱基、所在的小区段上下文），然后跟已知可信集合比对，算出每个特征组合下的经验错误率。这些错误率跟仪器报告的 Phred 质量分之间的差异，就是我们要校正的量。
2. 应用校正：根据上一步统计出来的偏差表，对 BAM 中每个碱基的原始质量分做加减调整，得到新的质量分。新质量分写入 BAM 的 QUAL 字段，旧分会被保留在 OQ 标签里方便溯源。

简单说，就是让 "Q 值 = -10 × log₁₀(实际错误率)" 这条等式不再被破坏。