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

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

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

CRAM的高壓縮是通過借助參考序列和對(duì)其他信息的進(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）比對(duì)數(shù)據(jù)和變異數(shù)據(jù)（VCF和BCF）的Python包。它的核心是htslib——一個(gè)高通量數(shù)據(jù)處理API（來自samtools和bwa的核心，基于C語言），開發(fā)者們用Python對(duì)它直接進(jìn)行輕量級(jí)包裝，因此能夠在Python中方便地進(jìn)行調(diào)用，并且保證了它與原生C-API功能上的高度一致。

為什么是Pysam

因?yàn)镻ysam可以說是最為官方的版本，有比較固定的開發(fā)者在維護(hù)，它的穩(wěn)定性和可靠性都很高。雖然還有一些其它的包同樣能夠處理BAM但其實(shí)它們大多繞不開對(duì)htslib的使用，但卻沒有pysam周全。而且Pysam還集成了tabix的接口，所以除了比對(duì)數(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版本。

如何使用Pysam

首先，要為我們的Python環(huán)境安裝這個(gè)包，如果已安裝過的話可以忽略這一步。有兩個(gè)方法，pip和bioconda，都比較簡單，我們這里以pip——Python的包管理工具來進(jìn)行：

$ pip install pysam

安裝完成之后我們就可以在Python程序中調(diào)用pysam了。

讀取BAM/CRAM/SAM文件

Pysam中的函數(shù)有很多，但是重要的讀取函數(shù)主要有：

• AlignmentFile：讀取BAM/CRAM/SAM文件
• VariantFile：讀取變異數(shù)據(jù)（VCF或者BCF）
• TabixFile：讀取由tabix索引的文件；
• FastaFile：讀取fasta序列文件；
• FastqFile：讀取fastq測序序列文件；

等以上幾個(gè)，其中尤以AlignmentFile和VariantFile為核心。需要我們注意到的地方是，Pysam中的有些函數(shù)由于歷史原因存在重復(fù)，比如名字上只有大小寫的差異，但功能卻是一樣的（比如下圖的TabixFile），有些則只是簡化了函數(shù)名，這些情況用的時(shí)候留個(gè)心眼就行了。

另外，這篇文章的目的是介紹如何處理比對(duì)文件，所以我打算只介紹AlignmentFile。

讀取比對(duì)文件前，我建議先使用samtools index為比對(duì)文件構(gòu)建好索引。當(dāng)然如果是SAM文件就不必了——它是文本文件，索引的作用是讓我們可以對(duì)文件進(jìn)行隨機(jī)讀取，而不必總是從頭開始。

下面我先用一個(gè)例子作為引子：

import pysam
bf = pysam.AlignmentFile('in.bam', 'rb')

我在這個(gè)例子里面，先在程序中導(dǎo)入pysam包，然后調(diào)用AlignmentFile函數(shù)讀取'in.bam'文件，并把句柄賦值給了bf，bf其實(shí)是一個(gè)迭代器——Python中的術(shù)語，意思就是適合在for循環(huán)中進(jìn)行遍歷的對(duì)象。

這樣我們就是可以通過bf獲取這份比對(duì)文件中的內(nèi)容了。比如我們想把in.bam中每一條read的比對(duì)位置（包含染色體編號(hào)和位置信息），比對(duì)質(zhì)量值和插入片段長度輸出（我們的in.bam來自PE測序數(shù)據(jù)的結(jié)果），那么可以這樣做：

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循環(huán)對(duì)其從頭到尾地遍歷，并把每個(gè)值都賦給r，r在這里代表的就是比對(duì)的read信息，它是一個(gè)對(duì)象（在Pysam由AlignedSegment定義），通過它就可以獲取所有的比對(duì)信息，比如上面例子中：

• r.reference_name代表read比對(duì)到的參考序列染色體id；
• r.pos代表read比對(duì)的位置；
• r.mapq代表read的比對(duì)質(zhì)量值；
• r.isize代表PE read直接的插入片段長度，有時(shí)也稱Fragment長度；

這里例子的結(jié)果如下：

chrM 160 50 235
chrM 161 30 -283
chrM 314 60 -207
...

另外，由于bf是一個(gè)迭代器，我們其實(shí)還可以用bf.next()一個(gè)一個(gè)地對(duì)其進(jìn)行訪問，而不必在for循環(huán)中遍歷，這在一些特殊的情況下，這個(gè)做法是非常有用且方便的。

當(dāng)然，上面這個(gè)例子其實(shí)非常簡單，實(shí)際上r變量中還有很多其它關(guān)于比對(duì)的信息，下面這個(gè)截圖，就是變量中能夠獲取到的所有比對(duì)相關(guān)的信息，有好幾十個(gè)。

眼尖的同學(xué)可能也發(fā)現(xiàn)了，這里面存在一些名字類似的變量，如：r.mapping_quality 和 r.mapq，它們其實(shí)都是比對(duì)質(zhì)量值。類似的也還有幾個(gè)，這都是上面我提到的歷史原因所致，不過這種多余變量隨著Pysam的維護(hù)也正在逐步變少。

此外，Pysam中的位點(diǎn)坐標(biāo)體系是0-base（意思是染色體的起始位置是從0而不是1開始算的）而不是1-base，所以上面的輸出的160，其實(shí)真實(shí)位置應(yīng)該要+1，也就是161。

還有，上文我也說過，AlignmentFile除了能夠讀/寫B(tài)AM之外，還同樣能夠讀/寫CRAM和SAM。區(qū)別就在于函數(shù)中的第二個(gè)參數(shù)，比如上面例子中的字符'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')

