SAM/BAM 格式与 samtools 实战

引言

BAM 是生信世界的中心格式：上游比对产出 BAM，下游变异检测、覆盖度分析、可视化全都消费 BAM。可以说，掌握 BAM 与 samtools，就掌握了生信数据处理的基本功。但很多人只会背几个命令（samtools sort、samtools index），对格式本身一知半解，遇到问题时束手无策。

SAM/BAM 格式的设计堪称工程典范：SAM 是纯文本、人类可读，便于调试；BAM 是二进制压缩版，体积小、支持随机访问。两者语义完全等价，可以无损互转。格式里承载的信息远超「序列 + 坐标」——FLAG 标志位编码了 12 种状态，CIGAR 描述了比对细节，可选标签（TAG）携带了 MAPQ 之外的质量信息。

工程上的难点有三个：一是过滤条件复杂，samtools view -F 3844 这种魔法数字没人记得住；二是大文件性能，一个 80 GB 的 BAM 随便扫一遍就是几十分钟，必须用索引做随机访问；三是并行化，samtools 的单线程瓶颈常被忽略。本文会把这些讲透。

本文按「格式解析 → FLAG → 索引 → 命令全景 → 过滤 → 统计 → 覆盖度 → 性能」的顺序展开。目标是让你读完能「看着 BAM 就知道它是否健康，遇到问题知道用哪个命令定位」。命令基于 samtools 1.19/1.20。

目录

1. SAM 格式逐字段解析
2. FLAG 标志位详解
3. BAM 编码与索引机制
4. samtools 命令全景
5. 过滤语法与子集提取
6. 统计与质量指标
7. 覆盖度分析
8. 格式转换与互操作
9. 大文件性能与并行

1. SAM 格式逐字段解析

SAM 文件由「头部（header）」和「比对记录（alignment）」两部分组成。头部以 @ 开头：

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:248956422
@RG	ID:sample1	SM:sample1	PL:ILLUMINA	LB:lib1
@PG	ID:bwa	PN:bwa	VN:0.7.17	CL:bwa mem ...

比对记录有 11 个必选字段加可选标签：

r001  99  chr1  7  60  8M2I4M1D3M  =  37  39  TTAGATAAAGGATACTG  *

可选标签中的 NM（编辑距离）、MD（错配详情）、AS（比对得分）、XS（次优得分）对下游工具至关重要。GATK 依赖 NM 和 MD 来判断读段与参考的差异。如果比对工具没写这些标签，下游会报错。

理解字段含义的实用价值：当 samtools view 输出一行看不懂时，逐字段对照表格就能解析。例如 TLEN=39 表示这对读段的插入片段长 39 bp（很短，提示片段异常）；MAPQ=0 表示多重比对。

2. FLAG 标志位详解

FLAG 是一个位掩码（bitmask），每一位代表一种状态。这是 SAM 格式最「劝退」的部分，但也是最有用的信息：

常见 FLAG 值速查：

99  = 1+2+32+64  → 配对+正确配对+配对负链+R1（正向主比对）
147 = 1+2+16+128 → 配对+正确配对+负链+R2（反向主比对）
4   = 未比对
16  = 单端负链
3844 = 4+2048+1024+512+256 → 未比对+补充+重复+QC失败+次要

过滤的关键：搞清「要什么」而非「排除什么」。samtools 的 -f（保留含任意位）和 -F（排除含任意位）是两个方向：

# 只要「正确配对的主比对读段」→ 排除未比对/次要/补充/重复/QC失败
samtools view -b -F 3844 aln.bam > clean.bam
# 3844 = 4(未比对) + 256(次要) + 512(QC失败) + 1024(重复) + 2048(补充)

# 只要 R1 正向主比对
samtools view -b -f 67 aln.bam > r1.bam   # 67 = 1+2+64

# 只要未比对读段（用于污染检测）
samtools view -b -f 4 aln.bam > unmapped.bam

# 排除重复但保留其他
samtools view -b -F 1024 aln.bam > nodup.bam

-F 3844 是「干净 BAM」的标准过滤，值得记住。注意这里用的是「排除」语义——因为「未比对」「次要」「补充」这些位在正常读段里都是 0，用 -F 排除最方便。

