R语言 具有起始位置和结束位置的重叠联接

o0lyfsai  于 2023-02-01  发布在  其他
关注(0)|答案(5)|浏览(128)

考虑下面的data.table s。第一个定义了一组区域,每个组"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

最后,我想提取"d2"中的行,其中"pos"在每个组x中的"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"已经丢失了。有没有办法保留所有带滚动连接的列,这样我就可以按需要报告"start","pos","end"了?

8ljdwjyq

8ljdwjyq1#

Overlap joins是在data.table v1.9.3中使用commit 1375实现的,并且在current stable release, v1.9.4中可用。函数名为foverlaps。从NEWS中:
29)Overlap joins#528现在终于来了!!除了type="equal"maxgapminoverlap参数之外,其他的都实现了。查看?foverlaps及其用法示例。这是data.table的一个主要特性。
让我们考虑x,一个定义为[a, b]的区间,其中a <= b,和y,另一个定义为[c, d]的区间,其中c <= d。区间y被认为 * 完全重叠 * x,当且仅当d >= ac <= b1。并且y完全包含 * 在 * x内,iff a <= c,d <= b2.对于实现的不同类型的重叠,请查看?foverlaps
您的问题是重叠联接的一个特例:在d1中,你有startend位置的真实物理间隔。另一方面,在d2中,只有位置(pos),而不是interval。为了能够进行重叠连接,我们还需要在d2中创建interval。这可以通过创建一个额外的变量pos2来实现,这与posd2[, pos2 := pos])相同。因此,我们现在在d2中有一个区间,尽管具有相同的 startend 坐标。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(x),如果d2不存在,则使用key(y)。但是d2不存在键,我们无法设置y中的列,因为它们的名称不相同。因此,我们显式设置by.x
重叠的 * 类型 * 是 within,我们希望有 * 所有 * 匹配,只有在有匹配的情况下。
注意:foverlaps在底层使用了data.table的二分查找特性(必要时沿着使用了roll),但一些函数参数(重叠类型、maxgap、minoverlap等)受到Bioconductor包IRanges中的findOverlaps()函数的启发,这是一个优秀的包(GenomicRanges也是如此,它为Genomics扩展了IRanges)。
那么优势是什么呢?
对上述代码进行基准测试,您的数据结果比Gabor的答案慢x1m42 n1倍(计时:Gabor的数据。表格解决方案= 0.004 vs foverlaps = 0.021秒)。但在这个粒度下,这真的很重要吗?
真正有趣的是看看它在 * 速度 * 和 * 内存 * 方面的伸缩性如何。在Gabor的答案中,我们基于键列x进行连接。* 然后 * 过滤结果。
如果d1有大约4万行,d2有10万行(或更多),那么对于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是420 MB。这里花费的大部分时间实际上是在子集上。你可以通过设置参数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。所以,它肯定比上面显示的使用更少的内存。
[The从ans1中删除pos2并在两个答案上设置密钥后,验证答案相同。]
请注意,此重叠连接在设计时会遇到一些问题,即d2不一定具有相同的起始坐标和结束坐标(例如:基因组学,我来自的领域,其中d2通常是大约30-150百万或更多行)。
foverlaps()是稳定的,但仍在开发中,这意味着一些参数和名称可能会更改。
NB:既然我在上面提到了GenomicRanges,它也完全有能力解决这个问题。它在引擎盖下使用interval trees,而且内存效率也相当高。在我对基因组数据的基准测试中,foverlaps()更快。但这是另一篇(博客)文章,改天吧。

vyswwuz2

vyswwuz22#

data.table v1.9.8+有一个新特性-- non-equi joins。有了这个特性,这个操作就变得更加简单了:

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
oogrdqng

oogrdqng3#

    • 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解。我们还使用数据表中的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
axr492tv

axr492tv4#

重叠连接在dplyr 1.1.0中通过函数join_by可用。
对于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
7kjnsjlb

7kjnsjlb5#

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

相关问题