📚 课程目录第4章变异注释工具选型

4. SAM/BAM处理：samtools与PICARD

变异解读的pipeline里，SAM/BAM文件是绕不开的中间产物。说白了，测序仪下机后的原始数据经过比对，得到的就是SAM或BAM格式。我个人习惯把SAM看成文本版，BAM看成压缩二进制版——两者内容一样，但BAM更省空间、读写更快。

一句话总结：后续的排序、去重、质量统计、变异检测，全都得先把数据转成BAM并处理好。

4.1 SAM/BAM格式速览

SAM文件长什么样？每行代表一条比对记录。我截一小段实际处理过的数据示例：

@SQ SN:chr1 LN:248956422
@SQ SN:chr2 LN:242193529
@PG ID:bwa PN:bwa VN:0.7.17
read001  99  chr1  10001  60  150M  =  10201  250  ATGC...  *  NM:i:0  MD:Z:150
read001 147  chr1  10201  60  150M  =  10001 -250  TACG...  *  NM:i:1  MD:Z:6A143

头部（@开头）包含参考序列名、程序信息等。比对行有11个必填字段：QNAME, FLAG, RNAME, POS, MAPQ, CIGAR, RNEXT, PNEXT, TLEN, SEQ, QUAL。嗯，这里要注意FLAG值，它用二进制位表示配对、链向、是否unmapped等信息，初看容易晕。

字段	含义	常见问题
FLAG	比对标志（整型）	用samtools flags解析，别自己算二进制
MAPQ	比对质量	低于20的read在变异检测中最好过滤掉
CIGAR	比对操作串	M=匹配/错配，I=插入，D=删除，S=软截断

为什么要有BAM？因为一个全基因组测序的SAM文件能上TB级，BAM压缩后一般只有200-300GB。我在项目里就吃过亏，初期图方便一直用SAM，结果磁盘爆了——从那以后，只要比对完成就立刻转BAM。

4.2 samtools核心命令

samtools是SAM/BAM处理的瑞士军刀。我个人用得最勤的就四个命令：view, sort, index, flagstat。下面一个个说。

view – 格式转换与过滤

# SAM → BAM
samtools view -bS input.sam > output.bam

# BAM → SAM（想看具体内容时用）
samtools view -h output.bam | less

# 只保留MAPQ≥30的比对
samtools view -q 30 -b input.bam > high_quality.bam

-b输出BAM，-S输入SAM（老版本需要，新版本自动判断）。-h表示输出头部，否则默认只输出比对行。我曾经漏掉头部导致下游工具报错，排查了半天……后来养成习惯，只要不是管道传数据，都会加-h。

sort – 按坐标排序

# 按染色体和位置排序（默认）
samtools sort input.bam -o sorted.bam

# 指定临时目录（大文件必备）
samtools sort -T /tmp/sort_temp -o sorted.bam input.bam

排序后的BAM才能建立索引，也是大多数变异检测软件的前提。-T指定临时文件前缀，我一般用/tmp/sort_$$来避免多任务冲突。排序本身挺占内存，记得设-m 2G控制每线程内存。

index – 建立索引

samtools index sorted.bam
# 会生成 sorted.bam.bai 文件

索引文件很小，但能实现快速随机访问。你想想看，几百GB的BAM，没有索引的话想截取chr1:100000-200000区域，得从头扫描——有了索引几乎瞬间定位。我在做CNV分析时，经常用samtools view指定区域来提取特定基因的覆盖情况。

flagstat – 比对统计

samtools flagstat input.bam

输出像这样：

