1

我有一些来自不同公司提供外显子组测序试剂盒的 .bed 文件。

我想要一个总结所有这些套件的所有目标区域的文件。.bed 文件具有由三列（chr#、Start、End）组成的基本结构。

我想得到一个输出表，显示哪些基因组区域仅被这些试剂盒之一覆盖，哪些区域被多个（以及哪些）覆盖。说明这一点的最好方法是通过一个例子：

床文件 1

字符# 开始 结尾
1 100 300

床文件 2

字符# 开始 结尾
1 150 350

床文件 3

字符# 开始 结尾
1 80 200

从这些文件中，我创建了一个包含所有目标区域的数据框，并按chr#和Start坐标对其进行排序。这是生成的数据框的样子：

概述数据框

我想合并和交叉文件以获得输出，该输出根据输入文件之间的重叠将区域划分为子区域。它应该看起来像这样：

字符# 开始 结尾 套件 1 套件 2 套件 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 上可能有这样的功能，但我不熟悉该库及其功能。

任何帮助，将不胜感激。

4

1 回答 1

1

更新

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,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
于 2021-03-20T18:13:41.707 回答