序列比对：BWA 与比对算法

引言

比对（alignment）是生信流程的中枢：把每条读段定位到参考基因组上的坐标，后续的变异检测、表达定量、覆盖度分析全都建立在这个坐标之上。比对错了，后面全错——这是生信里最「上游」的正确性来源。也正因为如此，比对算法的选择与参数调优值得深入理解。

从算法视角看，比对要解决的是一个「近似字符串匹配」问题：给定一条 150 bp 的读段和一份 3.1 Gbp 的参考基因组，找出读段最可能来自哪个位置，允许少量错配与插入缺失。暴力搜索显然不可行（3.1e9 × 1.5e8 读段的天文数字），必须借助索引结构与启发式搜索。

工程上的挑战有三：一是速度，一个 30x WGS 有上亿条读段，每条都要在秒级内找到位置；二是准确性，尤其是重复区、变异附近、测序错误附近的比对；三是可扩展性，工具要能利用多核并行，还要支持剪接比对（RNA-seq）、长读（三代）等变体场景。

本文按「问题定义 → 索引原理 → BWA 实战 → 打分与 CIGAR → 工具对比 → 后处理 → 评估」的顺序展开。重点讲清「BWA-MEM 为什么这样设计」以及「参数怎么调」，而不是罗列命令。读完你应该能判断一次比对的质量，并在结果异常时定位是数据问题还是参数问题。

目录

1. 比对问题的定义与挑战
2. 索引结构：FM-index 与 BWT
3. BWA 家族：aln、mem 与 bwasw
4. BWA-MEM 算法与参数调优
5. 比对打分、CIGAR 与 MAPQ
6. 替代比对工具对比
7. 比对后处理链
8. 比对率与比对质量评估
9. 特殊场景：剪接、长读与古 DNA

1. 比对问题的定义与挑战

比对的形式化定义：给定读段 R 与参考序列 G，找到 G 中与 R 最相似的子串位置，允许替换（substitution）、插入（insertion）、缺失（deletion）三类编辑操作。相似度用打分函数衡量，通常替换罚分、插入缺失罚分（gap penalty）不同。

难点来自四个方面：

规模。参考基因组 3.1 Gbp，读段数上亿。如果每条读段都要扫描整个参考，复杂度是 O(N×M)，完全不可行。索引技术把「扫描参考」变成「查表」，把复杂度降到接近 O(读段长度)。

重复序列。人类基因组约 50% 是重复序列（转座子、卫星 DNA、片段重复）。一条来自重复区的读段可能在基因组里有几十个「同样好」的位置，工具必须报告「比对不确定」（MAPQ=0）。这不是工具的缺陷，而是生物学事实——短读长本身无法区分。

变异与错误的混淆。读段与参考的差异可能来自真实变异、测序错误或比对错误。好的比对算法应倾向于把差异集中成少数几个变异位点（符合生物学），而非分散成多个错配（符合测序错误随机分布）。这个「简约性偏好」是 BWA-MEM 打分的隐含假设。

速度与准确性的权衡。精确算法（如 Smith-Waterman 动态规划）能保证最优解，但太慢。所有实用工具都用「种子 + 扩展」的启发式：先用快速索引找到候选位置（种子），再在候选位置附近做精细比对。启发式可能错过最优解，但换来了可接受的运行时间。

理解这些难点，才能理解为什么比对没有「完美工具」，只有「适合场景的工具」。

2. 索引结构：FM-index 与 BWT

BWA 的核心是 FM-index（Full-text index in Minute space），它基于 Burrows-Wheeler 变换（BWT） 构建。理解 BWT 是理解 BWA 为何如此快的关键。

BWT 的基本思想：把字符串的所有循环移位排序，取最后一列。这个变换是可逆的，且变换后的字符串高度可压缩（相似字符聚集）。更重要的是，它支持「后向搜索」：从读段的末端往前逐字符匹配，每一步只需 O(1) 的区间查询。

参考序列: banana$
循环移位排序后取末列 → BWT: annb$aa

后向搜索 "an":
  从 '$' 区间开始，逐个字符向前查
  每一步用 LF-mapping 更新区间
  最终区间大小 = 匹配次数

FM-index 在 BWT 基础上加了两样东西：检查点（checkpoint） 与 后缀数组采样，把「区间查询」和「定位到具体坐标」都做到高效。结果是：索引大小约为参考基因组的 1-1.5 倍（人类约 4-5 GB），查询复杂度与参考大小无关（只与读段长度相关）。

