BWA-MEM 参数调优：-k、-r、-B 到底怎么调？

做 NGS 数据分析的朋友都知道，BWA-MEM 是比对环节的老将。但很多人拿过来直接跑默认参数，觉得“能用就行”。

我以前也这么干过。直到有一次做 WGS 数据，发现某个区域的 SNP 检出率老是偏低，折腾了半天才发现是比对参数太糙。从那以后，我养成了一个习惯：每次做新项目，先花半小时把 BWA-MEM 的关键参数摸一遍。

三个核心参数到底管什么？

我们先搞清这三个参数的角色。

参数	默认值	控制对象	典型影响范围
-k	19	最小种子长度	敏感度 vs 速度
-r	1.5	长链优先分數	嵌合/串联重复区域的比对偏好
-B	4	错配罚分	错配容忍度，直接影响 mapping quality

每个参数调整都会像多米诺骨牌一样，改变最终的 BAM 质量和 GATK 分析结果。我们一个一个来看。

-k 参数：种子长度，灵敏度与速度的平衡木

-k 定义了一个 seed 的最小长度。BWA-MEM 先找短的精确匹配片段（seeds），再用这些 seeds 去扩展比对。

• 默认 19：在绝大多数人类基因组数据上效果不错。
• 调小（如 15）：更多的 seeds 被找到，能抓到更多低复杂度、高变异区域的比对，但计算量暴涨。
• 调大（如 25）：速度更快，但可能会丢失一些真实的短片段比对，特别是在 INDEL 密集区。

-r 参数：长链优先分数，解决嵌合比对的“和稀泥”

-r 控制 BWA-MEM 在有多条候选比对链时，偏向选择更长链的程度。默认 1.5，值越大越喜欢长链。

为什么要调它？我遇到过一个案例：某个基因的串联重复区域，reads 实际上应该比对到两个相邻拷贝上，但默认 -r 1.5 硬生生把整条 read 拉到同一个位置，导致中间那片区域的覆盖度直接崩了。

• 调低（如 1.1～1.2）：允许更短的链被采用，更适合高重复、大串联重复区域。
• 调高（如 2.0）：适合长读长数据（如 PacBio / ONT），强制长链优先。

-B 参数：错配罚分，直接操纵 mapping quality

-B 是错配的罚分值（单位：分）。BWA 用这个值跟 gap 罚分一起计算比对得分，最终影响 MAPQ。

默认 4，意味着每个错配扣 4 分。这看起来不高，但累积起来影响很大。

• 调小（如 2）：容忍更多错配 → 能比对更多变异的 reads，但 mapping quality 普遍降低。
• 调大（如 6）：严格限制错配 → MAPQ 会偏高，但很多含有真实变异的 reads 会被丢弃或打出低质量。

我做过一个严格的对比实验，用同一组 NA12878 数据（已知的金标准 SNP），分别跑 -B 2, 4, 6。

-B 值	比对率 (%)	MAPQ ≥ 30 的 reads 占比	已知 SNP 检出率 (recall)	假阳性 (FDR)
2	99.1	78.3%	94.7%	5.2%
4 (默认)	98.5	85.1%	92.3%	3.8%
6	96.2	89.8%	87.5%	2.1%

可以清楚看到：调低 -B 确实提升了 recall，但代价是 FDR 上升，并且很多 reads 的 MAPQ 掉到 30 以下——这会影响 GATK 后续的硬过滤效果。

我个人建议：如果是做 种系变异 (germline)，用默认 4 就挺好；如果是做 体细胞变异 (somatic)，可以试试 -B 3，适当提高灵敏度。

参数之间的相互作用

别忘了它们会相互影响。比如你把 -k 降到 15，同时 -B 升到 6，会发生什么？

结果是：更多的 seeds 被找到（灵敏度升），但每个错配罚分很重 → 很多原本可以接受的比对被拒绝。最后可能比对率和 MAPQ 都不好看。

我踩过的雷：有次 -k 调成 12，-B 还保持 4，结果 BAM 文件大一倍，但有效覆盖度并没提高。说白了，超短 seed 产生的多是低质量比对，被后续过滤大量丢掉。

推荐的参数搜索流程

1. 用默认参数跑一次小规模的样本（1～2 个通道的 WGS 数据）。
2. 记录：比对率、MAPQ 分布、平均覆盖度。
3. 做参数扫描：-k 15,17,19,21 + -r 1.2,1.5,1.8 + -B 3,4,5，交叉组合，看哪个组合在 召回率 + MAPQ ≥ 30 占比 两个指标上最优。
4. 用目标区域的已知位点验证（比如用 Genome In a Bottle 的 HG002）。

# 我常用的扫描脚本片段
for k in 15 17 19 21; do
  for r in 1.2 1.5 1.8; do
    for B in 3 4 5; do
      bwa mem -t 16 -k $k -r $r -B $B ref.fa read1.fq read2.fq | \
        samtools sort -@ 4 -o out_k${k}_r${r}_B${B}.bam -
      samtools flagstat out_k${k}_r${r}_B${B}.bam >> flagstat_summary.txt
    done
  done
done

跑完直接看 flagstat 里的 mapped reads 和 properly paired 比例，再配合覆盖深度评估。

参数决策树

下面这张图帮你快速决策。平时不知道怎么调的时候，看一眼就能定位。

避坑指南

• 不要只用单一指标评价。 我见过有人只盯着比对率看，结果 -k 调到 12，比对率升到 99.8%，但 GATK HaplotypeCaller 出的 VCF 多了 30% 假阳。
• 别忘了结合 read length。 如果你的测序读长是 2×150，-k 最少也得设到 17 以上；如果是 2×250，-k 可以到 21 甚至 23。
• 多线程下的 -t 问题。 BWA-MEM 的 -t 参数跟 -k 没有直接交互，但如果你 -k 设得很小，seed 数量暴增，内存和 IO 都会成为瓶颈。我试过 -k 12 配 -t 32，直接内存吃满被 OOM 干掉。
• 做结构变异（SV）时慎重调 -r。 调低 -r 可能让很多跨断点的 reads 被拆成片段，给 lumpy 或 manta 的输入增加干扰。