如何将重叠基因组区域拆分为带重叠计数的非重叠片段

给定如下R数据框，前三列存储染色体位置信息，第四列标记该行区间是否与其他行重叠：

df <- data.frame(
  chrom = c("chr1", "chr1", "chr1", "chr1", "chr1", "chr1"),
  start = c(10, 20, 25, 30, 90,  100),
  end = c(20, 30, 38, 40, 120,  200),
  count = c("no_overlap", "overlap", "overlap", "overlap", "overlap", "overlap")
)

我们需要拆分重叠区域，生成新的数据框，其中每个子区间的count表示有多少原始区间覆盖该子区间，预期输出如下：

chrom   start end count 
 <chr>   <dbl> <dbl> <int>     
  chr1       10    20 1
  chr1       20    24 1 
  chr1       25    38 2
  chr1       39    40 1
  chr1       90    99 1
  chr1      100   120 2
  chr1      121   200 1   

解决方案1：使用基因组区间专用包（推荐） #

借助Bioconductor的IRanges和GenomicRanges包可以高效处理这类区间问题：

# 安装依赖包（首次使用时执行）
# install.packages("IRanges")
# install.packages("GenomicRanges")

library(IRanges)
library(GenomicRanges)

# 转换为GRanges对象（基因组区间标准格式）
gr <- makeGRangesFromDataFrame(df, keep.extra.columns = FALSE)

# 提取所有区间断点（起始位+结束位+1），去重排序
breaks <- sort(unique(c(start(gr), end(gr) + 1)))

# 生成所有不重叠的子区间
sub_gr <- GRanges(
  seqnames = rep(seqnames(gr)[1], length(breaks)-1),
  ranges = IRanges(start = breaks[-length(breaks)], end = breaks[-1] - 1)
)

# 计算每个子区间被原始区间覆盖的次数
coverage_counts <- countOverlaps(sub_gr, gr)

# 转换回数据框并调整格式
result_df <- as.data.frame(sub_gr)
result_df$count <- coverage_counts
result_df <- result_df[, c("seqnames", "start", "end", "count")]
colnames(result_df)[1] <- "chrom"

print(result_df)

解决方案2：基础R实现（无需额外包） #

如果不想安装Bioconductor包，可以用基础R代码实现：

# 提取所有区间的断点（起始位+结束位+1），去重排序
breaks <- sort(unique(c(df$start, df$end + 1)))

# 生成所有子区间
sub_intervals <- data.frame(
  chrom = rep(df$chrom[1], length(breaks)-1),
  start = breaks[-length(breaks)],
  end = breaks[-1] - 1
)

# 遍历每个子区间，统计覆盖它的原始区间数量
sub_intervals$count <- apply(sub_intervals, 1, function(row) {
  sum(df$chrom == row["chrom"] & df$start <= row["end"] & df$end >= row["start"])
})

print(sub_intervals)

两种方法运行后都会得到符合预期的结果，其中专用包的方法在处理大规模数据时效率更高。

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