工程含义：

• 索引是一次性开销：bwa index 建人类索引需要约 1 小时、数 GB 磁盘。索引建好后可复用，应作为共享资源而非每次重建。
• 内存占用与索引相关：bwa mem 加载索引需要约 3-4 GB 内存（人类），这是每进程的固定开销，多进程并行时要注意总内存。
• BWT 不支持任意编辑：纯 FM-index 只支持精确匹配。要支持错配，BWA 用「多次后向搜索 + 回溯」来枚举允许的错配组合，这是种子阶段的核心。

理解了索引，就能理解 BWA 的运行特征：索引加载慢（一次性），但每条读段的查询极快。这也是为什么「读段数少但参考大」的场景下，索引开销占比高。

3. BWA 家族：aln、mem 与 bwasw

BWA 有三个子命令，对应不同的算法代际：

bwa aln 是初代算法，逐读段做后向搜索，允许有限错配。它的缺点是慢、且不支持长读和 gapped 比对（不能处理插入缺失）。现在只在处理古 DNA 的超短读（< 50 bp）时偶尔使用，因为 mem 对极短读的种子策略不友好。

bwa mem 是当前绝对主力。它的核心创新是 MEM（Maximal Exact Match）种子 + 链式（chaining）：先找出读段与参考的所有「极大精确匹配」，再把这些种子串成链（chain），最后对链做仿射间隙的 Smith-Waterman 扩展。这个策略对长读、含 indel 的读段都表现优异。

# BWA-MEM 标准用法（双端）
bwa mem -t 16 \
  -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA\tLB:lib1' \
  ref.fa R1.fq.gz R2.fq.gz \
  | samtools sort -@ 8 -m 2G -o aln.bam -

samtools index aln.bam

注意 -R 参数：它给所有读段加上**读组（Read Group）**标签。这是极易被忽略但极其重要的一步——GATK 等下游工具依赖 RG 区分样本，缺失 RG 会导致 BQSR 报错或样本混淆。SM（sample）字段尤其关键。

4. BWA-MEM 算法与参数调优

BWA-MEM 的工作流可以拆成四步：

1. 种子查找：对读段每个位置，找极大精确匹配（MEM）
2. 种子过滤：丢弃过短的种子（默认 -k 19），保留有区分度的
3. 链式：把同一对角线上相邻的种子串成链，按得分排序
4. 扩展：对最优链做仿射间隙 Smith-Waterman，得到最终比对与 CIGAR

关键参数：

几个实用调参场景：

高错误率数据（如古 DNA、某些 FFPE 样本）：测序错误多，标准罚分（-B 4）会把真实比对判为低分。可降低 -B 到 3，或增加 -L 容忍更多软剪切。

重复区密集的物种（如玉米）：种子太短会在重复区产生海量候选。可增大 -k（如 25）提高种子特异性，代价是敏感性略降。

比对率低：先别急着调参，用 -x 预设排查。BWA-MEM 提供几个预设：

bwa mem -x intractg   # 种内高变异（如不同菌株）
bwa mem -x pacbio     # PacBio 长读（已不推荐，用 minimap2）
bwa mem -x ont2d      # ONT 2D 读
bwa mem -x ava-ont    # ONT 组装纠错

对于剪接感知的 RNA-seq 比对，BWA-MEM 不是好选择——它不支持跨内含子的「大间隙」比对。应该用 STAR 或 HISAT2。BWA 只适合 DNA 比对。

5. 比对打分、CIGAR 与 MAPQ

比对结果的核心信息有三：位置（坐标）、CIGAR 串（比对详情）、MAPQ（比对质量）。

CIGAR（Compact Idiosyncratic Gapped Alignment Report） 用一串「长度+操作」描述读段与参考的对应关系：

CIGAR: 100M2I50M3D20M5S
  100M  → 100 个碱基匹配/错配（M = alignment match）
  2I    → 2 个碱基插入（读段有、参考无）
  50M   → 50 个匹配
  3D    → 3 个碱基缺失（参考有、读段无）
  20M   → 20 个匹配
  5S    → 5 个碱基软剪切（soft clip，读段末端未比对）

