在上一篇的文章里我詳細介紹了BAM（SAM/CRAM）的格式和一些需要注意的細節(jié)，還說了該如何使用samtools在命令行中對其進行操作。但是很多時候這些操作是不能滿足我們的實際需要的，比如統(tǒng)計比對率、計算在某個比對質量值之上的read有多少，或者計算PE比對的插入片段長度分布，甚至需要你根據(jù)實際情況編寫一個新的變異檢測算法等。這個時候往往難以直接通過samtools來實現(xiàn)【注】，而是需要編寫專門的程序進行計算。因此，在這一篇文章里我們就一起來學習應該如何在程序中借助Pysam來處理BAM文件。

【注】關于統(tǒng)計比對率其實是可以通過samtools stats計算獲得的。不過我們這篇文章不是為了爭辯samtools能做什么，不能做什么，而是要跟大家討論該如何編寫程序處理BAM。

不過，在開始之前我想稍微再補充一下上一節(jié)中提到的CRAM——我習慣將其稱為BAM的高壓縮格式，因為它和BAM/SAM的格式基本相同，但有四點我們需要注意一下：

CRAM的高壓縮是通過借助參考序列和對其他信息的進一步編碼來實現(xiàn)的，它相比于BAM有著更高的壓縮率，能夠節(jié)省30%-50%的空間；

CRAM目前的IO效率沒有BAM高（壓得密嘛），約慢30%，但在不斷進步，現(xiàn)在已經更新到了3.x版本了；
CRAM和BAM可以通過samtools或者picard方便地實現(xiàn)互轉；
CRAM一定會取代BAM，這話并不是我說的，而是bwa/samtools的作者lh3說的。

什么是Pysam

Pysam是一個專門用來處理（BAM/CRAM/SAM）比對數(shù)據(jù)和變異數(shù)據(jù)（VCF和BCF）的Python包。它的核心是htslib——一個高通量數(shù)據(jù)處理API（來自samtools和bwa的核心，基于C語言），開發(fā)者們用Python對它直接進行輕量級包裝，因此能夠在Python中方便地進行調用，并且保證了它與原生C-API功能上的高度一致。

為什么是Pysam

因為Pysam可以說是最為官方的版本，有比較固定的開發(fā)者在維護，它的穩(wěn)定性和可靠性都很高。雖然還有一些其它的包同樣能夠處理BAM但其實它們大多繞不開對htslib的使用，但卻沒有pysam周全。而且Pysam還集成了tabix的接口，所以除了比對數(shù)據(jù)之外，還能夠用于處理所有用tabix構建過索引的文件，總之就是全且可靠。

如果是文本格式的sam的話，其實也可以直接將其當作普通文本文件來處理，不需借助任何程序包（這在早期的數(shù)據(jù)分析中經?？吹竭@種操作），只是要麻煩很多（必須自己在程序中處理所有細節(jié)，包括解析FLAG和CIGAR信息，以前我也干過不少類似的事情），甚至我還看到有人直接在程序中調用samtools view把BAM轉換成SAM之后再處理的。。。這樣的做法實在不推薦。

所以，只要你用的是Python，那么Pysam真的是目前看來比較好的選擇。當然如果你用C/C++那么直接用htslib或者bamtools，如果是Java，那么直接使用htsjdk——htslib的java版本。