如何解决如何根据原始输入文件与python/pandas相交合并文件以将重叠划分为子区域?
我有一些来自提供外显子组测序套件的不同公司的 .bed 文件。
我想要一个文件来总结所有这些试剂盒的所有目标区域。 .bed 文件的基本结构由三列(chr#、Start、End)组成。
我想得到一个输出表,显示哪些基因组区域仅被这些试剂盒之一覆盖,哪些区域被多个(以及哪些)覆盖。 举例说明这一点的最佳方式是:
床档 1
chr# | 开始 | 结束 |
---|---|---|
1 | 100 | 300 |
床档 2
chr# | 开始 | 结束 |
---|---|---|
1 | 150 | 350 |
床档 3
chr# | 开始 | 结束 |
---|---|---|
1 | 80 | 200 |
从这些文件中,我创建了一个包含所有目标区域的数据框,并按 chr# 和 Start 坐标对其进行排序。这是结果数据框的样子:
我想合并和交叉文件以获得输出,该输出根据输入文件之间的重叠将区域划分为子区域。它应该看起来像这样:
chr# | 开始 | 结束 | Kit 1 | Kit 2 | Kit 3 |
---|---|---|---|---|---|
1 | 80 | 100 | 0 | 0 | 1 |
1 | 100 | 150 | 1 | 0 | 1 |
1 | 150 | 200 | 1 | 1 | 1 |
1 | 200 | 300 | 1 | 1 | 0 |
1 | 300 | 350 | 0 | 1 | 0 |
我知道 Bioconductor 的 GRanges 上可能有这样的功能,但我不熟悉该库及其功能。
任何帮助将不胜感激。
解决方法
更新
使用 bedtools 中的 multiIntersectBed
$ multiIntersectBed -i *.bed
1 80 100 1 3 0 0 1
1 100 150 2 1,3 1 0 1
1 150 200 3 1,2,3 1 1 1
1 200 300 2 1,2 1 1 0
1 300 350 1 2 0 1 0
还有一个 Python 接口:
http://daler.github.io/pybedtools/
>>> import pybedtools as bt
>>> from glob import glob
>>> print(bt.BedTool().multi_intersect(i=glob('*.bed')))
1 80 100 1 3 0 0 1
1 100 150 2 1,2 1 1 0
1 300 350 1 2 0 1 0
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。