发现数百万范围/间隔之间的重叠

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 -- 重叠长度

我在具有高达 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. 删除区间时，请查看集合中的其他区间。它们都与被移除的区间重叠，并且按照重叠的长度排序，最长的在前。遍历集合直到长度太短，并报告每个重叠。

Answer 2

pyranges作者在这里。感谢您试用我的图书馆。

您的数据有多大？当两个 PyRanges 都是 1e7 行时，在 200 GB 内存的繁忙服务器上使用 24 个内核，由 pyranges 完成的繁重工作大约需要 12 秒。

设置：

import pyranges as pr
import numpy as np
import mkl

mkl.set_num_threads(1)

### Create data ###
length = int(1e7)
minimum_overlap = 5

gr = pr.random(length)
gr2 = pr.random(length)

# add ids
gr.id1 = np.arange(len(gr))
gr2.id2 = np.arange(len(gr))

# add lengths
gr.length1 = gr.lengths()
gr2.length2 = gr2.lengths()
gr

# +--------------+-----------+-----------+--------------+-----------+-----------+
# | Chromosome   | Start     | End       | Strand       | id1       | length1   |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   | (int32)   |
# |--------------+-----------+-----------+--------------+-----------+-----------|
# | chr1         | 146230338 | 146230438 | +            | 0         | 100       |
# | chr1         | 199561432 | 199561532 | +            | 1         | 100       |
# | chr1         | 189095813 | 189095913 | +            | 2         | 100       |
# | chr1         | 27608425  | 27608525  | +            | 3         | 100       |
# | ...          | ...       | ...       | ...          | ...       | ...       |
# | chrY         | 21533766  | 21533866  | -            | 9999996   | 100       |
# | chrY         | 30105890  | 30105990  | -            | 9999997   | 100       |
# | chrY         | 49764407  | 49764507  | -            | 9999998   | 100       |
# | chrY         | 3319478   | 3319578   | -            | 9999999   | 100       |
# +--------------+-----------+-----------+--------------+-----------+-----------+
# Stranded PyRanges object has 10,000,000 rows and 6 columns from 25 chromosomes.

进行分析：

j = gr.join(gr2, new_pos="intersection", nb_cpu=24)
# CPU times: user 3.85 s, sys: 3.56 s, total: 7.41 s
# Wall time: 12.3 s
j.overlap = j.lengths()
out = j.df["id1 id2 length1 length2 overlap".split()]

out = out[out.overlap >= minimum_overlap]
out.head()

       id1     id2  length1  length2  overlap
1    2  485629      100      100       74
2    2  418820      100      100       92
3    3  487066      100      100       13
4    7  191109      100      100       31
5   11  403447      100      100       76