Samtools：SAM/BAM格式操作全解

1.1 SAM 文件长什么样#

SAM（Sequence Alignment Map）是纯文本格式，分两部分：

1
Header（以 @ 开头，元数据）
2
@HD VN:1.6  SO:coordinate          # 格式版本 + 排序方式
3
@SQ SN:chr1 LN:248956422           # 染色体名 + 长度
4
@RG ID:sample1 SM:sample1          # Read Group 信息
5

6
Body（tab分隔的11个必填字段）
7
read1  0  chr1  100  42  50M  *  0  0  ATCG...  IIII...

1.2 CIGAR 字符串——看懂比对怎么对齐的#

CIGAR 是你排查比对问题最常用到的字段：

例子：35M1I14M = 35bp匹配 → 1bp插入 → 14bp匹配。

CIGAR 和 SEQ 长度必须一致。 计算规则：SEQ长度 = 所有 M+I+S 的碱基数之和。如果你手动编辑 SAM 文件忘了更新 CIGAR，samtools 会直接报错。

1.3 FLAG 位标志——不用死记，用工具解码#

FLAG 这一列让人头疼——一个数字编码了 12 种信息：

$FLAG = \sum_{i} bit_i \times 2^i$

常用的几个位：

不用手算。 samtools flags 帮你解码：

1
samtools flags 83
2
# 输出：
3
# 0x53  83  PAIRED,PROPER_PAIR,MREVERSE,READ1
4
# 解读：双端、正确配对、mate在反链、这是read1

常见组合记忆：

• 99 = read1 + 双端 + 正确配对 + mate反链 → RNA-seq Read1 的正常比对
• 147 = read2 + 双端 + 正确配对 + 自身反链 → RNA-seq Read2 的正常比对
• 4 = 未比对 → 这条read没比对上，但保留在文件里

2.1 view——转换格式和过滤#

1
# SAM→BAM（-b = BAM输出，-S = SAM输入）
2
samtools view -bS input.sam > output.bam
3

4
# BAM→SAM（人类可读）
5
samtools view -h input.bam | head -20
6

7
# 过滤：只保留比对成功且MAPQ≥30的reads
8
samtools view -b -F 4 -q 30 input.bam > filtered.bam
9

10
# 提取特定染色体区域
11
samtools view -b input.bam chr1:1000-5000 > region.bam
12

13
# 统计比对的reads数
14
samtools view -c -F 4 input.bam    # 总比对reads
15
samtools view -c -f 2 input.bam    # 正确配对的reads

参数详解：

• -F 4：排除 FLAG 含 bit 4（未比对）的 reads → 只要比对上的
• -f 2：只要 FLAG 含 bit 2（正确配对）的 reads
• -q 30：MAPQ≥30（错误率 < 0.1%）

MAPQ 阈值怎么选？

$P_{correct} = 1 - 10^{-MAPQ/10}$

生信流程通常用 -q 20（RNA-seq）或 -q 30（WGS 变异检测）。

2.2 sort——按坐标排序#

不排序的 BAM 文件很多下游工具拒绝处理。 排序方式有两种：

1
# 按染色体坐标排序（默认，最常用）
2
samtools sort input.bam -o sorted.bam
3

4
# 按read名称排序（用于某些特殊工具如Picard）
5
samtools sort -n input.bam -o name_sorted.bam
6

7
# 指定临时目录 + 多线程加速（大文件必用）
8
samtools sort -@ 8 -T /tmp/samtools_tmp input.bam -o sorted.bam

-T 指定临时目录很重要。samtools sort 需要额外的磁盘空间做归并排序——如果 /tmp 满了会直接报错。

2.3 index——创建索引#

索引让你可以快速访问 BAM 文件的任意区域：

1
samtools index sorted.bam
2

3
# 生成 sorted.bam.bai（或 sorted.bam.csi）
4
# .bai 支持到 2^29 bp 的基因组
5
# .csi 支持到 2^37 bp（超大基因组用）
6
samtools index -c sorted.bam   # 强制生成 .csi