讀取特定比對(duì)區(qū)域內(nèi)的數(shù)據(jù)

有時(shí)候我們并不需要遍歷整一份BAM文件，我們可能只想獲得區(qū)中的某一個(gè)區(qū)域（比如chrM中301-310中的信息），那么這個(gè)時(shí)候可以用Alignmen模塊中的fetch函數(shù)：

import pysam
bf = AlignmentFile('in.bam', 'rb')
for r in bf.fetch('chrM', 300, 310)：
 print r
bf.close()

通過fetch函數(shù)就可以定位特定區(qū)域了，非常方便。不過，這個(gè)時(shí)候輸入文件in.bam就必須要有索引，不然無法實(shí)現(xiàn)這種讀取操作。最后用完了，要記得關(guān)閉文件（bf.close()）。

來個(gè)稍微難一點(diǎn)的例子

問題：如何輸出覆蓋在某個(gè)位置上，比對(duì)質(zhì)量值大于30的所有堿基？

這個(gè)問題包含兩個(gè)部分：

• 固定的某個(gè)位置（我們這里還是用chrM 301這個(gè)位置）
• read比對(duì)質(zhì)量值必須是大于30。

如何做呢？這個(gè)時(shí)候我們要用AlignmentFile模塊的另一個(gè)函數(shù)——pileups來協(xié)助解決，代碼如下：

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()

這段代碼看起來雖然簡單，但其實(shí)包含了很多信息?？偟膩碚f，就是通過pileup獲取了所有覆蓋到該位置的read，并將其存到pileupcolumn中。然后，對(duì)pileupcolumn調(diào)用pileups（注意多了一個(gè)s）獲得一條read中每個(gè)比對(duì)位置的信息（一條read那么長，并非只覆蓋了一個(gè)位置），然后通過判斷語句留下覆蓋到目標(biāo)位點(diǎn)（301）的堿基。代碼中的read.alignment是Pysam中AlignedSegment對(duì)象，它包含的內(nèi)容和上述其它例子中的r是一樣的。read.alignment.pos + 1還是0-base的原因。最后結(jié)果如下：

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
...

創(chuàng)建BAM/CRAM/SAM文件

最后這個(gè)例子，我想告訴大家該如何用Pysam輸出BAM/CRAM/SAM格式，具體還是看代碼吧，這里想輸出結(jié)果是BAM文件，所以輸出模式是“wb”，例子中我們只輸出一條比對(duì)結(jié)果作為說明。

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() # 定義一個(gè)AlignedSegment對(duì)象用于存儲(chǔ)比對(duì)信息
  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)

小結(jié)

我寫這篇文章的目的主要有兩個(gè)：第一，充實(shí)上一篇文章中關(guān)于如何操作BAM的內(nèi)容；第二，介紹Pysam這一個(gè)值得使用的包給大家。另外，我上面列舉的例子其實(shí)都比較偏于基礎(chǔ)操作，這可能和我自身對(duì)認(rèn)知的看法有關(guān)。我一直認(rèn)為，只有真正理解并靈活地應(yīng)用基礎(chǔ)操作，才可以靈活地解決一切復(fù)雜的問題。

而且，上面幾個(gè)例子中用到的模塊和函數(shù)其實(shí)都是比較常用的，所以我比較推薦優(yōu)先掌握它們。這些例子里面用到的數(shù)據(jù)我已放到github了，感興趣的同學(xué)可以在公眾號(hào)后臺(tái)回復(fù)“WGS”即可獲得，后續(xù)也會(huì)陸續(xù)有其它的代碼和數(shù)據(jù)可供參考。

最后，Pysam的內(nèi)容其實(shí)還有很多，我所介紹的也僅在對(duì)于比對(duì)數(shù)據(jù)的處理，其它很多的模塊和函數(shù)，包括對(duì)Fasta，F(xiàn)astq，VCF，BCF和Tabix文件的處理，我就不進(jìn)行一一介紹了，建議大家在使用的時(shí)候多看看它的完整文檔。

以上所述是小編給大家介紹的使用Python處理BAM的方法，希望對(duì)大家有所幫助，如果大家有任何疑問請(qǐng)給我留言，小編會(huì)及時(shí)回復(fù)大家的。在此也非常感謝大家對(duì)腳本之家網(wǎng)站的支持！

最新評(píng)論