pysam 模塊介紹！！！！

http://pysam.readthedocs.io/en/latest/index.html

?

?

在開發(fā)基因組相關(guān)流程或工具時，經(jīng)常需要讀取、處理和創(chuàng)建bam、vcf、bcf文件。目前已經(jīng)有一些主流的處理此類格式文件的工具，如samtools、picard、vcftools、bcftools，但此類工具集成的大多是標(biāo)準(zhǔn)功能，在編程時如果直接調(diào)用的話往往顯得不夠靈活。

本文介紹的是一個處理基因組數(shù)據(jù)的python模塊，它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能，能在編程時非常靈活的處理bam和bcf文件。

以下主要介紹pysam的安裝和使用方法：

1. 安裝

如果Linux上安裝了pip，可以一鍵安裝，在集群上的話，需要登錄安裝節(jié)點(diǎn)進(jìn)行安裝。

pip3 install pysam

檢查是否安裝成功

import pysam

?

2.讀取bam文件(pysam.AlignmentFile)

bam是sam的二進(jìn)制文件，因其占用空間少，所以都會使用bam進(jìn)行存儲和操作。

要讀取bam文件,必須先創(chuàng)建一個AlignmentFile對象.

path_in = ‘./test.bam‘ samfile = pysam.AlignmentFile(path_in, "rb")

之后就可以逐行讀取和處理bam文件了(順序讀取),以下打印出了bam的一行.

for line in samfile:print(line)break

?

但順序讀取還不夠靈活,我們有時需要隨機(jī)讀取(提示:sam不能隨機(jī)讀取),pysam的fetch方法提供了隨機(jī)讀取功能.

直接使用fetch會報(bào)錯

ValueError: fetch called on bamfile without index

提示我們需要建立(.bai)索引

samtools index corrected.bam

fetch返回的是一個迭代器(iterator),可以迭代讀取內(nèi)容.

for read in samfile.fetch(‘chr6‘, 28478220, 28478222): ... print(read)

fetch方法的API如下,chr6為參考序列,后面數(shù)字分別為讀取的起始和終止位置.

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

?

3.讀取vcf/bcf文件(pysam.VariantFile)

讀取方法同上,只是使用的是VariantFile方法:

gvcf = "./MHC.unified.g.vcf.gz" vcf_in = pysam.VariantFile(gvcf)

若想隨機(jī)讀取,仍然需要建立索引:

首先使用bgzip壓縮vcf

bgzip -c MHC.g.vcf > MHC.g.vcf.gz

然后用bcftools建立索引

bcftools index -c MHC.g.vcf.gz

使用fetch讀取

for rec in vcf_in.fetch(‘chr6‘, 28577796, 28577896): ... print(rec) ... break

?

4.創(chuàng)建并寫入到新的bam或vcf文件

pysam的核心功能是可以隨心所欲的讀取數(shù)據(jù),處理之后,寫入到一個新建的bam或bcf文件里.

我們可以完全自定義一些內(nèi)容,然后寫入到一個新的bam文件里,如下:

header = { ‘HD‘: {‘VN‘: ‘1.0‘},‘SQ‘: [{‘LN‘: 1575, ‘SN‘: ‘chr1‘},{‘LN‘: 1584, ‘SN‘: ‘chr2‘}] }with pysam.AlignmentFile(tmpfilename, "wb", header=header) as outf:a = pysam.AlignedSegment()a.query_name = "read_28833_29006_6945"a.query_sequence="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"a.flag = 99a.reference_id = 0a.reference_start = 32a.mapping_quality = 20a.cigar = ((0,10), (2,1), (0,25))a.next_reference_id = 0a.next_reference_start=199a.template_length=167a.query_qualities = pysam.qualitystring_to_array("<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<")a.tags = (("NM", 1),("RG", "L1"))outf.write(a)

同理,我們也可以讀取一個已有的bam文件,逐個修改以上的屬性,然后存儲到一個新的bam文件里.這里不再舉例.

上面設(shè)置header可能有點(diǎn)麻煩,容易出錯,但我們可以復(fù)制一個已有bam文件的header到一個新的bam文件里.

outf = pysam.AlignmentFile(path_out, "wb", template=samfile)

以上template參數(shù)指定了模板bam文件.

?

5. 關(guān)閉文件

outf.close()

?

總結(jié):

pysam模塊非常實(shí)用,有了pysam模塊,我們就可以非常靈活的操縱bam/bcf文件,而不必依賴于samtools或bcftools. pysam可以隨機(jī)讀取bam/bcf文件,也可以將處理后的內(nèi)容自定義輸出到bam/bcf文件.

?

以上只介紹了pysam最常見的功能,更多pysam功能請參照:http://pysam.readthedocs.io/en/latest/index.html

pysam - 多種格式基因組數(shù)據(jù)（sam/bam/vcf/bcf/cram/…）讀寫與處理模塊（python）

標(biāo)簽：class???style???log???http???it???使用???la???sp???文件???

原文：http://www.cnblogs.com/leezx/p/5908767.html

轉(zhuǎn)載于:https://www.cnblogs.com/nkwy2012/p/6558069.html

總結(jié)

以上是生活随笔為你收集整理的pysam - 多种格式基因组数据（sam/bam/vcf/bcf/cram/…）读写与处理模块（python)--转载...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。