在開發基因組相關流程或工具時，經常需要讀取、處理和創建bam、vcf、bcf文件。目前已經有一些主流的處理此類格式文件的工具，如samtools、picard、vcftools、bcftools，但此類工具集成的大多是標准功能，在編程時如果直接調用的話往往顯得不夠靈活。
本文介紹的是一個處理基因組數據的python模塊，它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能，能在編程時非常靈活的處理bam和bcf文件。
以下主要介紹pysam的安裝和使用方法：
1. 安裝
如果Linux上安裝了pip，可以一鍵安裝，在集群上的話，需要登錄安裝節點進行安裝。
pip3 install pysam
檢查是否安裝成功
import pysam
2.讀取bam文件(pysam.AlignmentFile)
bam是sam的二進制文件，因其占用空間少，所以都會使用bam進行存儲和操作。
要讀取bam文件,必須先創建一個AlignmentFile對象.
path_in = './test.bam' samfile = pysam.AlignmentFile(path_in, "rb")
之后就可以逐行讀取和處理bam文件了(順序讀取),以下打印出了bam的一行.
for line in samfile: print(line) break
但順序讀取還不夠靈活,我們有時需要隨機讀取(提示:sam不能隨機讀取),pysam的fetch方法提供了隨機讀取功能.
直接使用fetch會報錯
ValueError: fetch called on bamfile without index
提示我們需要建立(.bai)索引
samtools index corrected.bam
fetch返回的是一個迭代器(iterator),可以迭代讀取內容.
for read in samfile.fetch('chr6', 28478220, 28478222): ... print(read)
fetch方法的API如下,chr6為參考序列,后面數字分別為讀取的起始和終止位置.
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)
若想隨機讀取,仍然需要建立索引:
首先使用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
3.隨機讀取fasta文件（faidx建立索引）
讀取方法略有不同，fetch返回的本身就是一個字符串。
samtools faidx total_PacBio_reads.fasta
fasta_file = pysam.FastaFile(path) fasta_file.fetch("m160727_060737_42266_c101014182550000001823222610211695_s1_p0/110008/22268_22731")
4.創建並寫入到新的bam或vcf文件
pysam的核心功能是可以隨心所欲的讀取數據,處理之后,寫入到一個新建的bam或bcf文件里.
我們可以完全自定義一些內容,然后寫入到一個新的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 = 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文件,逐個修改以上的屬性,然后存儲到一個新的bam文件里.這里不再舉例.
上面設置header可能有點麻煩,容易出錯,但我們可以復制一個已有bam文件的header到一個新的bam文件里.
outf = pysam.AlignmentFile(path_out, "wb", template=samfile)
以上template參數指定了模板bam文件.
5. 關閉文件
outf.close()
總結:
pysam模塊非常實用,有了pysam模塊,我們就可以非常靈活的操縱bam/bcf文件,而不必依賴於samtools或bcftools. pysam可以隨機讀取bam/bcf文件,也可以將處理后的內容自定義輸出到bam/bcf文件.
以上只介紹了pysam最常見的功能,更多pysam功能請參照:http://pysam.readthedocs.io/en/latest/index.html