补充说明 256（secondary）与 2048（supplementary）的区别：secondary 是「读段的另一个可能位置」（多重比对），supplementary 是「读段被拆成多段分别比对」（嵌合或剪接）。变异检测通常两者都排除，但结构变异检测需要保留 supplementary（因为 SV 会让读段拆分）。

3. BAM 编码与索引机制

BAM 是 SAM 的二进制压缩版，采用 BGZF（Blocked GZIP Format）压缩。BGZF 的关键特性是「分块」：每 64 KB 一个块，每块独立压缩。这带来两个好处：

1. 支持随机访问：可以只解压需要的块，而非整个文件；
2. 支持并行：多个块可以并行解压，这是 samtools 多线程的基础。

但随机访问需要「知道某个坐标的数据在哪个块」，这就是**索引（.bai 或 .csi）**的作用。索引本质是一个「参考区间 → 虚拟文件偏移」的映射表。

samtools index aln.bam          # 生成 aln.bam.bai（默认）
samtools index -c aln.bam       # 生成 .csi（支持超长染色体，> 512 Mb）

# 有索引后可以按区间快速提取（无需扫描全文件）
samtools view aln.bam chr1:1000000-2000000 > region.bam
samtools view -c aln.bam chr7:117480000-117550000   # 数 EGFR 区域读段数

索引的两个前提：BAM 必须坐标排序，且索引必须与 BAM 同步更新（BAM 改了索引就失效）。忘记建索引会导致下游工具报错（“index not found”）或退化为全文件扫描。

.bai 与 .csi 的选择：.bai 是传统格式，坐标上限约 512 Mb；.csi 支持更大坐标。人类染色体都小于 512 Mb，用 .bai 即可；但某些物种（如某些植物、蝾螈）染色体超大，必须用 .csi。混用会导致「染色体太长无法索引」的错误。

4. samtools 命令全景

samtools 有 20 多个子命令，按用途分类记忆：

最常用的几个组合：

samtools quickcheck *.bam              # 检查所有 BAM 是否完整（CI 必用）
samtools flagstat aln.bam              # 快速看比对率
samtools idxstats aln.bam              # 按染色体看读段分布（查性别、污染）
samtools sort -@ 8 -m 2G -o s.bam in.bam
samtools view -b -q 20 -F 3844 aln.bam > filtered.bam
samtools depth -a -q 20 aln.bam | awk '$3>=10' | wc -l   # 统计 ≥10x 的位点数

samtools quickcheck 值得强调：它快速检查 BAM 的 EOF 标记（BGZF 块完整性）。如果 BAM 在传输或写入过程中被截断，quickcheck 会报错。把这一步放进流程的入口检查，能拦住「文件损坏」这类隐蔽问题。

5. 过滤语法与子集提取

samtools view 的过滤参数是日常使用频率最高的：

实用的子集提取示例：

# 提取特定样本的读段（多样本混合 BAM）
samtools view -b -r sample1 merged.bam > sample1.bam

# 提取目标区域（如 panel 基因）
samtools view -b -L panel.bed aln.bam | samtools sort -o panel.bam -

# 取 1% 读段做快速测试（固定种子保证可复现）
samtools view -b -s 0.01 --subsample-seed 42 aln.bam > sample.bam

# 排除特定染色体（如 decoy contig）
samtools view -b aln.bam $(samtools idxstats aln.bam | cut -f1 | grep -v decoy | tr '\n' ' ')

一个常见需求是「提取特定基因的读段用于可视化或验证」。组合 -L（区域）与 -q（质量）过滤即可。注意 -L 保留的是「与区间重叠」的读段，包括部分重叠的，这对 IGV 可视化是对的，但做精确计数时可能高估。

6. 统计与质量指标

samtools 的统计命令输出大量指标，关键是知道看哪些：

# flagstat：高层汇总
samtools flagstat aln.bam
# 42263840 + 0 in total (QC-passed reads + QC-failed reads)
# 40214730 + 0 mapped (95.15% : N/A)          ← 比对率
# 38802650 + 0 properly paired (91.81% : N/A)  ← 正确配对率
# 5120340 + 0 duplicates (12.11% : N/A)        ← 重复率

