具有开始和结束位置的重叠连接

Posted 2023-02-23

【中文标题】具有开始和结束位置的重叠连接

【英文标题】：Overlap join with start and end positions

【发布时间】：2014-08-20 05:29:26

【问题描述】：

考虑以下data.tables。第一个定义了一组区域，其中每个组“x”的开始和结束位置：

library(data.table)

d1 <- data.table(x = letters[1:5], start = c(1,5,19,30, 7), end = c(3,11,22,39,25))
setkey(d1, x, start)

#    x start end
# 1: a     1   3
# 2: b     5  11
# 3: c    19  22
# 4: d    30  39
# 5: e     7  25

第二个数据集具有相同的分组变量“x”，并在每个组内定位“pos”：

d2 <- data.table(x = letters[c(1,1,2,2,3:5)], pos = c(2,3,3,12,20,52,10))
setkey(d2, x, pos)

#    x pos
# 1: a   2
# 2: a   3
# 3: b   3
# 4: b  12
# 5: c  20
# 6: d  52
# 7: e  10

最终，我想在每个组x 中提取“d2”中“pos”在“start”和“end”定义的范围内的行。想要的结果是

#    x pos start  end
# 1: a   2     1    3
# 2: a   3     1    3
# 3: c  20    19   22
# 4: e  10     7   25

任何组x 的开始/结束位置永远不会重叠，但可能存在不在任何区域的值间隙。

现在，我相信我应该使用滚动联接。据我所知，我不能在连接中使用“结束”列。

我试过了

d1[d2, roll = TRUE, nomatch = 0, mult = "all"][start <= end]

得到了

#    x start end
# 1: a     2   3
# 2: a     3   3
# 3: c    20  22
# 4: e    10  25

这是我想要的正确行集；然而“pos”变成了“start”，原来的“start”已经丢失。有没有办法通过滚动连接保留所有列，以便我可以根据需要报告“开始”、“位置”、“结束”？

【参考方案1】：

Overlap joins 在data.table v1.9.3 中使用commit 1375 实现，并且在current stable release, v1.9.4 中可用。该函数称为foverlaps。来自NEWS：

29) Overlap joins #528 现在终于来了！！除了 type="equal" 和 maxgap 和 minoverlap 参数外，其他所有内容都已实现。查看?foverlaps 及其用法示例。这是data.table 的一个主要功能。

让我们考虑定义为[a, b] 的区间x，其中a <= b，以及定义为[c, d] 的另一个区间y，其中c <= d。如果 d >= a 和 c <= b 1，则区间 y 完全重叠 x。并且 y 完全包含在内 x，如果 a <= c,d <= b 2。对于实现的不同类型的重叠，请查看?foverlaps。

您的问题是重叠连接的一种特殊情况：在d1 中，您有start 和end 位置的真实物理间隔。另一方面，在d2 中，只有位置（pos），没有间隔。为了能够进行重叠连接，我们还需要在d2 中创建间隔。这是通过创建一个额外的变量pos2 来实现的，它与pos (d2[, pos2 := pos]) 相同。因此，我们现在在d2 中有一个区间，尽管具有相同的 start 和 end 坐标。然后可以在foverlap 中使用d2 中的这个“虚拟、零宽度间隔”与d1 进行重叠连接：

require(data.table) ## 1.9.3
setkey(d1)
d2[, pos2 := pos]
foverlaps(d2, d1, by.x = names(d2), type = "within", mult = "all", nomatch = 0L)
#    x start end pos pos2
# 1: a     1   3   2    2
# 2: a     1   3   3    3
# 3: c    19  22  20   20
# 4: e     7  25  10   10

by.y 默认是key(y)，所以我们跳过了它。 by.x 如果存在则默认采用key(x)，如果不存在则采用key(y)。但是d2 的键不存在，并且我们不能设置来自y 的列，因为它们的名称不同。所以，我们明确地设置了by.x。

重叠的类型是within，我们希望所有匹配，前提是存在匹配。 p>

注意：foverlaps 在底层使用了 data.table 的二分搜索功能（必要时还使用了roll），但一些函数参数（重叠类型、maxgap、minoverlap 等）受到函数 @ 的启发987654371@ 来自 Bioconductor 包IRanges，一个优秀的包（GenomicRanges 也是如此，它扩展了IRanges 用于基因组学）。

---

那么有什么优势呢？