1000000 + 0 in total (QC-passed reads + QC-failed reads)
990000 + 0 primary
5000 + 0 secondary
5000 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
970000 + 0 mapped (97.00% : N/A)
960000 + 0 primary mapped (96.97%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
...

这个统计是pipeline质量监控的第一道防线。我会重点看mapped比例和duplicates比例——如果mapped低于85%或者duplicates高于20%，一般就说明样本或建库有问题，得停下来检查。

小技巧：用samtools stats能得到更详细的指标，包括插入片段分布、碱基质量分布等，适合写入报告。

4.3 PICARD工具在变异解读中的应用

PICARD是Broad Institute维护的一套Java工具集，专门处理SAM/BAM/CRAM的元数据操作。在变异解读场景下，我最常用的是MarkDuplicates和FixMateInformation。

MarkDuplicates – 标记PCR重复

java -jar picard.jar MarkDuplicates \
  I=sorted.bam \
  O=dedup.bam \
  M=dedup_metrics.txt \
  REMOVE_DUPLICATES=false \
  ASSUME_SORTED=true

PCR重复会高估变异支持的深度，变异检测时必须标记出来。REMOVE_DUPLICATES=false表示只标记不删除，这样GATK能自动跳过重复读段。metrics文件里可以看到PERCENT_DUPLICATION，是文库复杂度的关键指标。我曾经遇到一个样本重复率高达45%，后来发现是建库时PCR循环数太多，直接要求重做。

注意：MarkDuplicates要求输入的BAM已经按坐标排好序。并且对单端read只能通过比对位置判断重复，对双端read同时检查位置和方向，更准确。

其他常用PICARD工具

工具	功能	使用场景
SortSam	排序	等价于samtools sort，但输入输出支持更多格式
BuildBamIndex	建立索引	等价于samtools index
CollectInsertSizeMetrics	插入片段统计	双端测序文库质量评估
ValidateSamFile	格式校验	检查BAM文件是否损坏或不符合规范
FixMateInformation	修复配对信息	某些比对软件可能没有正确设置配对标记

你可能想问：既然samtools能排序、索引，为什么还要用PICARD？其实两者可以混用，我自己的做法是：快速处理（比如只看比对统计）用samtools，需要记录详细元数据或做严格质量控制时用PICARD。

4.4 实战中常见的坑

讲几个我亲手踩过的坑：

排序后忘记建索引：好多工具（IGV、GATK）会直接报错“index not found”。建议养成写pipeline时自动在sort后紧跟index的习惯。
去重参数搞错：REMOVE_DUPLICATES=true会直接删掉重复read，这样后续统计的数据量就不对了，除非你非常清楚自己在做什么，否则用false只标记。
samtools view过滤过于严格：我一开始习惯猛加-q 30 -F 1024，结果丢失了很多低频变异信息。后来改为先做完整BAM，到GATK的VQSR阶段再自动过滤质量。
PICARD版本冲突：不同版本的PICARD对Java版本有要求。我曾在老服务器上装了新picard.jar，折腾半天才发现是Java版本不对。建议固定使用官方推荐的LTS版本组合。

避坑指南：我曾经在一个高深度靶向测序项目中，没做MarkDuplicates就直接跑GATK HaplotypeCaller，结果所有常见SNP位点的Allele Balance都严重偏离0.5，深度也比预期高一倍。加上MarkDuplicates后，所有指标回归正常。所以请记住——变异解读之前，BAM必须去重。

4.5 什么时候用samtools，什么时候用PICARD？

任务	推荐工具	理由
快速查看BAM头部	samtools view -H	不用等Java启动
格式转换 (SAM↔BAM)	samtools view	轻量、管道友好
排序	samtools sort	Python版本早就不推荐PICARD的SortSam了
标记重复	PICARD MarkDuplicates	行业标准，GATK官方推荐
插入片段统计	PICARD CollectInsertSizeMetrics	输出直方图数据，可直接用于质控报告
文件完整性校验	PICARD ValidateSamFile	检查BAM是否符合规范，比samtools -t更详细

其实在实际pipeline里，我一般把samtools放在前面做“粗加工”，PICARD放在中间做“精处理”。两者互不冲突，用好各自的优势就行。

个人建议：如果你是新手，先熟练samtools的四个常用命令；等pipeline跑顺了再加PICARD去重。别一上来就Java -jar，会把你劝退的。