用开始/结束窗口滚动连接
考虑下面的data.table
。 第一个定义了一组具有开始和结束位置的区域
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
第二个代表每个组的观察结果
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中的行在d1匹配x
值的区域中的行。 期望的结果是
# 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=T, 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”已经成为“开始”,原来的“开始”已经失去。 有没有办法保留所有与滚连接的列,所以我可以根据需要报告“开始”,“pos”,“结束”?
data.table v1.9.3中的 提交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 <= b
1 ,则间隔y据说与 x完全重叠 。 如果a <= c,d <= b
2 a <= c,d <= b
y完全包含在 x中。 对于实现的不同types的重叠,请看看?foverlaps
。
你的问题是一个重叠连接的特殊情况:在d1
你有真实的物理间隔的start
和end
位置。 另一方面,在d2
,只有位置( pos
),而不是间隔。 为了能够进行重叠连接,我们还需要在d2
创build间隔。 这是通过创build一个与pos
( d2[, pos2 := pos]
)相同的附加variablespos2
来实现的。 因此,我们现在在d2
有一个间隔,尽pipe起始和结束坐标相同。 然后可以在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
重叠的types 在内 ,只有匹配时我们才会有所有匹配。
注意: foverlaps
使用了data.table的二进制searchfunction(在需要的地方连同roll
),但是一些函数参数(重叠,maxgap,minoverlap等types)受到来自Bioconductor软件包IRanges
的函数IRanges
,一个很好的软件包( GenomicRanges
也是GenomicRanges
,它扩展了基因组学的IRanges
)。
那么优势是什么?
上面的代码的基准testing结果foverlaps()
比Gabor的回答慢(Timings:Gabor的data.table解= 0.004 vs foverlaps = 0.021秒)。 但在这个粒度上真的很重要吗?
真正有趣的是从速度和记忆的angular度来看,它的规模有多大。 在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))
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
是ans1
。 大部分时间都在这个子集上。 您可以通过设置参数verbose=TRUE
来检查它。
Gabor的解决scheme:
## 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。 所以,它肯定比上面显示的使用更less的内存。
[从ans1
删除pos2
并在两个答案中设置键后,答案被证实是相同的。]
请注意,这个重叠连接的devise存在问题,其中d2
不一定具有相同的开始和结束坐标(例如:基因组学,我来自哪里,其中d2
通常约为30-150百万行或更多行)。
foverlaps()
是稳定的,但仍在开发中,这意味着一些参数和名称可能会改变。
注意:由于我上面提到了GenomicRanges
,所以它也完全有能力解决这个问题。 它使用引擎盖下的间隔树 ,并且具有相当高的内存效率。 在基因组学数据的基准testing中, foverlaps()
更快。 但这是另一个(博客)的post,其他时间。
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
请注意,这有一个缺点,即形成SQL可能不会执行的可能较大的中间结果d1[d2]
。 其余的解决scheme也可能有这个问题。
3)dplyr这表明了相应的dplyr解决scheme。 我们也使用来自data.table:
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
补充注意这里的数据表解决scheme比另一个答案中的要快得多:
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 v1.9.8+
有一个新function – 非 data.table v1.9.8+
join。 这样,这个操作变得更直接:
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