在上一篇的文章裡我詳細介紹了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的方法，希望對大家有所幫助，如果大家有任何疑問請給我留言，小編會及時回覆大家的。在此也非常感謝大家對指令碼之家網站的支援！