# stats：详细统计（适合程序解析）
samtools stats aln.bam > aln.stats
grep ^SN aln.stats | cut -f2- | head -30

# idxstats：按染色体分布（查性别、污染）
samtools idxstats aln.bam
# chr1  248956422  4000000  100000
# chrX  156040895  200000   80000     ← X 覆盖高、Y 覆盖低 → 可能是女性
# chrY  57227415   5000     200

idxstats 是排查「样本调包」「性别不符」的利器：人类女性样本的 X 染色体覆盖度约为常染色体的一半（因为女性有两条 X），Y 染色体几乎无覆盖；男性则 X、Y 都是单拷贝。如果元数据记录是男性但 idxstats 显示 Y 无覆盖，说明样本可能调包或污染。

samtools stats 输出的 error rate、average length、insert size average 等指标也很关键。特别是插入片段（insert size）：正常 PE 文库应在 300-600 bp，若显示 100 bp 或双峰，提示建库异常。

7. 覆盖度分析

覆盖度（coverage/depth）是判断数据能否支撑分析的核心指标。两种工具：

# samtools depth：逐碱基深度（简单但慢，大 BAM 慎用）
samtools depth -a -q 20 -Q 20 aln.bam > depth.txt
#  -a: 输出零覆盖位点   -q: 最小 MAPQ   -Q: 最小碱基质量

# mosdepth：更快，支持按区间汇总
mosdepth --by 1000 --fast-mode -t 4 sample aln.bam
# 输出 sample.mosdepth.summary.txt（总览）
# 输出 sample.per-base.bed.gz（逐碱基）
# 输出 sample.regions.bed.gz（按 1kb 窗口）

mosdepth 比 samtools depth 快数倍（用更高效的数据结构），是大规模覆盖度分析的首选。它输出的 summary 包含关键指标：

chrom   length    bases    mean   min   max
total   3099922541  95000000000  30.65  0  1234

用覆盖度做质控的实用脚本：

# 计算「覆盖 ≥ 10x 的碱基比例」
zcat sample.per-base.bed.gz | awk '{n++; if($4>=10) c++} END {print c/n}'

# 按外显子区间统计（WES 关键指标）
mosdepth --by targets.bed wes_sample aln.bam
zcat wes_sample.regions.bed.gz | awk '{if($4>=30) c++; n++} END {print c/n}'

覆盖均匀性是比平均深度更重要的指标。平均 30x 但分布极不均匀（部分区域 100x、部分 0x）的数据，实际可用性远不如均匀的 25x。检测方法是看深度的标准差，或直接看「覆盖 ≥ 10x 的位点比例」（人类 WGS 应 ≥ 90%）。

8. 格式转换与互操作

BAM 与各种格式的转换是日常操作：

# BAM ↔ SAM（文本）
samtools view -h aln.bam > aln.sam        # BAM → SAM
samtools view -b aln.sam > aln.bam        # SAM → BAM

# BAM → FASTQ（重新提取原始读段）
samtools fastq -1 R1.fq -2 R2.fq aln.bam  # PE 分开
samtools fastq aln.bam > reads.fq         # SE

# BAM ↔ CRAM（归档压缩）
samtools view -C -T ref.fa aln.bam > aln.cram   # BAM → CRAM
samtools view -b -T ref.fa aln.cram > aln.bam   # CRAM → BAM

# BAM → BED（区间）
bedtools bamtobed -i aln.bam > aln.bed

# BAM → 覆盖度轨道（可视化）
bamCoverage -b aln.bam -o coverage.bw --binSize 25   # deepTools

CRAM 的取舍：比 BAM 小 40-50%，但解码需要参考序列（-T ref.fa）。这带来一个工程问题——CRAM 无法独立存在，必须配套参考基因组。归档时用 CRAM 省空间，但分析时通常转回 BAM（因为很多工具对 CRAM 支持不完整）。参见 Lustre 并行文件系统 了解大文件存储策略。