上述代码对您的数据的基准测试导致foverlaps() 比 Gabor 的答案慢（计时：Gabor 的 data.table 解决方案 = 0.004 与 foverlaps = 0.021 秒）。但在这种粒度上真的重要吗？

真正有趣的是看看它的扩展性如何 - 就速度和内存而言。在 Gabor 的回答中，我们基于键列 x 加入。 然后过滤结果。

如果d1 有大约 40K 行，而d2 有 100K 行（或更多）怎么办？对于d2 中与d1 中的x 匹配的每一行，所有 这些行将被匹配并返回，仅在以后过滤。这是您的 Q 仅略微缩放的示例：

生成数据：

require(data.table)
set.seed(1L)
n = 20e3L; k = 100e3L
idx1 = sample(100, n, TRUE)
idx2 = sample(100, n, TRUE)
d1 = data.table(x = sample(letters[1:5], n, TRUE), 
                start = pmin(idx1, idx2), 
                end = pmax(idx1, idx2))

d2 = data.table(x = sample(letters[1:15], k, TRUE), 
                pos1 = sample(60:150, k, TRUE))

重叠：

system.time(
    setkey(d1)
    d2[, pos2 := pos1]
    ans1 = foverlaps(d2, d1, by.x=1:3, type="within", nomatch=0L)
)
# user  system elapsed 
#   3.028   0.635   3.745 

这总共占用了大约 1GB 的内存，其中ans1 是 420MB。大部分时间都花在子集上。您可以通过设置参数verbose=TRUE来检查它。

Gabor 的解决方案：

## new session - data.table solution
system.time(
    setkey(d1, x)
    ans2 <- d1[d2, allow.cartesian=TRUE, nomatch=0L][between(pos1, start, end)]
)
#   user  system elapsed 
# 15.714   4.424  20.324 

这总共花费了大约 3.5GB。

我刚刚注意到 Gabor 已经提到了中间结果所需的内存。所以，试试sqldf：

# new session - sqldf solution
system.time(ans3 <- sqldf("select * from d1 join 
            d2 using (x) where pos1 between start and end"))
#   user  system elapsed 
# 73.955   1.605  77.049 

总共占用了 ~1.4GB。因此，它使用的内存肯定比上面显示的要少。

[在从ans1 中删除pos2 并在两个答案上设置密钥后，验证答案相同。]

请注意，此重叠连接的设计问题是 d2 不一定具有相同的开始和结束坐标（例如：基因组学，我来自的领域，d2 通常约为 30-1.5 亿或更多行）。

---

foverlaps() 稳定，但仍在开发中，这意味着某些参数和名称可能会更改。

NB：既然我上面提到了GenomicRanges，它也完全能够解决这个问题。它在后台使用interval trees，并且内存效率也很高。在我的基因组数据基准测试中，foverlaps() 更快。但那是另一篇（博客）帖子，其他时间。

【讨论】：

你写过那篇博文吗？我试图用一组范围重叠人类基因组中的几乎每个核苷酸（为相关指标评分），我在大约四个小时后杀死了 Granges 的 findOverlaps 以从 foverlaps 切线开始。今天，有了这篇带有说明的现成博客文章，可以为我节省很多时间。 :D

【参考方案2】：

data.table v1.9.8+ 有一个新功能 - non-equi 加入。这样，这个操作就变得更加简单了：

require(data.table) #v1.9.8+
# no need to set keys on `d1` or `d2`
d2[d1, .(x, pos=x.pos, start, end), on=.(x, pos>=start, pos<=end), nomatch=0L]
#    x pos start end
# 1: a   2     1   3
# 2: a   3     1   3
# 3: c  20    19  22
# 4: e  10     7  25

【讨论】：

您是否介意为此方法添加内存使用和计时

这里需要注意的是，x 指的是两个东西。第一个 x 是“x”列，x.pos 中的第二个 x 指的是 d2。

【参考方案3】：

1) sqldf 这不是 data.table，但复杂的连接条件很容易在 SQL 中以直接的方式指定：

library(sqldf)

sqldf("select * from d1 join d2 using (x) where pos between start and end")

给予：

  x start end pos
1 a     1   3   2
2 a     1   3   3
3 c    19  22  20
4 e     7  25  10

2) data.table 对于 data.table 的答案，试试这个：

library(data.table)

setkey(d1, x)
setkey(d2, x)
d1[d2][between(pos, start, end)]

