在 R 中识别连续重叠的片段

Question

在 R 中识别连续重叠的片段

rdata.tablegroupingoverlaplocf

8

我需要将重叠的片段聚合成一个范围涵盖所有连接的片段。

请注意，简单的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的答案比David的评论略好。

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按start排序），但后者更快。
即时生成grp的速度要快得多，不幸的是我需要分配grp。
cumsum(cummax(shift(end, fill = 0)) < start)比cumsum(c(0, start[-1L] > cummax(head(end, -1L))))快得多。
我没有测试过GenomicRanges包的解决方案。

- rluech

4

x[, .(start[1L], end[.N]), by = .(grp = cumsum(c(0, start[-1L] > cummax(head(end, -1L)))))] 可以工作。基本上是一个数据表版本的我的解决方案在这里。 - David Arenburg

2个回答

5

您可以尝试使用 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碱基。因此，必须从片段末端减去1以获得正确的图表。

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)

- Roman

1

这个群组正是他正在寻找的。 - Onyambu

在你的编辑中你仍然在使用 x$grp=..... 这是手动输入组的方法。这不应该是这种情况。 - Onyambu

@Jimbou，你的第一条语句明确说明了“#add the grouping”，但我们没有这个分组...我们应该创建/获取该分组。请查看Uwe提供的答案。 - Onyambu

@Onyambu 我删除了我的第一个答案。现在你只需要传递 start 和 end，然后你就会得到预期的结果。对此满意吗？ - Roman

这解决了这个问题。喜欢它。 - Onyambu

@Jimbou 太好了！感谢您更新图表。由于它们不再需要，我将删除我的评论。 - Uwe

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Uwe · Accepted Answer

这位楼主请求将重叠的片段合并为一个包含所有连接片段的单个片段。

以下是另一种解决方案，它使用 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

声明：我在SO上看过一个聪明的方法，但我不记得确切的位置。

编辑：

正如David Arenburg指出的那样，在by =参数中可以即时完成创建grp变量的工作，无需单独创建。

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)]