4. 比对：从BWA-MEM到分析就绪的BAM
聊到GATK流程，比对这一步是绕不开的。很多人觉得拿个工具跑一下就行，但这里面的门道，说白了就是“快”和“准”的平衡。我早期做全外显子组项目时，就因为参数没调好，多花了三天时间重跑。嗯，这一节我们把BWA-MEM的底细、SAM/BAM/CRAM的差异、还有samtools那几条关键命令彻底搞明白。
4.1 BWA-MEM：比对引擎的“多面手”
BWA-MEM，全称Burrows-Wheeler Aligner's Maximal Exact Match。听名字挺唬人，其实就是把短读长（100bp-1Mbp）映射到参考基因组上的算法。我个人习惯把它当成“默认引擎”，因为它在处理长读长、嵌合体read和结构变异时比老版BWA-backtrack稳得多。
核心思路分三步：
- 种子查找：先把read切成小段，在参考基因组上用BWT索引找精确匹配的“种子”。
- 链式延伸：把这些种子在基因组上连成更长的匹配链，允许间隙和错配。
- 配对与评分：对双端read分开比对，再用SW算法做局部优化，最后按MAPQ（比对质量）打分。
这里有个经验：BWA-MEM默认对短读长效果不错，但如果你跑的是150bp以上的长读长，记得加-M（标记短拆分）和-T（降低阈值）。我曾经在一条250bp的模拟数据上没加-M，结果下游GATK报错说拆分标记不对……排查了一下午。
常用比对命令
bwa mem -t 16 -M -R '@RG\tID:sample1\tLB:lib1\tPL:ILLUMINA\tSM:sample1' \
hg38.fa sample1_R1.fastq.gz sample1_R2.fastq.gz > sample1.sam-t 16表示用16线程。别超CPU物理核数，别问我怎么知道的……
我补充一个避坑指南：Read Group（@RG）一定要写全。ID、LB、PL、SM四个字段缺一不可。早期GATK版本会直接报错，新版虽然能过但会给你留一堆坑。
注意：BWA-MEM输出的是SAM格式，如果中间文件太大，可以用samtools view -bS直接转BAM，省下一半磁盘空间。
4.2 SAM/BAM/CRAM：三种格式的取舍
比对完拿到的SAM（Sequence Alignment/Map）是纯文本，一行一条记录，包含11个必填字段和若干可选tag。你想想看，一个全基因组30×的测序数据，SAM文件轻松上几百GB。所以没人会拿SAM做分析，大家用BAM或CRAM。
| 格式 | 类型 | 压缩比 | 随机访问 | 兼容性 |
|---|---|---|---|---|
| SAM | 文本 | 1× | 无 | 最好 |
| BAM | 二进制 | ~5×-8× | 需索引 | 好 |
| CRAM | 二进制 | ~10×-20× | 需索引+参考基因组 | 稍差 |
我个人的建议：长期存储用CRAM，日常分析用BAM。CRAM依赖参考基因组文件来压缩碱基质量，能省下不少空间。但如果你经常在不同机器之间传文件，BAM的兼容性更省心。
说到SAM的必填字段，我整理一下常见的那11列：
- QNAME：read名字，一对read也一样
- FLAG：用数字表示比对状态，比如1是配对、2是正确匹配、4是自己没比上
- RNAME：参考序列名，比如chr1
- POS：比对起点的1-index坐标
- MAPQ：比对质量，越高越可靠
- CIGAR：描述比对形状，比如150M表示完美匹配150bp
剩下的我就不列全了，但FLAG那个数字你最好记住怎么看。我以前写过一个小脚本，每次跑完都自动解析FLAG，防止配对比对信息不对。
4.3 排序与索引：没有排序的BAM就是一堆废数据
从BWA-MEM出来的BAM是按read名字排的，而GATK要求按基因组坐标排序。为什么？因为下游的IndelRealigner、HaplotypeCaller都需要在连续基因组区间内快速查找read。你想想看，如果不排序，每次找一个小区域都得扫一遍全文件——那速度，能慢到让你抓狂。
排序的命令很简单：
samtools sort -@ 8 -m 4G -o sample1.sorted.bam sample1.bam-@ 8表示8线程并行排序，-m 4G表示每个线程允许4G内存。我建议根据机器物理内存调这个值，比如128G内存，可以设-m 8G -@ 12，排序会快很多。
经验之谈：如果磁盘IO是瓶颈，可以指定一个临时目录用SSD：-T /scratch/tmp。我曾经在机械硬盘上排序500GB的BAM，耗时8小时；换SSD后不到2小时。
排序完后必须建索引，不然没法随机访问：
samtools index sample1.sorted.bam索引文件.sorted.bam.bai就是坐标和文件偏移量的映射表。GATK每请求一个区间，samtools通过这个索引直接跳到对应位置，不用遍历全文件。
我记得有次我忘了做samtools index，直接丢进GATK流程，结果报“BAM index file not found, cannot parse”。排查半天才发现——嗯，就少了一行命令。
4.4 比对统计：samtools flagstat
好不好，看数据。比对完后第一件事就是盯着flagstat的输出看几秒。它能告诉你比对的整体质量，帮你快速判断是不是样本污染、接头没去干净、或者参考基因组选错了。
命令行：
samtools flagstat sample1.sorted.bam输出示例（我手改过数值）：
13367344 + 0 in total (QC-passed reads + QC-failed reads)
13041792 + 0 primary
0 + 0 secondary
325552 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
13143216 + 0 mapped (98.32% : N/A)
12817664 + 0 primary mapped (98.28%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)我重点关注这几项：
- mapped (%)：正常在95%-99%。如果低于90%，可能样本有问题或者参考基因组错了。
- properly paired (%)：双端正确配对的比例。全基因组数据应该在90%以上。我跑过一个全外显子，properly paired只有75%，后来发现是因为文库构建时片段太长，BWA-MEM默认参数不匹配。
- secondary + supplementary：这两个标记多说明read在基因组上有多个比对位置，可能是重复区域或结构变异。太多的话，下游GATK会剔除一部分，你得心里有数。
我用个实战场景：某次跑单细胞RNA-seq数据，flagstat显示mapped只有80%。一开始以为细胞污染，后来发现是参考基因组里转了重复序列masked的版本，导致大量read比不上。换回unmasked版本，mapped立刻升到95%。这种事，嗯，我经历过不止一次。
下面这张图展示从原始FASTQ到分析就绪BAM的完整流程：
理解了这张图的每一步，比对环节就算吃透了。排序和索引虽然是“体力活”，但每次漏掉都会让你多花半小时排查。我建议你在脚本里把这几个步骤写成一个管道，一键产出sorted.bam + .bai + flagstat报告，省心又标准化。