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

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

不過，在開始之前我想稍微再補(bǔ)充一下上一節(jié)中提到的CRAM——我習(xí)慣將其稱為BAM的高壓縮格式，因?yàn)樗虰AM/SAM的格式基本相同，但有四點(diǎn)我們需要注意一下：

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

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

什么是Pysam

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

為什么是Pysam

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

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

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