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

Question

考虑以下的 data.table。第一个定义了一组带有起始和结束位置的区域 '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 的起始/结束位置永远不会重叠，但可能存在没有任何区域的值间隙。

现在，我认为我应该使用滚动连接。从我所看到的情况来看，我无法在连接中使用“end”列。

我已经尝试过。

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

and got

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

哪一组行是我想要的正确集合; 然而“pos”已经变成了“start”，原始的“start”已经丢失了。是否有一种方法可以通过滚动连接来保留所有列，以便我可以报告所需的“start”，“pos”，“end”？

- MrFlick

Answer 1

Overlap joins在data.table v1.9.3中通过commit 1375实现，并且在当前稳定版本v1.9.4中可用。该函数名为foverlaps。来自NEWS：

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

让我们考虑一个区间x，定义为[a，b]，其中a <= b，以及另一个区间y，定义为[c，d]，其中c <= d。当且仅当d>= a 并且 c <= b1时，称区间y与x重叠。当且仅当a <= c，d <= b2时，称y完全包含在x中。有关实施的不同类型的重叠，请查看?foverlaps。

您的问题是重叠连接的一个特殊情况：在d1中，您具有真实的物理间隔，具有start和end位置。另一方面，在d2中只有位置（pos），没有间隔。为了能够进行重叠连接，我们需要在d2中创建间隔。这可以通过创建一个额外的变量pos2来实现，该变量与pos相同（d2[, pos2 := pos]）。因此，现在我们在d2中具有一个间隔，尽管其start和end坐标相同。然后可以使用d2中的这个“虚拟的零宽度间隔”来在foverlap中与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，我们希望在只有匹配项时获得 所有 匹配项。

注意：在底层，foverlaps使用了data.table的二分查找功能（在必要时还使用了roll），但是一些函数参数（如重叠类型、最大间隔、最小重叠等）受到了IRanges包中的findOverlaps()函数的启发，这是一个很棒的包（同样GenomicRanges也是优秀的，它为基因组学扩展了IRanges）。

---

那么有什么优势呢？

对上述代码的基准测试在您的数据上进行，结果表明foverlaps()比Gabor的答案慢（时间：Gabor的data.table解决方案= 0.004，而foverlaps = 0.021秒）。但是，在这个粒度下真的很重要吗？

真正有趣的是看它在速度和内存方面的扩展性如何。在Gabor的答案中，我们基于关键列x进行连接。然后过滤结果。

如果d1有大约40K行，而d2有100K行（或更多）呢？对于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))

foverlaps：

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 

Took a total of ~1.4GB. So, it definitely uses less memory than the one shown above.

[在删除ans1中的pos2并在两个答案上设置键后，已验证答案相同。]

请注意，此重叠连接是针对问题设计的，其中d2的起始和结束坐标不一定相同（例如：基因组学，我来自这个领域，其中d2通常约为3000万或更多行）。

---

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

NB：由于我上面提到了GenomicRanges，它也完全能够解决这个问题。它在内部使用interval trees，并且非常节省内存。在我的基因组数据基准测试中，foverlaps()更快。但这是另一个（博客）帖子，以后再说。

- Arun

---

你是否写过那篇博客文章？我正试图将人类基因组中几乎每个核苷酸与一组范围重叠（根据相关度量进行评分），但在 GRanges 的 findOverlaps 执行了四个小时后，我决定转而使用 foverlaps。如果有一篇包含指令的现成博客文章，今天就可以省下很多时间 :D - GenesRus

Answer 2

data.table v1.9.8+推出了一个新的功能——非等值连接。有了这个功能，该操作变得更加简单：

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

- Arun

---

你介意为这种方法添加内存使用和时间吗？ - Rentrop

这里需要注意的是，x指的是两个不同的东西。第一个x是列“x”，而x.pos中的第二个x指的是d2。 - Ryan Ward Valverde

@RyanWardValverde 为什么需要包括 pos=x.pos？假设在 d2 中有更多的列，例如 'a' 和 'b'。当选择这些列时，a 和 b 不需要以 x 为前缀，但是 pos 需要。此外，如果您不使用 x 作为 pos 的前缀，则 pos==start。您知道为什么会出现这种情况吗？ - Jamie

似乎是操作顺序的副产品。在j位置中删除括号之间的所有内容，即.(...)会给我们带来d2[d1, on = .(x, pos>=start, pos<=end), nomatch = NULL]，并产生名为pos和pos.1的列，它们对应于start和end，这对我来说是意外的行为。这似乎表明连接首先完成，并且要访问x位置中的data.table，需要显式调用它。 - Ryan Ward Valverde

Answer 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

- G. Grothendieck

---

我很欣赏您的全面性。 我只是想要利用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")，当您可以向表添加索引时非常好。 我认为这可能避免在内存中创建完整外部联接（但未测试）。 - MrFlick

然而，请注意，这里的data.table代码更短且更快。 - G. Grothendieck

你是对的，between 更快。我也用你的测试工具测试了 d1[d2, roll=T, nomatch=0, mult="all"][start<=end]，它与 between 很接近。这说明你必须测试所有东西，因为你永远不知道哪个更快。感谢你抽出时间来检查这些内容。非常有趣。也许当 @Arun在 data.table 中实现“范围连接”时，他可以让它更快。 - MrFlick

请问您能解释一下"d1[d2][between(pos, start, end)]"是如何工作的吗？ - skan

d1[d2] 沿着 x 连接 d1 和 d2，而 [between(...)] 则选择那些 between(...) 为 TRUE 的行。 - G. Grothendieck

Answer 4

通过函数join_by，可以在dplyr 1.1.0中使用重叠连接。

使用join_by，可以使用between或手动使用>=和<=进行重叠连接。

library(dplyr)
inner_join(d2, d1, by = join_by(x, between(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

inner_join(d2, d1, by = join_by(x, pos >= start, pos <= 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

- Maël

Answer 5

使用 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

- Ronak Shah