操作符含义：
  M/I/D  → 消耗读段或参考（比对内）
  S      → soft clip（读段保留但不比对）
  H      → hard clip（读段直接截断，不在序列中）
  N      → 跳过分隔（spliced，跨内含子）
  =/X    → 明确的匹配/错配（GATK 的 HaplotypeCaller 输出）

MAPQ（Mapping Quality） 是「这个比对位置正确的置信度」，用 Phred 尺度表示：

MAPQ = 60  → 唯一比对，几乎确定正确
MAPQ = 0   → 多重比对，无法确定位置
MAPQ = 1-59 → 有一定竞争位置，置信度递减

MAPQ 的计算依赖「有多少个同样好的备选位置」。如果一个读段在基因组里有 3 个同样好的位置，MAPQ 会显著降低（甚至为 0）。这是后续过滤的关键依据：变异检测通常只用 MAPQ ≥ 20 的读段，避免重复区产生的假变异。

理解 CIGAR 和 MAPQ 的工程价值：当变异检测出现假阳性时，第一步就是看该位点的读段 CIGAR 与 MAPQ。如果大量读段是软剪切或 MAPQ=0，说明该区域比对不可靠，变异不可信。

6. 替代比对工具对比

BWA 不是唯一选择，不同场景有更合适的工具：

选择逻辑：

• DNA 短读 → BWA-MEM（或 BWA-MEM2 求速度）；
• RNA-seq → STAR（准确、全）或 HISAT2（省内存）；
• 长读 → minimap2（配 -ax map-hifi 等预设）；
• ChIP-seq/ATAC-seq → Bowtie2（内存友好）；
• 高计算量、追求速度 → 考虑 GPU 加速版本，参见 GPU kernel 优化 。

BWA-MEM2 值得特别一提：它在算法上完全兼容 BWA-MEM（输出一致），但用 SIMD 指令重写了关键循环，速度提升 1.5-3 倍。代价是索引更大（约 2 倍）。如果比对是流程瓶颈且有足够内存，切换到 BWA-MEM2 是「零风险」的加速手段。

7. 比对后处理链

比对输出的是 SAM，要变成可用的 BAM 还需要一串后处理。标准链条：

# 1. 比对 + 排序（管道衔接，避免中间文件）
bwa mem -t 16 -R '@RG\tID:s1\tSM:s1\tPL:ILLUMINA' ref.fa R1.fq R2.fq \
  | samtools sort -@ 8 -m 2G -o sorted.bam -

# 2. 标记重复（PCR 重复）
gatk MarkDuplicates -I sorted.bam -O dedup.bam \
  -M metrics.txt --CREATE_INDEX true

# 3. 碱基质量重校准（BQSR）
gatk BaseRecalibrator -I dedup.bam -R ref.fa \
  --known-sites dbsnp.vcf --known-sites indels.vcf \
  -O recal.table
gatk ApplyBQSR -I dedup.bam -R ref.fa --bqsr-recal-file recal.table -O recal.bam

# 4. 索引
samtools index recal.bam

每一步的作用：

• 排序（sort）：按坐标排序，是几乎所有下游工具的前提。排序本身很吃内存（需要在内存中缓冲），-m 控制每线程内存，-@ 控制线程数。用 -T 把临时文件放本地盘。
• 标记重复（MarkDuplicates）：PCR 重复会虚高覆盖度、产生假变异。注意是「标记」（加 flag）而非删除，因为某些重复是真实的。
• BQSR（碱基质量重校准）：这是 GATK 最佳实践的争议点。它通过已知变异位点（dbSNP）学习「系统性的质量值偏差」，重新校准质量值。对于现代测序仪，这个收益已不大，且需要大量已知位点。是否启用取决于流程传统——GATK 最佳实践推荐启用，但很多新流程已省略。
• 索引：生成 .bai，供按区间随机访问。

这套后处理的每一步都产出中间 BAM，是磁盘占用的主要来源。用流程引擎管理它们的生命周期，避免堆积。

8. 比对率与比对质量评估

比对完成后必须评估质量。核心指标：

# 比对统计
samtools flagstat aln.bam
# 输出关键行：
#   mapped: 95.2%        → 比对率
#   properly paired: 92.1% → 正确配对率
#   duplicates: 12.3%     → 重复率

# 更详细的统计
samtools stats aln.bam | grep ^SN | head -20
# 覆盖度（按区间）
mosdepth --by 1000 sample aln.bam

