python - 寻找数百万范围/间隔之间的重叠

Question

我试图找到至少重叠用户设置的某个最小重叠长度的间隔对。间隔来自这个 pandas 数据框:

import pandas as pds

print(df1.head())
print(df1.tail())

   query_id  start_pos  end_pos  read_length orientation
0   1687655          1     4158         4158           F
1   2485364          1     7233         7233           R
2   1412202          1     3215         3215           R
3   1765889          1     3010         3010           R
4   2944965          1     4199         4199           R
         query_id  start_pos   end_pos  read_length orientation
3082467    112838   27863832  27865583         1752           F
3082468    138670   28431208  28431804          597           R
3082469    171683   28489928  28490409          482           F
3082470   2930053   28569533  28569860          328           F
3082471   1896622   28589281  28589554          274           R

start_pos 是间隔开始的位置，end_pos 是间隔结束的位置。 read_length 为区间长度。

数据按start_pos排序。

程序应具有以下输出格式:

query_id1 -- query_id2 -- read_length1 -- read_length2 -- overlap_length

我在具有高达 512gb RAM 和 4x Intel Xeon E7-4830 CPU(32 核)的计算节点上运行该程序。

我尝试运行自己的代码来查找重叠部分，但运行时间太长。

这是我试过的代码，

import pandas as pds

overlap_df = pds.DataFrame()

def create_overlap_table(df1, ovl_len):

...
(sort and clean the data here)
...

    def iterate_queries(row):
        global overlap_df
        index1 = df1.index[df1['query_id'] == row['query_id']]
        next_int_index = df1.index.get_loc(index1[0]) + 1

        if row['read_length'] >= ovl_len:
            if df1.index.size-1 >= next_int_index:
                end_pos_minus_ovlp = (row['end_pos'] - ovl_len) + 2

                subset_df = df1.loc[(df1['start_pos'] < end_pos_minus_ovlp)]
                subset_df = subset_df.loc[subset_df.index == subset_df.index.max()]
                subset_df = df1.iloc[next_int_index:df1.index.get_loc(subset_df.index[0])]
                subset_df = subset_df.loc[subset_df['read_length'] >= ovl_len]

                rows1 = pds.DataFrame({'read_id1': np.repeat(row['query_id'], repeats=subset_df.index.size), 'read_id2': subset_df['query_id']})
                overlap_df = overlap_df.append(rows1)

    df1.apply(iterate_queries, axis=1)

    print(overlap_df)

再次，我在计算节点上运行了这段代码，但在我最终取消作业之前它运行了几个小时。

我也尝试过使用两个包来解决这个问题——PyRanges，以及一个名为 IRanges 的 R 包，但它们的运行时间也太长。我看过关于区间树和一个名为 pybedtools 的 python 库的帖子，我正计划下一步研究它们。

非常感谢任何反馈

---

编辑:对于最小重叠长度，例如 800，前 5 行应如下所示，

query_id1 query_id2 read_length1 read_length2 overlap_length
1687655    2485364       4158         7233          4158
1687655    1412202       4158         3215          3215
1687655    1765889       4158         3010          3010             
1687655    2944965       4158         4199          4158
2485364    1412202       7233         3215          3215

因此，query_id1 和 query_id2 不能相同。此外，没有重复(即 A 和 B 之间的重叠不应在输出中出现两次)。

Answer 1

这是一个算法。

1. 准备一组按起点排序的区间。最初集合是空的。
2. 对所有起点和终点进行排序。
3. 遍历点。如果遇到起点，则将相应的区间添加到集合中。如果遇到终点，则从集合中移除相应的区间。
4. 删除一个区间时，查看集合中的其他区间。它们都与被删除的区间重叠，并且它们按重叠的长度排序，最长的在前。遍历集合直到长度太短，并报告每个重叠。