samtools fastq 的一个重要用途是「读段回收」：当发现比对参数不对想重跑时，可以从 BAM 提取原始 FASTQ（前提是 BAM 未做质量重校准、且保留了原始序列）。但如果 BAM 经过了 BQSR，质量值已被修改，提取的 FASTQ 不等于原始数据。这是「中间产物依赖」的典型陷阱。

9. 大文件性能与并行

一个 80 GB 的 BAM，随便一条命令就是几十分钟。性能优化要点：

一是善用索引做随机访问。全文件扫描（samtools view aln.bam）是 O(文件大小)，而区间查询（samtools view aln.bam chr1:...）用索引只读相关块，快几个数量级。凡是「只关心某区域」的操作，一定加区间。

二是合理设置线程。samtools 的多数命令支持 -@：

samtools sort -@ 8 -m 2G -o s.bam in.bam
#  -@ 8: 8 个线程
#  -m 2G: 每线程 2GB 内存（总内存 = 线程数 × m）
#  排序内存需求 = 数据量 × 常数，内存不足会写临时文件变慢

排序是最吃资源的操作。-m 设大能减少临时文件，但内存叠加要控制——-@ 8 -m 2G 意味着最多用 16 GB。用 -T /scratch/tmp 把临时文件放本地 NVMe，避免打到共享存储。

三是避免管道中的单线程瓶颈。bwa mem | samtools sort 这样的管道，如果 sort 比 mem 慢，会拖慢整体。给 sort 分配足够线程（-@ 8）能缓解。更激进的优化是用 sambamba（多线程更充分）替代 samtools sort。

四是分染色体并行。对于可分割的任务（如按染色体做变异检测），可以按染色体拆分 BAM 并行处理：

for chr in $(samtools idxstats aln.bam | cut -f1 | grep -v '^\*'); do
  samtools view -b aln.bam $chr | samtools sort -o ${chr}.bam - &
done
wait

这种「按染色体 scatter-gather」是生信并行的经典模式，流程引擎（Nextflow/Snakemake）对它有原生支持。

权衡取舍

常见坑清单

1. 忘记建索引：下游工具报 “index not found” 或退化全扫描；BAM 生成后立即 samtools index。
2. BAM 未排序就索引：索引失败或结果错误；索引前必须坐标排序（@HD SO:coordinate）。
3. BAM 改动后索引失效：过滤/合并后忘了重建索引；每次生成新 BAM 都重建。
4. -F 3844 记不住：混淆过滤方向导致漏掉或保留错误读段；记住 3844 = 4+256+512+1024+2048。
5. 混淆 secondary 与 supplementary：SV 检测需保留 supplementary，误删会丢结构变异信号。
6. samtools depth 用于大 BAM：慢到不可接受；改用 mosdepth。
7. 忽略 idxstats 的性别校验：样本调包到分析后期才发现；早期用 X/Y 覆盖比例核对。
8. BQSR 后提取 FASTQ：质量值已被修改，不等于原始数据；需要原始数据应从 FASTQ 备份取。
9. 排序内存溢出：-m 与 -@ 组合超出节点内存；按「线程数 × m」预估总占用。
10. CRAM 脱离参考：归档 CRAM 时忘了记录参考版本，日后无法解码；把参考路径写入元数据。

小结

SAM/BAM 是生信数据处理的通用语言，samtools 是这门语言的「标准库」。理解格式字段（尤其是 FLAG 与 CIGAR）、索引机制与过滤语法，是从「会敲命令」到「能解决问题」的分水岭。一个能熟练用 samtools 定位问题的人，在生信工程里永远不愁没活干。

工程上最该建立的三个习惯：一是「索引随 BAM 走」，任何新生成的 BAM 立即建索引；二是「先统计后分析」，用 flagstat/idxstats/stats 摸清数据底细再往下走；三是「能随机访问就不全扫描」，用区间查询把大文件的处理成本降到最低。

掌握了 BAM 的操作，你就可以进入 变异检测与 VCF 处理 ——从比对记录到变异位点的关键一跃。而在整个流程层面，如何把这些 samtools 步骤编排成可复现的管线，则是 生信流程编排 的主题。