有了索引才能用 samtools view chr1:1000-5000 这种区域查询，否则需要扫描整个文件。

2.4 flagstat——快速统计看数据质量#

这是我拿到比对文件后跑的第一个命令：

1
samtools flagstat sorted.bam

输出解读：

1
100000 + 0 in total (QC-passed reads + QC-failed reads)
2
 95000 + 0 primary mapped (95.00%)        # 比对率
3
 90000 + 0 properly paired (90.00%)       # 正确配对率

关键指标：

• 比对率（mapped rate） > 70% 正常，< 60% 可能参考基因组不对或污染
• 正确配对率（properly paired）接近比对率才正常——差异大说明文库插入片段异常
• duplicates 看 PCR 重复程度

2.5 stats——更详细的统计#

1
samtools stats sorted.bam > stats.txt
2

3
# 用 plot-bamstats 画图（samtools 自带的小工具）
4
plot-bamstats -p stats_plot stats.txt

输出里有几个关键指标：

1
SN  raw total sequences:    100000
2
SN  reads mapped:           95000
3
SN  reads properly paired:  90000
4
SN  insert size average:    250.5
5
SN  insert size standard deviation: 80.3

插入片段长度是质控指标之一——标准偏差太大说明文库质量不稳定。

2.6 idxstats——每个染色体多少 reads#

1
samtools idxstats sorted.bam

输出四列：染色体名、染色体长度、比对上的 reads 数、未比对的 reads 数。

做 RNA-seq 的注意： 如果你看到 Y 染色体有大量 reads 而样本是女性——说明可能有比对偏差或污染。

提取特定区域并统计深度#

1
# 看BRCA1基因区域的覆盖度
2
samtools depth -r chr17:43044295-43125483 sorted.bam | \
3
    awk '{sum+=$3; count++} END {print "Mean depth:", sum/count}'

筛出低质量比对用于排查#

1
# 找出MAPQ=0的reads（多重比对），看它们在哪些区域聚集
2
samtools view -q 0 -Q 1 sorted.bam | cut -f3 | sort | uniq -c | sort -rn | head

批量统计多个BAM的比对率#

1
for bam in *.bam; do
2
    mapped=$(samtools flagstat "$bam" | grep "mapped (" | head -1 | cut -d' ' -f1)
3
    echo "$bam: $mapped reads mapped"
4
done

坑1：samtools sort 报 “No space left on device”#

症状：sort 到一半崩溃。

原因：samtools sort 需要临时磁盘空间，大约是 BAM 文件的 1.5-2 倍。/tmp 在根分区且空间不够。

1
# 检查 /tmp 剩余空间
2
df -h /tmp
3

4
# 指定大磁盘作为临时目录
5
samtools sort -T /data/tmp_sort -@ 8 input.bam -o sorted.bam

坑2：BAM 文件很大但 flagstat 显示 0 mapped#

原因：SAM 转 BAM 时忘了加 -h 参数，header 丢了。没有 SQ header，samtools 不知道有哪些染色体。

1
# 正确做法：加 -h 保留 header
2
samtools view -bSh input.sam > output.bam
3

4
# 检查 header
5
samtools view -H output.bam | head

坑3：idxstats 显示某染色体 reads 数是 0 但 IGV 里能看到#

原因：BAM 文件里的染色体名跟 FASTA 索引的不一致。比如 BAM 里写 chr1，但索引里是 1（没前缀）。

1
# 看看你的 BAM 用了什么染色体名
2
samtools idxstats sorted.bam | head -5
3

4
# 看看参考基因组的染色体名
5
grep "^>" ref.fa | head -5

如果确实不一致，用 sed 统一或重新比对。

坑4：samtools view 过滤后 reads 数不符合预期#

比如 samtools view -c -F 4 -f 2 想统计”正确配对的reads”，结果比 -F 4 少很多。

原因：-F 4 和 -f 2 是不同级别的过滤——很多 reads 虽然比对上了（-F 4），但 mate 的比对质量差导致不满足”proper pair”条件。这是正常的，关键是看比例是否合理。