📚 课程目录第4章 SAMtools处理

第4章 SAMtools处理：排序、去重、索引

拿到原始的SAM文件，说实话，那基本是一锅粥。GATK的HaplotypeCaller虽然能直接吃SAM，但效率低得吓人。

我个人习惯是，先把SAM转成BAM，然后快速做个排序和去重。这一步做不好，后面的变异检测质量会大打折扣。

4.1 SAM与BAM：为什么必须转二进制

SAM是文本格式，BAM是它的二进制压缩版。区别是什么？举个简单的例子：一个30X覆盖度的全外显子组SAM文件，动不动就几十个GB。BAM呢？大概只有它的五分之一大。

我刚开始做生物信息时，就用SAM直接跑流程。结果服务器硬盘爆了，老板差点让我走人。从那以后，我再也没碰过纯文本的比对文件。

转换命令很简单：

samtools view -bS sample.sam > sample.bam

这条命令里，-b指定输出BAM格式，-S告诉软件输入是SAM。我建议你把这个命令记牢，它几乎天天用。

  核心思路： 文本转二进制，节省空间 + 提升随机读取速度。
  这个转换几乎是免费的，没有信息损失，只有好处。

4.2 排序：为索引和下游分析铺路

未经排序的BAM文件，记录是按比对顺序排列的。这种顺序对于需要按基因组位置查询的工具来说，简直是噩梦。

你想啊，HaplotypeCaller要处理某个区域，它得把整个文件扫一遍。这效率，谁能受得了？

排序之后，所有记录按照染色体和位置排列。这样，GATK就能快速定位到指定区域，然后干活。

排序命令：

samtools sort -o sample.sorted.bam sample.bam

这里有个坑。我以前在集群上跑任务，内存给得少，结果排序跑到一半就挂掉了。为什么？因为数据太大，内存不够用。

后来我学会了用-m参数控制每个线程的最大内存，再用-@控制线程数：

samtools sort -m 2G -@ 4 -o sample.sorted.bam sample.bam

这个配置的意思是：每个线程用2GB内存，用4个线程并行排序。这样写，适合大多数场景。

注意： 排序后的文件名一定要记得加.sorted标记。别问为什么，问就是我在几十个BAM文件里找原始文件找到想哭。

4.3 去重：扔掉PCR的痕迹

PCR扩增是文库构建的必备步骤，但它有个麻烦：它会产生重复的read。这些重复read在变异检测时会造成假阳性。

嗯，这里要讲清楚。GATK的标记重复（MarkDuplicates）和SAMtools的去重（rmdup）道理一样，但用途不太一样。

我个人习惯用GATK的MarkDuplicates，因为它更精细，能处理更多边缘情况。不过，如果你只想快速去除PCR重复，SAMtools也够用：

samtools rmdup sample.sorted.bam sample.dedup.bam

但是！注意了，这个命令只能处理单端数据。如果是双端数据，得加-s参数：

samtools rmdup -s sample.sorted.bam sample.dedup.bam

为什么双端要加-s？因为双端read的左右两条链位置不同，去重的逻辑更复杂。不加这个参数，软件会错误地把很多真正的read当成重复丢掉。

我在一个项目中，就吃过这个亏。当时用了默认参数去重，结果造成大量假阴性——该检出的变异没了。后来重新跑了一遍，才发现是去重参数不对。

避坑指南： 如果你准备用HaplotypeCaller做变异检测，我推荐你用GATK的MarkDuplicates替代SAMtools rmdup。它会在BAM文件中标记重复，而不是直接删除。这样，HaplotypeCaller拿到BAM时，可以自己决定怎么处理它们。这更安全。

4.4 索引：快速定位的钥匙

排序之后，最后一步就是建立索引。索引文件是个小文件，它记录了BAM中每个bin的偏移量。

有什么用？有了索引，GATK就可以直接跳到指定区域读取数据，而不需要从头扫描整个文件。

说白了，就是给基因组装了个“书签”。

建立索引的命令：

samtools index sample.dedup.bam

这条命令会生成一个名为sample.dedup.bam.bai的索引文件。注意，索引文件必须和BAM文件在同一个目录，而且名字必须完全对应。

为什么会这样？我猜GATK是按文件名自动查找索引的。如果找不到，它就会报错：

[E::hts_open_format] fail to open file sample.dedup.bam.bai

这时候别慌。检查一下文件名是不是匹配，是不是有拼写错误。我见过有人把.bai写成了.bam.bai，结果怎么也查不到。

  一句话总结：
  排序、去重、索引。这三个步骤，每一步都决定了你后续的HaplotypeCaller能不能跑得快、跑得准。
  别偷懒。

4.5 完整流程与SVG图

上面每一步都讲清楚了，现在我把它们串成一条完整的命令链。你可以直接复制到脚本里用：

samtools view -bS sample.sam > sample.bam
samtools sort -m 2G -@ 4 -o sample.sorted.bam sample.bam
samtools rmdup -s sample.sorted.bam sample.dedup.bam
samtools index sample.dedup.bam

这样跑完之后，你会得到三个文件：sample.dedup.bam和sample.dedup.bam.bai。这两个文件，就是GATK HaplotypeCaller的标准输入。

我把这个流程画成了SVG图，方便你理解数据是怎么一步步被“加工”好的：

从SVG图里你能看到，整个流程是线性的。从原始SAM出发，每次处理都生成一个中间文件。最终，只有排序、去重并索引过的BAM，才能交给HaplotypeCaller。

  我的建议：
  如果你在跑大批样本，可以考虑用samtools merge先把多份BAM合并，再统一走排序去重流程。
  这样能省下不少IO开销。