Pysam 處理bam文件

Pysam可用來處理bam文件

安裝：

用 pip 或者 conda即可

 

使用：

Pysam的函數有很多，主要的讀取函數有：

• AlignmentFile：讀取BAM/CRAM/SAM文件

• VariantFile：讀取變異數據（VCF或者BCF）

• TabixFile：讀取由tabix索引的文件；

• FastaFile：讀取fasta序列文件；

• FastqFile：讀取fastq測序序列文件

一般常用的是第一個和第二個。

 

例子：

import pysam

bf = pysam.AlignmentFile("in.bam","rb");  其中r = read， b：binary.  二進制文件。   bam文件index

bf是一個迭代器，可以next（）或者for讀取

 for i in bf:
    print i.reference_name,i.pos,i.mapq,i.isize

結果：

ctg000331_np121 144935 27 -284
ctg000331_np121 144940 48 291
ctg000331_np121 144941 48 309
ctg000331_np121 144944 48 255
ctg000331_np121 144946 27 -370
ctg000331_np121 144947 27 -346

• i.reference_name代表read比對到的參考序列染色體id；

• i.flag bam的flag值

• i.pos代表read比對的位置；

• i.mapq代表read的比對質量值；

• i.isize代表PE read直接的插入片段長度，有時也稱Fragment長度；

 很多功能見下圖：

 

** pysam中的坐標位點是0開始，染色體起始位置為0，不是1

## sam 文件依次對應的12列
r.qname:  reads 名
r.flag ：Flag
r.reference_name: 比對到的染色體
r.pos+1： 比對位置，必須得加一
r.mapq： 比對質量
r.cigarstring： CIGAR
r.next_reference_name：另外一條reads比對的參考基因組，若和第一條相同，則輸出=
r.mpos+1： 比對的位置，必須得加1
r.isize： 插入片段長度
r.seq：reads seq
r.qual： reads 質量

 

 

 

一些功能：

• check_index() 

          檢測index文件是否存在存在即為true

bf.check_index()
True

• close（）

          用完記得關閉

bf.close()

• count(self，contig=None, start=None, stop=None, region=None, until_eof=False, read_callback='nofilter', reference=None，end=None)

             計算目標區域內比對上的reads數目

bf.count(contig="ctg000331_np121", start=1, stop=6000)
24

• count_coverage(self, contig=None, start=None, stop=None, region=None, quality_threshold=15, read_callback='all', reference=None, end=None)

         計算目標區域內的覆蓋度。返回1個4維的array，代表ACGT的覆蓋度，而每個維度的array長度為100，里面的數字代表該鹼基在各個位置上的覆蓋度。

  bf.count_coverage(contig="ctg000331_np121",start=1,stop=100)

• fetch(self, contig=None, start=None, stop=None, region=None, tid=None, until_eof=False, multiple_iterators=False, reference=None, end=None)

          提取出比對到目標區域內的全部reads。返回的是一個迭代器，可以通過for循環或者next函數從中取出reads，我們使用next()函數取出第一條reads，reads是用         AlignedSegment對象表示，可以通過該對象的內置方法再對這條reads進行一些查詢操作。

allreads=bf.fetch(contig="ctg000331_np121",start=1,stop=10000)
是一個迭代器，可以用for循環獲得

• get_index_statistics(self)
  通過index統計該BAM文件中在各個染色體上mapped/unmapped的reads個數

bf.get_index_statistics()

• fetch函數定位特定區域

有時候我們並不需要遍歷整一份BAM文件，我們可能只想獲得區中的某一個區域（比如chr1中301-310中的信息），那么這個時候可以用Alignmen模塊中的fetch函數：

bam文件必須要index

for r in bf.fetch('chr1', 300, 310)：  
    print r
bf.close()

 

 

 

關注下方公眾號可獲得更多精彩

 

 

參考

1、如何使用Pysam處理BAM

2、使用Pysam操作BAM文件