比对率低的常见原因（按概率排序）：

1. 参考基因组不对：用了错误的物种或版本，这是最致命的，比对率会掉到 10% 以下；
2. 接头/污染残留：质控没做干净，未比对读段含大量接头；
3. 物种差异过大：样本与参考的进化距离远（如不同菌株），需要 -x intractg；
4. 测序质量问题：读段太短或错误率太高。

排查方法是提取未比对读段看内容：

samtools view -f 4 aln.bam | head -1000 | awk '{print $10}' | head
# 若看到 AGATCGGAAGAGC 开头 → 接头残留
# 若看到随机序列 → 可能是污染

9. 特殊场景：剪接、长读与古 DNA

剪接比对（RNA-seq）：mRNA 成熟过程中内含子被切除，所以 RNA 读段在基因组上可能跨越两个相距很远的外显子。这需要「剪接感知」的比对，CIGAR 里用 N 表示跳过内含子。BWA 不支持，必须用 STAR 或 HISAT2：

# STAR 剪接比对
STAR --runThreadN 16 --genomeDir star_index \
  --readFilesIn R1.fq.gz R2.fq.gz --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix sample_

长读比对：长读的错误模式（插入缺失为主）与短读不同，BWA 会误判。用 minimap2 的预设：

minimap2 -ax map-hifi -t 16 ref.fa reads.fq.gz | samtools sort -o hifi.bam -
minimap2 -ax map-ont  -t 16 ref.fa reads.fq.gz | samtools sort -o ont.bam -

古 DNA（aDNA）：古 DNA 高度降解，读段极短（平均 50-70 bp）且有特征性的「C→T 脱氨基」损伤。比对时需要特殊处理：用 bwa aln 的宽松参数（-n 0.01 -o 2）而非 mem，并考虑损伤模式（可用 mapDamage 分析）。这是「短读 + 高错误」的极端场景，标准流程不适用。

这三个场景的共同启示是：没有万能比对工具，场景决定工具。理解每种场景的生物学特征（剪接、长读错误模式、古 DNA 损伤），才能选对工具与参数。

权衡取舍

常见坑清单

1. 忘记 -R 读组标签：下游 GATK 报错或样本混淆；比对时必设 -R，SM 字段尤其重要。
2. 参考版本不匹配：用 hg19 索引配 hg38 序列，比对率暴跌；索引与序列必须同源。
3. 比对率低却盲目调参：先排查参考、污染、接头，再考虑参数；多数低比对率不是参数问题。
4. RNA-seq 用 BWA：不支持跨内含子比对，导致大量未比对；必须用 STAR/HISAT2。
5. 长读用 BWA-MEM：错误模式不匹配，假阳性高；用 minimap2 的长读预设。
6. 排序临时文件放共享存储：小文件风暴拖垮 Lustre；用 -T /scratch 重定向到本地盘。
7. 忽略 MAPQ 过滤：重复区低 MAPQ 读段产生假变异；变异检测前过滤 MAPQ < 20。
8. 索引重复重建：每次运行都 bwa index，浪费数小时；索引应作为共享只读资源。
9. 内存估算不足：多进程并行时索引内存叠加，导致 OOM；按「进程数 × 4GB」预留。
10. 混淆 CIGAR 的 S 与 H：软剪切保留序列、硬剪切截断，影响变异检测的坐标；理解区别再看结果。

小结

序列比对是生信流程的坐标基石，它的正确性直接决定下游分析的可信度。BWA 的 FM-index 与 MEM 种子策略是「速度与准确性权衡」的典范：用索引换速度，用启发式换可行性，用打分函数编码生物学先验。理解这些设计，你就能判断「什么时候该调参、什么时候该换工具」。

工程上，比对环节最该建立的习惯是「评估先行」：比对完先看 flagstat 的比对率、配对率、重复率，再看 MAPQ 分布，最后才进入下游。一个 90% 比对率的样本和一个 98% 的样本，后续分析的信任度完全不同。把比对质量作为样本的「健康档案」，是保证结论可靠的基础。

干净的数据加上正确的比对，下一步就是 SAM/BAM 与 samtools 实战 ——你将学会如何高效地查询、过滤、统计这些比对结果，为 变异检测 做好准备。如果你的比对率异常，回头检查 FASTQ 质控 那一环，往往能找到答案。