4. 比对到参考基因组
群体联合基因分型的第一步——把测序 reads 放到参考基因组上。这一步要是做不好，后面所有变异检测都是空中楼阁。我个人习惯用 BWA-MEM，它对于 Illumina 双端测序数据来说，速度和准确度都算得上首选。你想想看，一个完整 WGS 样本运行时间在几小时内，而且对复杂区域的 handling 也比旧版 BWA 好。
我的经验：新手往往只看比对率，其实比对的均匀性、插入片段分布同样关键。我曾经遇到过一个样本比对率 98%，但覆盖深度严重不均，最后基因分型结果一塌糊涂。
4.1 BWA-MEM 核心参数与性能调优
BWA-MEM 的命令行结构很清晰，但参数调整直接影响结果。
| 参数 | 说明 | 推荐值 / 说明 |
|---|---|---|
-t |
线程数 | 4-16（根据服务器核心数，我一般设物理核数的一半） |
-M |
将 split reads 标记为 supplementary | 必须加，否则后续 MarkDuplicates 会出问题 |
-R |
Read group 信息（@RG 标签） | 必须，格式：@RG\\tID:样本名\\tSM:样本名\\tLB:文库名\\tPL:ILLUMINA |
-T |
最小比对分数阈值 | 默认 30，对长 reads 可以降到 20 |
-k |
种子长度 | 默认 19，reads 短时可减小 |
-w |
gap 开放罚分位移 | 默认 5，一般不动 |
下面是我项目里常用的完整命令：
bwa mem -t 8 -M -R '@RG\tID:NA12878\tSM:NA12878\tLB:lib1\tPL:ILLUMINA' \
/ref/hg38.fa \
sample_R1.fastq.gz sample_R2.fastq.gz \
> sample_aln.sam
避坑指南：我曾经漏加
-M 参数，结果在去重复时发现大量 reads 被误标为 duplicate。另外 -R 里的 \\t 一定要转义成制表符，否则 read group 解析会失败。你可以在命令行里用 echo 测试一下格式。
4.2 SAM/BAM 格式关键字段
SAM（Sequence Alignment/Map）是文本格式，BAM 是其压缩二进制版本。你得看懂几个核心列，否则后续过滤无从下手。
| 列号 | 字段 | 含义 |
|---|---|---|
| 1 | QNAME | read 名称 |
| 2 | FLAG | 位标志，描述比对状态（paired, mapped, proper pair 等） |
| 3 | RNAME | 参考序列名称（染色体/contig） |
| 4 | POS | 1-based 比对起始位置 |
| 5 | MAPQ | 比对质量（Phred 分数，越高越好） |
| 6 | CIGAR | 比对操作字符串（如 150M 表示完全匹配 150bp） |
FLAG 的解析我经常用 samtools flags 命令快速查看。举个例子：FLAG=99 表示“双端测序、第一端已比对、插入片段方向正确”。
4.3 从 SAM 到 BAM：排序与索引
原始 SAM 文件太大，而且未排序。我们需要把它转成 BAM 并排序、建索引，这几乎是标准流水线：
# SAM -> BAM
samtools view -bS sample_aln.sam > sample_aln.bam
# 排序（按坐标）
samtools sort -@ 4 -m 4G -o sample_sorted.bam sample_aln.bam
# 建索引
samtools index sample_sorted.bam
# 一步到位（省空间）
bwa mem -t 8 -M -R '@RG\tID:NA12878\tSM:NA12878\tLB:lib1\tPL:ILLUMINA' \
/ref/hg38.fa R1.fq R2.fq \
| samtools sort -@ 4 -m 4G -o sample_sorted.bam
samtools index sample_sorted.bam
我的习惯：我会直接用管道把
bwa mem 的输出接到 samtools sort，省掉中间磁盘写 SAM 的步骤。注意 -m 4G 表示每个线程 4GB 内存，乘以线程数就是总内存，别把服务器撑爆了。
4.4 比对质量的评估与过滤
拿到排序后的 BAM，先跑一遍 flagstat 和 idxstats，快速判断整体健康度：
samtools flagstat sample_sorted.bam
samtools idxstats sample_sorted.bamflagstat 输出示例（关键行）：
1000000 + 0 in total (QC-passed reads + QC-failed reads)
980000 + 0 mapped (98.00% : N/A)
910000 + 0 properly paired (91.00% : N/A)
50000 + 0 with mate mapped to a different chr
30000 + 0 with mate mapped to a different chr (mapQ>=5)我通常关注的几个阈值：
- 总体比对率 ≥ 95% （全基因组测序）
- properly paired ≥ 90%
- 不同染色体配对 ≤ 5% （过高提示样本污染或参考基因组问题）
- duplicate rate ≤ 15% （用
MarkDuplicates统计）
实战技巧：如果比对上的是线粒体、核糖体等重复区域，比对率会虚高。我习惯再用
samtools view -q 30 过滤低质量比对后重新统计，才更能反映真实情况。
下面这张图概括了从原始 FASTQ 到排序索引 BAM 的整个比对流程：
整个流程走下来，你手里就有了一份干净的排好序的 BAM 文件和对应的 .bai 索引文件。这个 BAM 就是后续 GVCF 联合分析的原料了。
一个常见的坑：很多人忘记建索引就直接跑 GATK，结果报错“Index not found for bam file”。另外如果 BAM 文件很大（>100GB），建议用
samtools index -c 生成 CSI 索引，否则普通 .bai 索引只能支持参考序列长度小于 2^29。我在处理人类全基因组数据时遇到过这个问题。
好了，这一章我们把 BWA-MEM 比对、SAM/BAM 转换、排序索引和基本质控都过了一遍。这些命令你只要写进脚本，配合好参数，基本不会出大问题。下一阶段就是把这些 BAM 输入到 GATK HaplotypeCaller 中生成 GVCF 文件了。但我们先把这一步做扎实——地基不稳，楼盖得再高也白搭。