如何在R的data.table中实现无第三方依赖的区间合并？

我们有如下data.table数据：

library(data.table)
df <- data.table(
  seq = c('A', 'A', 'A', 'B', 'B'),
  start = c(100, 150, 900, 3000, 6000),
  end = c(200, 300, 1000, 7000, 7000),
  id1 = c('X', 'X', 'X', 'Y', 'Y'), 
  id2 = c('A', 'A', 'A', 'B', 'C')
)

原始数据输出：

seq start   end    id1    id2
   <char> <num> <num> <char> <char>
1:      A   100   200      X      A
2:      A   150   300      X      A
3:      A   900  1000      X      A
4:      B  3000  7000      Y      B
5:      B  6000  7000      Y      C

需求说明 #

每行包含id1、id2两个分组列，以及start-end区间数据。需要按seq、id1、id2分组，合并每组内重叠或连续的区间，去除冗余。预期输出如下：

df <- data.table(
  seq = c('A', 'A', 'B', 'B'),
  start = c(100, 900, 3000, 6000),
  end = c(300, 1000, 7000, 7000),
  id1 = c('X', 'X', 'Y', 'Y'), 
  id2 = c('A', 'A', 'B', 'C')
)

对应的结果输出：

seq start   end    id1    id2
   <char> <num> <num> <char> <char>
1:      A   100   300      X      A
2:      A   900  1000      X      A
3:      B  3000  7000      Y      B
4:      B  6000  7000      Y      C

已尝试的方法 #

• 使用外部工具bedtools：需要读写外部文件后重新加载数据，流程繁琐
• 使用GRanges和Reduce函数：依赖Bioconductor包，不属于纯data.table方案

纯data.table实现方案 #

可以通过自定义函数结合data.table的分组操作完成，无需依赖第三方工具。以下是实现代码：

reduce_intervals <- function(dt) {
  setorder(dt, start)
  if (nrow(dt) <= 1) {return(dt)}
  merged_starts <- c()
  merged_ends <- c()
  current_start <- dt$start[1]
  current_end <- dt$end[1]
  for (i in 2:nrow(dt)) {
    if (dt$start[i] <= current_end) {
      current_end <- max(current_end, dt$end[i])
    } else {
      merged_starts <- c(merged_starts, current_start)
      merged_ends <- c(merged_ends, current_end)
      current_start <- dt$start[i]
      current_end <- dt$end[i]
    }
  }
  merged_starts <- c(merged_starts, current_start)
  merged_ends <- c(merged_ends, current_end)
  return(data.table(start = merged_starts, end = merged_ends))
}

result <- df[, reduce_intervals(.SD), by = .(seq, id1, id2)]

代码说明 #

1. reduce_intervals函数负责单组内的区间合并：
   • 先按start排序区间，确保处理顺序正确
   • 遍历区间，若当前区间与已合并的最后一个区间重叠/连续，则更新合并区间的end为两者最大值；否则将当前合并区间存入结果，开始新的合并区间
2. 通过data.table的by = .(seq, id1, id2)分组，对每组应用区间合并函数，最终得到去重后的区间数据

内容的提问来源于stack exchange，提问作者zhang