pybedtools:对bedtools的封装和扩展

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

1. BedTool

2. feature

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

1. BedTool

构建BedTool对象的代码如下

>>> import pybedtools
>>> a = pybedtools.BedTool('a.bed')
bedtools这个软件相关的子命令,都在pybedtools中进行了实现,以取交集为例,对应的用法如下
>>> 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子命令的相关选项,则以函数参数的形式封装在对应的函数中。

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

>>> a_and_b.saveas('intersect.bed', trackline='track name')
>>> a_and_b.moveto('intersect.bed')

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

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

>>> a.intersect(b, u=True).saveas('a-with-b.bed').merge().saveas('a-with-b-merged.bed')

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

>>> x1 = a.intersect(b, u=True).merge()
>>> x2 = (a + b).merge()
>>> x1 == x2
True

2. feature

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

>>> a[0]
Interval(chr1:1-100)
>>> feature = a[0]
>>> print(feature)
chr1 1 100 feature1 0 +

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

>>> feature[0]
'chr1'
>>> feature[1]
'1'

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

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

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

>>> feature.fields
[u'chr1', u'1', u'100', u'feature1', u'0', u'+']

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

>>> a.filter(lambda x:len(x)> 100)

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

>>> 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·
(0)

相关推荐

  • 8 比对及找变异步骤的质控

    使用qualimap对wes的比对bam人家总结测序深度和覆盖度ls -lh *raw.vcf-rwxrwxrwx 1 root root 184M Jun 7 10:58 SRR7696207_ra ...

  • 使用bedtools根据染色体上的起止位置拿到基因symbol

    大家在进行各种组学研究的过程中都可能会遇到想要看看包含了某一段染色体的基因有哪些的情况.如果只有几个位置坐标,那可以很方便的在USUC.NCBI或IGV中进行查看,但如果有很多位置坐标,需要批量得到基 ...

  • SNV突变(96种)频谱的制作

    昨天我们学习了正常情况下,6种SNV(C>A, C>G, C>T, T>A, T>C, T>G)突变频谱的制作,但是如果考虑到突变的上下文,就可以变成96种(如下图 ...

  • 根据比对后的bam文件来计算外显子的RPKM

    bam2rpkm by bedtools 比对好的bam文件一般需要根据gtf文件来根据 genomic features 进行计数,但是htseq-counts或者featureCounts这样的软 ...

  • Harvard FAS Informatics出品的ATAC

    Harvard FAS Informatics出品的ATAC-seq测序指南 github链接:harvardinformatics/ATAC-seq 参考文献:ATAC-seq: A Method ...

  • 生信编程18.把文件内容按照染色体分开写出

    有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物.最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的. ...

  • R语言版本的bedtools--bedtoolsr

    bedtools是一款非常强大的用于处理bed,vcf,gff等格式数据的工具,该软件由犹他大学的Quinlan实验室开发.但是目前bedtools主要提供的是在linux,unxi等操作系统环境下的 ...

  • 对参考基因组按照200k分区间统计测序深度及GC含量

    以前是自己写脚本: [直播]我的基因组47:测序深度和GC含量的关系 可能是太复杂,大多数读者表示看不懂,所以我重新使用已有的轮子来做这件事. 下载hg38参考基因组 直接谷歌搜索即可: 拿到下载链接 ...

  • 生信编程​9.bedtools从bam文件中提取目的序列

    有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导,而且他们都有了自己的公众号,知乎号,也算是一番人物.最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师,挺有意思的. ...

  • lncRNA组装流程的软件介绍之bedtools

    咱们<生信技能树>的B站有一个lncRNA数据分析实战,缺乏配套笔记,所以我们安排了100个lncRNA组装案例文献分享,以及这个流程会用到的100个软件的实战笔记教程! 下面是100个l ...