前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >pybedtools:对bedtools的封装和扩展

pybedtools:对bedtools的封装和扩展

作者头像
生信修炼手册
发布2020-12-28 11:40:46
7720
发布2020-12-28 11:40:46
举报
文章被收录于专栏:生信修炼手册生信修炼手册

bedtools是区间操作最常用的软件,pybedtools对其进行了封装,可以在python编程环境中灵活使用bedtools,而且进一步拓展出了很多有用的功能。

在pybedtools中, 核心是以下两个对象

1. BedTool

2. feature

BedTools表示一个区间文件,支持BED, GTF, GFF, BAM等多种格式以及对应的压缩文件,feature表示文件的每一行,包含了染色体,起始,终止位置等相关属性和信息。

1. BedTool

构建BedTool对象的代码如下

代码语言:javascript
复制
>>> import pybedtools
>>> a = pybedtools.BedTool('a.bed')

bedtools这个软件相关的子命令,都在pybedtools中进行了实现,以取交集为例,对应的用法如下

代码语言:javascript
复制
>>> a = pybedtools.BedTool('a.bed')
>>> b = pybedtools.BedTool('b.bed')
>>> a.head()
chr1 1 100 feature1 0 +
chr1 100 200 feature2 0 +
chr1 150 500 feature3 0 -
chr1 900 950 feature4 0 +
>>> b.head()
chr1 155 200 feature5 0 -
chr1 800 901 feature6 0 +
>>> a_and_b = a.intersect(b)
>>> a_and_b.head()
chr1 155 200 feature2 0 +
chr1 155 200 feature3 0 -
chr1 900 901 feature4 0 +

取交集对应的函数为intersect, 对于bedtools子命令的相关选项,则以函数参数的形式封装在对应的函数中。

进行运算之后,结果可以通过如下两种方式来保存

代码语言:javascript
复制
>>> a_and_b.saveas('intersect.bed', trackline='track name')
>>> a_and_b.moveto('intersect.bed')

saveas操作支持添加trackline, 在基因组浏览器中展示,而moveto方法只是简单的保存结果。本质上,pybedtools的运算结果,都以临时文件的形式存在,saveas相当于拷贝临时文件,生成新文件,moveto相当于重命名临时文件,速度更快。

在pybedtools中,大部分函数的返回结果都是BedTools对象,因此可以接受连缀写法,示例如下

代码语言:javascript
复制
>>> a.intersect(b, u=True).saveas('a-with-b.bed').merge().saveas('a-with-b-merged.bed')

与此同时,考虑到区间操作的特性,还进行了运算符重载,比如+号表示并集,-号表示差集,示例如下

代码语言:javascript
复制
>>> x1 = a.intersect(b, u=True).merge()
>>> x2 = (a + b).merge()
>>> x1 == x2
True

2. feature

feature对象是文件中的一行,可以通过下标来访问,示例如下

代码语言:javascript
复制
>>> a[0]
Interval(chr1:1-100)
>>> feature = a[0]
>>> print(feature)
chr1 1 100 feature1 0 +

对于feature中的每一列信息,可以通过数字下标,属性值,key这3种方式来访问,示例如下

代码语言:javascript
复制
>>> feature[0]
'chr1'
>>> feature[1]
'1'

>>> feature.start
1
>>> feature.stop
100

>>> feature['name']
'feature1
>>> feature['chrom']
chr1

feature还拥有fields属性,以列表的形式存储了对应的值,可以用于遍历操作,示例如下

代码语言:javascript
复制
>>> feature.fields
[u'chr1', u'1', u'100', u'feature1', u'0', u'+']

对于feature而言,pybedtools提供了each和filter操作,进一步扩展其功能。filter操作则用于过滤行操作,过滤掉长度小于100feature的代码,示例如下

代码语言:javascript
复制
>>> a.filter(lambda x:len(x)> 100)

each操作对每一行进行遍历,主要用于统计一些新的指标;统计每个区间长度的代码,示例如下

代码语言:javascript
复制
>>> a_and_b = a.intersect(b)
>>> print(a_and_b)
chr1 155    200    feature2 0    +
chr1 155    200    feature3 0    -
chr1 900    901    feature4 0    +

>>> def cal_len(feature):
...     feature.score = len(feature)
...     return feature
...
>>> print(a_and_b.each(cal_len))
chr1 155    200    feature2 45    +
chr1 155    200    feature3 45    -
chr1 900    901    feature4 1    +

以上只是pybedtools的基本用法,更多的用法可以查看官方的API文档。

·end·

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-12-23,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信修炼手册 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档