在开发基因组相关流程或工具时，经常需要读取、处理和创建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("<<<<<<<<<<<<<<<<<<<<<:>

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 - 多种格式基因组数据(sam&sol;bam&sol;vcf&sol;bcf&sol;cram&sol;…)读写与处理模块(python&rpar;--转载

pysam 模块介绍!!!! http://pysam.readthedocs.io/en/latest/index.html 在开发基因组相关流程或工具时,经常需要读取.处理和创建bam.vcf.b ...

beego——多种格式的数据输出

beego当初设计的时候就考虑了API功能的设计,而我们在设计API的时候经常是输出JSON或者XML数据,那么beego提供了这样的方式直接输出: 1.JSON格式输出 func (this *Ad ...

Edit Distance编辑距离(NM tag)- sam&sol;bam格式解读进阶

sam格式很精炼,几乎包含了比对的所有信息,我们平常用到的信息很少,但特殊情况下,我们会用到一些较为生僻的信息,关于这些信息sam官方文档的介绍比较精简,直接看估计很难看懂. 今天要介绍的是如何通过b ...

sam&sol;bam格式

1)Sam (Sequence Alignment/Map) ------------------------------------------------- 1) SAM 文件产生背景 随着Ill ...

python多种格式数据加载、处理与存储

多种格式数据加载.处理与存储 实际的场景中,我们会在不同的地方遇到各种不同的数据格式(比如大家熟悉的csv与txt,比如网页HTML格式,比如XML格式),我们来一起看看python如何和这些格式的数 ...

mismatch位置(MD tag)- sam&sol;bam格式解读进阶

这算是第二讲了,前面一讲是:Edit Distance编辑距离(NM tag)- sam/bam格式解读进阶 MD是mismatch位置的字符串的表示形式,貌似在call SNP和indel的时候会用 ...

SAMTOOLS使用 SAM BAM文件处理

[怪毛匠子 整理] samtools学习及使用范例,以及官方文档详解 #第一步:把sam文件转换成bam文件,我们得到map.bam文件 system"samtools view -bS m ...

SAM&sol;BAM文件处理

当测序得到的fastq文件map到基因组之后,我们通常会得到一个sam或者bam为扩展名的文件.SAM的全称是sequence alignment/map format.而BAM就是SAM的二进制文件 ...

文件格式——Sam&amp&semi;bam文件

Sam&bam文件 SAM是一种序列比对格式标准, 由sanger制定,是以TAB为分割符的文本格式.主要应用于测序序列mapping到基因组上的结果表示,当然也可以表示任意的多重比对结果.当 ...

随机推荐

linux GO语言配置安装

1.下载地址 https://golang.org/dl/ 2.解压 解压到/usr/local/go目录下 cd /usr/local/go bin/go version 执行如上命令出现go的版本 ...

Java面试题之Class&period;forName的作用

按参数中指定的字符串形式的类名去搜索并加载相应的类,如果该类字节码已经被加载过,则返回代表该字节码的Class实例对象,否则,按类加载器的委托机制去搜索和加载该类,如果所有的类加载器都无法加载到该类, ...

Mybatis 元素内容必须由格式正确的字符数据或标记组成

一个web应用,框架为SpringMVC Spring Mybatis ,昨天写了一下午的代码,因为逻辑较大,期间也没测,打算写完这个功能点在进行测试,谁知道写完的时候,tomcat根本启动不起来了, ...

Flask信号源码流程

1. appcontext_pushed = _signals.signal('appcontext-pushed'# 请求app上下文push时执行 return RequestContext(se ...

关于 Glassfish

GlassFish 是一款强健的商业兼容应用服务器,达到产品级质量,可免费用于开发.部署和重新分发.开发者可以免费获得源代码,还可以对代码进行更改 GlassFish 是用于构建 Java EE 5应 ...

代码动态设置edittext输入类型为密码类型

开发中常常会用到EditText输入框,要将他的输入类型设置为密码输入,但是直接在布局文件中设置时,hint字体风格就会不一样 例如,在布局文件中直接设置是这样的(如下图),字体风格明显跟上一行的不一 ...

android 蓝牙通讯编程 备忘

1.启动App后: 判断->蓝牙是否打开(所有功能必须在打牙打开的情况下才能用) 已打开: 启动代码中的蓝牙通讯Service 未打开: 发布 打开蓝牙意图(系统),根据Activity返回进场 ...

counting the buildings - 第一类斯特灵数

2017-08-10 21:10:08 writer:pprp //TLE #include #include #include &lt ...

BZOJ 1263 整数划分&lpar;数学&plus;高精度&rpar;

我们不妨考虑可以划分为实数的情况,设划分为x份实数,使得总乘积最大. 易得当每一份都相等时乘积最大.即 ans=(n/x)^x. 现在只需要求出这个函数取得最大值的时候x的取值了. 两边取对数,则有l ...

Python&lpar;进程池与协程&rpar;

1.进程池与线程池: 为什么要用“池”:池子使用来限制并发的任务数目,限制我们的计算机在一个自己可承受的范围内去并发地执行任务 池子内什么时候装进程:并发的任务属于计算密集型 池子内什么时候装线程:并 ...