在上一篇的文章里我詳細介紹了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版本。
如何使用Pysam
首先，要為我們的Python環境安裝這個包，如果已安裝過的話可以忽略這一步。有兩個方法，pip和bioconda，都比較簡單，我們這里以pip――Python的包管理工具來進行：
$ pip install pysam
安裝完成之后我們就可以在Python程序中調用pysam了。
讀取BAM/CRAM/SAM文件
Pysam中的函數有很多，但是重要的讀取函數主要有：
- AlignmentFile：讀取BAM/CRAM/SAM文件
- VariantFile：讀取變異數據（VCF或者BCF）
- TabixFile：讀取由tabix索引的文件；
- FastaFile：讀取fasta序列文件；
- FastqFile：讀取fastq測序序列文件；
等以上幾個，其中尤以AlignmentFile和VariantFile為核心。需要我們注意到的地方是，Pysam中的有些函數由於歷史原因存在重復，比如名字上只有大小寫的差異，但功能卻是一樣的（比如下圖的TabixFile），有些則只是簡化了函數名，這些情況用的時候留個心眼就行了。
另外，這篇文章的目的是介紹如何處理比對文件，所以我打算只介紹AlignmentFile。
讀取比對文件前，我建議先使用samtools
index為比對文件構建好索引。當然如果是SAM文件就不必了――它是文本文件，索引的作用是讓我們可以對文件進行隨機讀取，而不必總是從頭開始。
下面我先用一個例子作為引子：
import pysam
bf = pysam.AlignmentFile('in.bam', 'rb')
我在這個例子里面，先在程序中導入pysam包，然后調用AlignmentFile函數讀取'in.bam'文件，並把句柄賦值給了bf，bf其實是一個迭代器――Python中的術語，意思就是適合在for循環中進行遍歷的對象。
這樣我們就是可以通過bf獲取這份比對文件中的內容了。比如我們想把in.bam中每一條read的比對位置（包含染色體編號和位置信息），比對質量值和插入片段長度輸出（我們的in.bam來自PE測序數據的結果），那么可以這樣做：
import pysam
bf = pysam.AlignmentFile('in.bam', 'rb')
for r in bf:
print r.reference_name, r.pos, r.mapq, r.isize
是不是很簡單！打開in.bam文件之后，用for循環對其從頭到尾地遍歷，並把每個值都賦給r，r在這里代表的就是比對的read信息，它是一個對象（在Pysam由AlignedSegment定義），通過它就可以獲取所有的比對信息，比如上面例子中：
- r.reference_name代表read比對到的參考序列染色體id；
- r.pos代表read比對的位置；
- r.mapq代表read的比對質量值；
- r.isize代表PE read直接的插入片段長度，有時也稱Fragment長度；
這里例子的結果如下：
chrM 160 50 235
chrM 161 30 -283
chrM 314 60 -207
...
另外，由於bf是一個迭代器，我們其實還可以用bf.next()一個一個地對其進行訪問，而不必在for循環中遍歷，這在一些特殊的情況下，這個做法是非常有用且方便的。
當然，上面這個例子其實非常簡單，實際上r變量中還有很多其它關於比對的信息，下面這個截圖，就是變量中能夠獲取到的所有比對相關的信息，有好幾十個。
眼尖的同學可能也發現了，這里面存在一些名字類似的變量，如：r.mapping_quality 和
r.mapq，它們其實都是比對質量值。類似的也還有幾個，這都是上面我提到的歷史原因所致，不過這種多余變量隨着Pysam的維護也正在逐步變少。
此外，Pysam中的位點坐標體系是0-base（意思是染色體的起始位置是從0而不是1開始算的）而不是1-base，所以上面的輸出的160，其實真實位置應該要+1，也就是161。
還有，上文我也說過，AlignmentFile除了能夠讀/寫BAM之外，還同樣能夠讀/寫CRAM和SAM。區別就在於函數中的第二個參數，比如上面例子中的字符'b'就是用於明確指定BAM文件，'r'字符代表“只讀”模式（read首字母）。如果要打開CRAM文件，只需要把b換成c（代表CRAM）就行了，如下：
import pysam
cf = pysam.AlignmentFile('in.cram', 'rc')
那么，如果是SAM文件呢？去掉b或c即可：
import pysam
sf = pysam.AlignmentFile('in.sam', 'r')
讀取特定比對區域內的數據
有時候我們並不需要遍歷整一份BAM文件，我們可能只想獲得區中的某一個區域（比如chrM中301-310中的信息），那么這個時候可以用Alignmen模塊中的fetch函數：
import pysam
bf = AlignmentFile('in.bam', 'rb')
for r in bf.fetch('chrM', 300, 310)：
print r
bf.close()
通過fetch函數就可以定位特定區域了，非常方便。不過，這個時候輸入文件in.bam就必須要有索引，不然無法實現這種讀取操作。最后用完了，要記得關閉文件（bf.close()）。
來個稍微難一點的例子
問題：如何輸出覆蓋在某個位置上，比對質量值大於30的所有鹼基？
這個問題包含兩個部分：
- 固定的某個位置（我們這里還是用chrM 301這個位置）
- read比對質量值必須是大於30。
如何做呢？這個時候我們要用AlignmentFile模塊的另一個函數――pileups來協助解決，代碼如下：
import pysam
bf = pysam.AlignmentFile("in.bam", "rb" )
for pileupcolumn in bf.pileup("chrM", 300, 301):
for read in [al for al in pileupcolumn.pileups if al.alignment.mapq>30]:
if not read.is_del and not read.is_refskip:
if read.alignment.pos + 1 == 301:
print read.alignment.reference_name,\
read.alignment.pos + 1,\
read.alignment.query_sequence[read.query_position]
bf.close()
這段代碼看起來雖然簡單，但其實包含了很多信息。總的來說，就是通過pileup獲取了所有覆蓋到該位置的read，並將其存到pileupcolumn中。然后，對pileupcolumn調用pileups（注意多了一個s）獲得一條read中每個比對位置的信息（一條read那么長，並非只覆蓋了一個位置），然后通過判斷語句留下覆蓋到目標位點（301）的鹼基。代碼中的read.alignment是Pysam中AlignedSegment對象，它包含的內容和上述其它例子中的r是一樣的。read.alignment.pos
- 1還是0-base的原因。最后結果如下：
chrM 301 A
chrM 301 A
chrM 301 A
chrM 301 C
chrM 301 C
chrM 301 C
chrM 301 C
chrM 301 C
chrM 301 C
chrM 301 C
...
創建BAM/CRAM/SAM文件
最后這個例子，我想告訴大家該如何用Pysam輸出BAM/CRAM/SAM格式，具體還是看代碼吧，這里想輸出結果是BAM文件，所以輸出模式是“wb”，例子中我們只輸出一條比對結果作為說明。
import pysam
header = {'HD': {'VN': '1.0'},
'SQ': [{'LN': 1575, 'SN': 'chr1'},
{'LN': 1584, 'SN': 'chr2'}]
}
tmpfilename = "out.bam"
with pysam.AlignmentFile(tmpfilename, "wb", header=header) as outf:
a = pysam.AlignedSegment() # 定義一個AlignedSegment對象用於存儲比對信息
a.query_name = "read_28833_29006_6945"
a.query_sequence="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"
a.flag = 99
a.reference_id = 0
a.reference_start = 32
a.mapping_quality = 20
a.cigar = ((0,10), (2,1), (0,25))
a.next_reference_id = 0
a.next_reference_start=199
a.template_length=167
a.query_qualities = pysam.qualitystring_to_array("<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<")
a.tags = (("NM", 1),
("RG", "L1"))
outf.write(a)
小結
我寫這篇文章的目的主要有兩個：第一，充實上一篇文章中關於如何操作BAM的內容；第二，介紹Pysam這一個值得使用的包給大家。另外，我上面列舉的例子其實都比較偏於基礎操作，這可能和我自身對認知的看法有關。我一直認為，只有真正理解並靈活地應用基礎操作，才可以靈活地解決一切復雜的問題。
而且，上面幾個例子中用到的模塊和函數其實都是比較常用的，所以我比較推薦優先掌握它們。這些例子里面用到的數據我已放到github了，感興趣的同學可以在公眾號后台回復“WGS”即可獲得，后續也會陸續有其它的代碼和數據可供參考。
最后，Pysam的內容其實還有很多，我所介紹的也僅在對於比對數據的處理，其它很多的模塊和函數，包括對Fasta，Fastq，VCF，BCF和Tabix文件的處理，我就不進行一一介紹了，建議大家在使用的時候多看看它的
完整文檔 。
以上所述是小編給大家介紹的使用Python處理BAM的方法，希望對大家有所幫助，如果大家有任何疑問請給我留言，小編會及時回復大家的。在此也非常感謝大家對腳本之家網站的支持！
