使用Python处理BAM的方法

在上一篇的文章里我详细介绍了BAM（SAM/CRAM）的格式和一些需要注意的细节，还说了该如何使用samtools在命令行中对其进行操作。但是很多时候这些操作是不能满足我们的实际需要的，比如统计比对率、计算在某个比对质量值之上的read有多少，或者计算PE比对的插入片段长度分布，甚至需要你根据实际情况编写一个新的变异检测算法等。这个时候往往难以直接通过samtools来实现【注】，而是需要编写专门的程序进行计算。因此，在这一篇文章里我们就一起来学习应该如何在程序中借助Pysam来处理BAM文件。

【注】关于统计比对率其实是可以通过samtools stats计算获得的。不过我们这篇文章不是为了争辩samtools能做什么，不能做什么，而是要跟大家讨论该如何编写程序处理BAM。

不过，在开始之前我想稍微再补充一下上一节中提到的CRAM——我习惯将其称为BAM的高压缩格式，因为它和BAM/SAM的格式基本相同，但有四点我们需要注意一下：

CRAM的高压缩是通过借助参考序列和对其他信息的进一步编码来实现的，它相比于BAM有着更高的压缩率，能够节省30%-50%的空间；

CRAM目前的IO效率没有BAM高（压得密嘛），约慢30%，但在不断进步，现在已经更新到了3.x版本了；
CRAM和BAM可以通过samtools或者picard方便地实现互转；
CRAM一定会取代BAM，这话并不是我说的，而是bwa/samtools的作者lh3说的。

什么是Pysam

Pysam是一个专门用来处理（BAM/CRAM/SAM）比对数据和变异数据（VCF和BCF）的Python包。它的核心是htslib——一个高通量数据处理API（来自samtools和bwa的核心，基于C语言），开发者们用Python对它直接进行轻量级包装，因此能够在Python中方便地进行调用，并且保证了它与原生C-API功能上的高度一致。

为什么是Pysam

因为Pysam可以说是最为官方的版本，有比较固定的开发者在维护，它的稳定性和可靠性都很高。虽然还有一些其它的包同样能够处理BAM但其实它们大多绕不开对htslib的使用，但却没有pysam周全。而且Pysam还集成了tabix的接口，所以除了比对数据之外，还能够用于处理所有用tabix构建过索引的文件，总之就是全且可靠。

如果是文本格式的sam的话，其实也可以直接将其当作普通文本文件来处理，不需借助任何程序包（这在早期的数据分析中经常看到这种操作），只是要麻烦很多（必须自己在程序中处理所有细节，包括解析FLAG和CIGAR信息，以前我也干过不少类似的事情），甚至我还看到有人直接在程序中调用samtools view把BAM转换成SAM之后再处理的。。。这样的做法实在不推荐。

所以，只要你用的是Python，那么Pysam真的是目前看来比较好的选择。当然如果你用C/C++那么直接用htslib或者bamtools，如果是Java，那么直接使用htsjdk——htslib的java版本。