python处理bam/sam文件利器pysam

在python中读取、处理文件可以用pysam这个包。以下简单介绍一下这个包的使用。

1. 读取文件

import pysam

samfile = pysam.AlignmentFile("ENCFF191HCE.sort.bam", "rb")

仅读取某条染色体某个区域的reads：

# 这里bam文件必须先index
for read in samfile.fetch('chr1', 904920, 904930):
    print(read)

这里返回了符合的read：
HWI-ST1293:36:D11AHACXX:6:2110:12478:90259 0 0 904928 37 20M -1 -1 20 TGAGCGTCGCCGGACTCCGC array('B', [34, 34, 31, 37, 37, 37, 35, 37, 39, 39, 39, 39, 39, 40, 41, 41, 41, 41, 41, 41]) [('X0', 1), ('X1', 0), ('MD', '20'), ('PG', 'MarkDuplicates.2'), ('XG', 0), ('NM', 0), ('XM', 0), ('XO', 0), ('XT', 'U')]

2. 写入文件

# 将双端reads写入新文件
pairedreads = pysam.AlignmentFile("allpaired.bam", "wb", template=samfile)
for read in samfile.fetch():
    if read.is_paired:
        pairedreads.write(read)

pairedreads.close()
samfile.close()

3. sort

pysam.sort("-o", "output.bam", "ENCFF191HCE.sort.bam")
# 相当于
samtools sort -o output.bam ENCFF191HCE.sort.bam

4. 判断read

# 是否双端
read.is_paired

#  是否没有map上
read.is_unmapped

# 是否read1/read2
read.is_read1
read.is_read2

此外还有read.is_duplicate；read.is_proper_pair；read.is_reverse；read.is_supplementary；read.isize；read.is_qcfail；read.is_secondary。

5. flag

read.flag

此外返回0，即单端的forward read且map上了。

6. 得到标签

read.get_tag('NM') 
read.get_tag('MD')

会返回NM标签0（即没有任何mismatch和indel），MD标签的20。

7. 染色体名称

read.reference_name

返回”chr1“。

此外，还有很多别的功能，并且还可以读取操作VCF/BCF文件。详细可参阅手册：https://buildmedia.readthedocs.org/media/pdf/pysam/v0.11.2.2/pysam.pdf

欢迎关注公众号！

生信编程日常