序列比对：从 FASTQ 到 SAM/BAM 的核心流程

序列比对，说白了就是把你手上那堆高通量测序的短读长（或者长读长）片段，怼到参考基因组上，看看每个 read 到底是从哪里来的。这一步做不好，后面的变异检测全是瞎猜。我个人习惯把比对看作整个 GATK 管线的基石——石头不稳，楼就塌。

BWA-MEM 算法：我为什么偏爱它

BWA-MEM 是目前最主流的短读比对工具之一，特别适合 Illumina 平台的数据。它用 BWT（Burrows–Wheeler 变换）+ FM 索引来做 seed-and-extend，能处理可变长度读长，也能应付一定程度的 gap 和软剪裁。我刚开始做 NGS 项目时试过其他工具，后来发现 BWA-MEM 在速度和敏感度之间平衡得最好——尤其对于人的全基因组数据，跑一次大概 4–8 小时（看机器配置）。

你要做的第一步，是先建索引。参考基因组 FASTA 文件需要先跑 bwa index：

bwa index -a bwtsw hg38.fa
# -a bwtsw 适合大基因组（>2G），小基因组可以用 is
# 生成 .amb .ann .bwt .pac .sa 等文件

然后比对命令很直接：

bwa mem -t 16 -R '@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA' \
    hg38.fa sample1_R1.fastq.gz sample1_R2.fastq.gz > sample1.sam

这里 -t 16 是开 16 个线程。-R 添加 Read Group 信息——这玩意儿千万不能省，GATK 后面要用它来区分样本和文库。我在项目里碰到过几次，有人没加 RG 标签，结果后面 HaplotypeCaller 报错说找不到样本，只能重跑，白白浪费两天。

关键参数记忆点：

• -k：最小 seed 长度，默认 19；对于长读长可以适当调大（比如 21–25）
• -M：将 split read 的辅助匹配标记为 secondary，Picard MarkDuplicates 处理时更干净
• -p：如果是单端或已合并的双端数据，用这个

SAM / BAM 格式：必须读懂的“身份证”

比对结果最原始的是 SAM 格式（文本，大得要命）。BAM 是它的二进制压缩版，也是实际管线里的标准格式。格式标准由 SAMspec 定义，但说实话，你不需要背全部字段，只要记住几个关键部分就好。

每一行代表一条 read 的比对记录。列字段顺序固定，一共 11 列是必选项：

FLAG 这一列是个 bitwise 的整数，每个比特代表一种特性。举个例子：

你想想看，一个 FLAG 值比如 99，怎么分析？

99 = 64 (read1) + 32 (mate是负链) + 2 (所有比对正确) + 1 (paired-end) = 0x63。查工具 samtools flags 可以帮你解码：

samtools flags 99
# 输出：0x63  99  PAIRED,PROPER_PAIR,MREVERSE,READ1

CIGAR 字符串也很重要，比如 150M 表示 150 个碱基完全匹配；50M10D100M 表示 50 匹配，10 个碱基删除，100 匹配。软剪裁（S）在 CIGAR 开头或结尾出现，表示那一段 read 比对不上但没被丢弃——常见于跨 junction 的 reads。

排序与压缩：samtools sort 与 samtools view

BWA-MEM 产出的 SAM 文件是按 read 顺序排列的，而 GATK 的后续工具（比如 MarkDuplicates、BaseRecalibrator）需要按坐标排序的 BAM。所以排序是格式化流水线里绕不过去的一步。

排序用 samtools sort，基本用法：

samtools sort -@ 8 -m 4G -o sample1.sorted.bam sample1.sam

• -@ 8：使用 8 个线程排序
• -m 4G：每个线程最多使用 4GB 内存（实际总内存约 32G，但还要给后面的处理留余量）
• -o：输出文件（自动检测格式，.bam 就是 BAM）

这里要留个心眼：排序会生成临时文件，默认放在当前目录，也可以用 -T /tmp/pfx 指定临时前缀，避免写满 /tmp。我清楚记得有次管线的 cron 任务炸了，就是因为 /tmp 被临时文件撑满——从那以后我习惯在排序命令里加 -T /data/tmp/xxx。

压缩级别方面，samtools 默认用 level 1（速度优先），你也可以用 --output-fmt BAM,level=5 来调。不过对于管线中间的中间文件，我通常保留 level 1；只有到最后存档时才用 level 9 重压缩一次。

排序完成后，一般还会用 samtools index 建立 .bai 索引，方便快速按坐标访问：

samtools index sample1.sorted.bam

顺便提一句：samtools view 可以随时从 BAM 里抽读子集，比如只看某个染色体的数据：

samtools view sample1.sorted.bam chr22:1-1000000

本章知识体系图

这张图把整个章节的核心流程串起来了。从左到右：原始 FASTQ 读长需要参考基因组索引（提前建好），然后跑 BWA-MEM 得到 SAM；SAM 再通过 samtools sort 按坐标排序，压缩成 BAM；同时建立 .bai 索引供快速随机访问。经过这一步后，我们就能得到一个干净的、有序的比对文件，直接喂给下一个环节——MarkDuplicates 去重。