识别R中连续重叠的段

Question

我需要将重叠的段聚合成一个范围内的所有连接段.

请注意,简单的foverlaps无法检测非重叠但已连接的段之间的连接,请参阅示例以进行说明.如果在我的情节中我的部分下雨,我正在寻找干燥的地面.

到目前为止,我通过迭代算法解决了这个问题,但我想知道这个问题是否有更优雅,更直接的前进方式.我肯定不是第一个面对它的人.

我正在考虑非平等滚动连接,但是要实现它

library(data.table)
(x <- data.table(start = c(41,43,43,47,47,48,51,52,54,55,57,59), 
                  end = c(42,44,45,53,48,50,52,55,57,56,58,60)))

#     start end
#  1:    41  42
#  2:    43  44
#  3:    43  45
#  4:    47  53
#  5:    47  48
#  6:    48  50
#  7:    51  52
#  8:    52  55
#  9:    54  57
# 10:    55  56
# 11:    57  58
# 12:    59  60

setorder(x, start)[, i := .I] # i is just a helper for plotting segments
plot(NA, xlim = range(x[,.(start,end)]), ylim = rev(range(x$i)))
do.call(segments, list(x$start, x$i, x$end, x$i))

x$grp <- c(1,3,3,2,2,2,2,2,2,2,2,4) # the grouping I am looking for
do.call(segments, list(x$start, x$i, x$end, x$i, col = x$grp))
(y <- x[, .(start = min(start), end = max(end)), k=grp])

#    grp start end
# 1:   1    41  42
# 2:   2    47  58
# 3:   3    43  45
# 4:   4    59  60

do.call(segments, list(y$start, 12.2, y$end, 12.2, col = 1:4, lwd = 3))

编辑:

这很棒,谢谢,cummax和cumsum完成这项工作,Uwe的答案稍微好于戴维斯评论.

• end[.N]可以得到错误的结果,尝试x下面的示例数据. max(end)在所有情况下都是正确的,并且更快.

  x <- data.table(start = c(11866, 12696, 13813, 14011, 14041), end = c(13140, 14045, 14051, 14039, 14045))

• min(start)并start[1L]给出相同的(x按开始排序),后者更快.
• 飞行中的grp明显更快,不幸的是我需要分配grp.
• cumsum(cummax(shift(end, fill = 0)) < start)明显快于cumsum(c(0, start[-1L] > cummax(head(end, -1L)))).
• 我没有测试GenomicRanges包的解决方案.

Answer 1

OP已要求将重叠段聚合成一个范围内的所有连接段.

这是另一种使用cummax()和cumsum()识别重叠或相邻片段组的解决方案:

x[order(start, end), grp := cumsum(cummax(shift(end, fill = 0)) < start)][
  , .(start = min(start), end = max(end)), by = grp]

   grp start end
1:   1    41  42
2:   2    43  45
3:   3    47  58
4:   4    59  60

免责声明:我已经看到了其他地方的聪明方法,但我不记得到底在哪里.

编辑:

正如David Arenburg指出的那样,没有必要grp单独创建变量.这可以在参数中即时完成by =:

x[order(start, end), .(start = min(start), end = max(end)), 
  by = .(grp = cumsum(cummax(shift(end, fill = 0)) < start))]

可视化

可以修改OP的图以显示聚合的段(快速和脏):

plot(NA, xlim = range(x[,.(start,end)]), ylim = rev(range(x$i)))
do.call(segments, list(x$start, x$i, x$end, x$i))
x[order(start, end), .(start = min(start), end = max(end)), 
  by = .(grp = cumsum(cummax(shift(end, fill = 0)) < start))][
    , segments(start, grp + 0.5, end, grp + 0.5, "red", , 4)]

Answer 2

你可以尝试一种GenomicRanges方法.在输出中,每行都是一个组.

library(GenomicRanges)
x_gr <- with(x, GRanges(1, IRanges(start, end)))
as.data.table(reduce(x_gr, min.gapwidth=0))[,2:3]
   start end
1:    41  42
2:    43  45
3:    47  58
4:    59  60

并且可以使用视觉检查Gviz.在这里,我们必须知道该包是为生物学家和遗传信息而建的.背后的模式是DNA碱基.因此,必须减去段末端中的一个以获得正确的绘图.

library(Gviz)
ga <- Gviz::GenomeAxisTrack()
xgr <- with(x, GRanges(1, IRanges(start, end = end - 1)))
xgr_red <- reduce(xgr, min.gapwidth=1)
ga <- GenomeAxisTrack()
GT <- lapply(xgr, GeneRegionTrack)
GT_red <- lapply(xgr_red, GeneRegionTrack, fill = "lightblue")
plotTracks(c(ga, GT, GT_red),from = min(x$start), to = max(x$start)+2)