由于bam文件是二進(jìn)制的壓縮文件皆看，正常情況下我們無法用普通編輯器打開查看仓坞，一般情況下我們借助samtools軟件來對文件進(jìn)行處理。但很多小伙伴一定遇到過這樣的問題：如果我想對bam文件進(jìn)行個性化的處理腰吟，需要怎么做呢无埃？常規(guī)的做法是利用samtools將bam文件先轉(zhuǎn)化成可讀寫的sam文件，再進(jìn)行后續(xù)的處理蝎困。這樣做一方面犧牲了大量的時間和存儲，另一方面在編程時需要調(diào)用shell命令十分繁瑣倍啥，今天介紹的pysam可以說是一條龍服務(wù)完美解決問題禾乘。
01 安裝pysam
pysam是python的一個庫，安裝好python以后可以利用python自帶的pip安裝pysam虽缕。在后續(xù)使用中請注意始藕，samtools軟件的坐標(biāo)是1-base，而pysam對于基因組坐標(biāo)的處理是0-base：
pip install pysam
02 讀入bam文件
pysam可讀入的文件格式有bam氮趋、sam伍派、cram，以test前綴的文件作為示例剩胁，讀入方法如下诉植，讀入后生成對象bf，后續(xù)對文件的操作都基于對bf的訪問：
import pysam
bf = pysam.AlignmentFile("test.bam", "rb") #對bam文件的讀入
bf = pysam.AlignmentFile("test.sam", "r") #對sam文件的讀入
bf = pysam.AlignmentFile("test.cram", "rc") #對cram文件的讀入
03
檢查bam文件是否存在index：
bf.check_index() #返回True則存在索引昵观，返回False則索引不存在
04
pysam封裝了samtools晾腔，因此可以避開命令行調(diào)用直接對bam文件進(jìn)行sort：
pysam.sort("-o", "output.bam", "test.bam")
05
想要提取bam文件的頭文件，可以這樣操作：
head=bf.header #返回一個迭代器啊犬，可以利用循環(huán)來訪問
06
想要提取bam文件的前10行比對結(jié)果灼擂，可以這樣操作：
out=bf.head(10) #該命令會以迭代器的形式返回bam文件的前10行
07
如果想提取基因組某個區(qū)域的比對結(jié)果，比如現(xiàn)在我想提取bam文件中比對到血紅蛋白基因上的結(jié)果觉至，可以這樣操作：
hemoglobin_region=bf.fetch(contig="chr11",start=5246694,stop=5250625) #結(jié)果同樣是以迭代器的形式返回
08
想要計(jì)算落在某個區(qū)間的reads數(shù)目剔应，可以這樣操作：
readsNumber=bf.count(contig="chr1",start=1,stop=1000000) #返回reads數(shù)目值
09
想要計(jì)算某個區(qū)間的每個堿基上的覆蓋深度，可以這樣操作：
baseCov=bf.count_coverage(contig="chr1",start=1,stop=1000000) #返回一個數(shù)組，數(shù)組的長度和區(qū)間長度一致峻贮，表示每個堿基的覆蓋深度
10
對參考基因組的信息進(jìn)行統(tǒng)計(jì)席怪，可以這樣操作：
bf.get_index_statistics() #返回一個元組，每個元素記錄一條染色體比對上的reads數(shù)目和未比對上的reads數(shù)目
bf.lengths #返回一個元組月洛，每個元素記錄一條染色體長度
bf.references #以元組的形式返回參考基因組每條染色體的名稱
bf.nreferences #返回參考基因組的染色體條數(shù)
11
對比對信息進(jìn)行統(tǒng)計(jì)何恶，可以這樣操作：
bf.mapped #返回比對上的reads數(shù)目
bf.unmapped #返回未比對上的reads數(shù)目
12
對于bf對象可以利用for循環(huán)進(jìn)行訪問，每個循環(huán)是一個read嚼黔，對于read的各種處理操作匯總?cè)缦拢?/p>
for read1 in bf:
read2=bf.mate(read1) #獲取與read1配對的read2的比對信息
read1.cigarstring #字符串形式返回read1的cigar值
read1.is_duplicate #返回布爾型變量细层，判斷read1是否為PCR重復(fù)
read1.is_paired #返回布爾型變量，判斷read1是否為成對reads
read1.is_read1 #返回布爾型變量唬涧，判斷read1是否為左端reads
read1.is_read2 #返回布爾型變量疫赎，判斷read1是否為右端reads
read1.is_reverse #返回布爾型變量，判斷read1是否為反向比對
read1.is_secondary #返回布爾型變量碎节，判斷read1是否為次級比對
read1.is_unmapped #返回布爾型變量捧搞，判斷read1是否比對上參考基因組
read1.is_qcfail #返回布爾型變量，判斷read1是否qc質(zhì)控失敗
read1.mapping_quality #返回?cái)?shù)值型變量狮荔，代表read1的比對質(zhì)量
read1.mate_is_reverse #返回布爾型變量胎撇，判斷read1的配對read2是否為反向比對
read1.get_blocks() #返回?cái)?shù)值型元組，表示read1在參考基因組比對的起始坐標(biāo)和終止坐標(biāo)
read1.get_overlap(5246694,5250625) #返回?cái)?shù)值殖氏，表示read1在基因組上的坐標(biāo)和輸入?yún)?shù)坐標(biāo)的重疊區(qū)間長度
read1.get_reference_sequence() #字符串形式返回read1所在參考基因組位置的參考序列
小結(jié)：
pysam庫的使用可以大大減少在腳本中對bash命令行的調(diào)用晚树，一方面美化代碼便于后期維護(hù)更新，另一方面可以個性化調(diào)取reads信息雅采。本篇代碼涵蓋了pysam的常用方法爵憎，其他參數(shù)詳解可見官網(wǎng)https://pysam.readthedocs.io/en/latest/usage.html。
