搜索
查看: 2524|回复: 3

[Python] pysam - 多种数据格式读写与处理模块(python)

[复制链接]

2

主题

8

帖子

284

积分

信息监察员

人见人爱的

Rank: 9Rank: 9Rank: 9

积分
284
发表于 2016-9-26 13:07:13 | 显示全部楼层 |阅读模式
在开发基因组相关流程或工具时,经常需要读取、处理和创建bam、vcf、bcf文件。目前已经有一些主流的处理此类格式文件的工具,如samtools、picard、vcftools、bcftools,但此类工具集成的大多是标准功能,在编程时如果直接调用的话往往显得不够灵活。
本文介绍的是一个处理基因组数据的python模块,它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能,能在编程时非常灵活的处理bam和bcf文件。
以下主要介绍pysam的安装和使用方法:

1. 安装
如果Linux上安装了pip,可以一键安装,在集群上的话,需要登录安装节点进行安装。
[Python] 纯文本查看 复制代码
pip3 install pysam

检查是否安装成功
[Python] 纯文本查看 复制代码
import pysam


2.读取bam文件(pysam.AlignmentFile)
bam是sam的二进制文件,因其占用空间少,所以都会使用bam进行存储和操作。
要读取bam文件,必须先创建一个AlignmentFile对象.
[Python] 纯文本查看 复制代码
path_in = './test.bam'
samfile = pysam.AlignmentFile(path_in, "rb")

之后就可以逐行读取和处理bam文件了(顺序读取),以下打印出了bam的一行.
[Python] 纯文本查看 复制代码
for line in samfile:
    print(line)
    break

但顺序读取还不够灵活,我们有时需要随机读取(提示:sam不能随机读取),pysam的fetch方法提供了随机读取功能.
直接使用fetch会报错
[Bash shell] 纯文本查看 复制代码
ValueError: fetch called on bamfile without index

提示我们需要建立(.bai)索引
[Bash shell] 纯文本查看 复制代码
samtools index corrected.bam

fetch返回的是一个迭代器(iterator),可以迭代读取内容.
[Python] 纯文本查看 复制代码
for read in samfile.fetch('chr6', 28478220, 28478222):
...     print(read)

fetch方法的API如下,chr6为参考序列,后面数字分别为读取的起始和终止位置.
[Python] 纯文本查看 复制代码
fetch(self, reference=None, start=None, end=None, region=None, tid=None, until_eof=False, multiple_iterators=False)


3.读取vcf/bcf文件(pysam.VariantFile)
读取方法同上,只是使用的是VariantFile方法:
[Python] 纯文本查看 复制代码
gvcf = "./MHC.unified.g.vcf.gz"
vcf_in = pysam.VariantFile(gvcf)

若想随机读取,仍然需要建立索引:
首先使用bgzip压缩vcf
[Bash shell] 纯文本查看 复制代码
bgzip -c MHC.g.vcf > MHC.g.vcf.gz

然后用bcftools建立索引
[Bash shell] 纯文本查看 复制代码
bcftools index -c MHC.g.vcf.gz

使用fetch读取
[Python] 纯文本查看 复制代码
for rec in vcf_in.fetch('chr6', 28577796, 28577896):
...     print(rec)
...     break


4.创建并写入到新的bam或vcf文件
pysam的核心功能是可以随心所欲的读取数据,处理之后,写入到一个新建的bam或bcf文件里.
我们完全可以自定义一些内容,然后写入到一个新的bam文件里,如下:
[Python] 纯文本查看 复制代码
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文件里.
[Python] 纯文本查看 复制代码
outf = pysam.AlignmentFile(path_out, "wb", template=samfile)

以上template参数指定了模板bam文件.

5.读取fasta/fastq文件(faidx建立索引)
读取方法略有不同,fetch返回的本身就是一个字符串。
[Bash shell] 纯文本查看 复制代码
samtools faidx total_PacBio_reads.fasta

[Python] 纯文本查看 复制代码
fasta_file = pysam.FastaFile(path)
fasta_file.fetch("m160727_060737_42266_c101014182550000001823222610211695_s1_p0/110008/22268_22731")


6. 关闭文件
[Python] 纯文本查看 复制代码
outf.close()


总结:
pysam模块非常实用,有了pysam模块,我们就可以非常灵活的操纵bam/bcf文件,而不必依赖于samtools或bcftools. pysam可以随机读取bam/bcf文件,也可以将处理后的内容自定义输出到bam/bcf文件.

以上介绍了pysam最常见的功能,更多pysam功能请参照:http://pysam.readthedocs.io/en/latest/index.html




上一篇:【月亮-Perl练习题2】使用Perl定长切割字符串
下一篇:[ROSALIND上的题]考虑兔子死亡的兔子繁殖问题
回复

使用道具 举报

7

主题

44

帖子

186

积分

版主

Rank: 7Rank: 7Rank: 7

积分
186
发表于 2016-9-26 23:08:33 | 显示全部楼层
好详细!感谢楼主
微博:bioinfoHeTing 欢迎来访 ^ ^
回复 支持 反对

使用道具 举报

633

主题

1172

帖子

3947

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3947
发表于 2016-9-27 07:32:25 来自手机 | 显示全部楼层
我比较纠结这个板块是有个放在这个板块还是计算机基础板块
回复 支持 反对

使用道具 举报

0

主题

6

帖子

83

积分

注册会员

Rank: 2

积分
83
发表于 2017-12-11 19:04:51 | 显示全部楼层
我一直装不上这个包,大神指点一下
在Windows上使用pip install pysam一直报错,这个包可以再win上装吗
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2017-12-19 06:01 , Processed in 0.121266 second(s), 32 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.