给予：

   x start end pos
1: a     1   3   2
2: a     1   3   3
3: c    19  22  20
4: e     7  25  10

请注意，这确实具有形成可能较大的中间结果 d1[d2] 的缺点，而 SQL 可能不会这样做。剩下的解决方案也可能有这个问题。

3) dplyr 这建议了相应的 dplyr 解决方案。我们还使用 data.table 中的between：

library(dplyr)
library(data.table) # between

d1 %>% 
   inner_join(d2) %>% 
   filter(between(pos, start, end))

给予：

Joining by: "x"
  x start end pos
1 a     1   3   2
2 a     1   3   3
3 c    19  22  20
4 e     7  25  10

4) 合并/子集 仅使用 R 的基数：

subset(merge(d1, d2), start <= pos & pos <= end)

给予：

   x start end pos
1: a     1   3   2
2: a     1   3   3
3: c    19  22  20
4: e     7  25  10

已添加请注意，此处的数据表解决方案比其他答案中的解决方案要快得多：

dt1 <- function() 
 d1 <- data.table(x=letters[1:5], start=c(1,5,19,30, 7), end=c(3,11,22,39,25))
 d2 <- data.table(x=letters[c(1,1,2,2,3:5)], pos=c(2,3,3,12,20,52,10))
 setkey(d1, x, start)
 idx1 = d1[d2, which=TRUE, roll=Inf] # last observation carried forwards

 setkey(d1, x, end)
 idx2 = d1[d2, which=TRUE, roll=-Inf] # next observation carried backwards

 idx = which(!is.na(idx1) & !is.na(idx2))
 ans1 <<- cbind(d1[idx1[idx]], d2[idx, list(pos)])


dt2 <- function() 
 d1 <- data.table(x=letters[1:5], start=c(1,5,19,30, 7), end=c(3,11,22,39,25))
 d2 <- data.table(x=letters[c(1,1,2,2,3:5)], pos=c(2,3,3,12,20,52,10))
 setkey(d1, x)
 ans2 <<- d1[d2][between(pos, start, end)]


all.equal(as.data.frame(ans1), as.data.frame(ans2))
## TRUE

benchmark(dt1(), dt2())[1:4]
##     test replications elapsed relative
##  1 dt1()          100    1.45    1.667  
##  2 dt2()          100    0.87    1.000  <-- from (2) above

【讨论】：

感谢您的全面性。我只是真的想要一些东西来利用 data.table 滚动连接中实现的优化树搜索。感谢sqldf 的建议。我习惯于在 SQL Server 中编写像 sqldf("select * from d1 join d2 on d1.x==d2.x and d2.pos>=d1.start and d2.pos<=d1.end") 这样的连接，当您可以向表中添加索引时，这很好。我认为这可以避免在内存中创建完整的外部连接（但没有测试）

但是，注意这里的data.table代码更短更快。

你说得对，between 更快。我还用你的安全带测试了d1[d2, roll=T, nomatch=0, mult="all"][start<=end]，并且接近于两者之间。它只是表明您必须测试所有内容，因为您永远不知道什么会更快。感谢您花时间检查这些。很有意思。也许当@Arun 在 data.table 中实现“范围连接”时，他可以让它更快。

您能解释一下“d1[d2][between(pos, start, end)]”的工作原理吗？

d1[d2] 连接 d1 和 d2 以及 x 和 [between(...)] 选择那些 between(...) 为 TRUE 的行。

【参考方案4】：

使用fuzzyjoin：

result <- fuzzyjoin::fuzzy_inner_join(d1, d2, 
                           by = c('x', 'pos' = 'start', 'pos' = 'end'),
                           match_fun = list(`==`, `>=`, `<=`))
result

#  x.x     pos x.y   start   end
#  <chr> <dbl> <chr> <dbl> <dbl>
#1 a         2 a         1     3
#2 a         3 a         1     3
#3 c        20 c        19    22
#4 e        10 e         7    25

由于fuzzyjoin 返回所有列，我们可能需要进行一些清理以保留我们想要的列。

library(dplyr)
result %>% select(x = x.x, pos, start, end)

# A tibble: 4 x 4
#  x       pos start   end
#  <chr> <dbl> <dbl> <dbl>
#1 a         2     1     3
#2 a         3     1     3
#3 c        20    19    22
#4 e        10     7    25