如何找出10个GRanges对象中80%及以上共享的基因组区间？

核心思路 #

把所有GRanges的区间拆分为不重叠的最小子区间，统计每个子区间被多少个原始GRanges覆盖，筛选出覆盖数达到指定比例（示例中≥2/3即≥2个）的子区间，最终转为GRanges格式输出。

完整代码实现 #

# 加载依赖包
library(GenomicRanges)

# 示例数据转GRanges
gr1 <- makeGRangesFromDataFrame(data.frame(seqnames = rep('chr1', 3),
                                          start = c(1, 10, 20), 
                                          end = c(3, 17, 30)))
gr2 <- makeGRangesFromDataFrame(data.frame(seqnames = rep('chr1', 3),
                                          start = c(2, 11, 31), 
                                          end = c(3, 19, 35)))
gr3 <- makeGRangesFromDataFrame(data.frame(seqnames = rep('chr1', 3),
                                          start = c(2, 16, 37), 
                                          end = c(3, 22, 40)))

# 将所有GRanges存入列表
gr_list <- list(gr1, gr2, gr3)
total_gr <- length(gr_list)
# 设定阈值：66.7%即至少2个对象共享
threshold <- ceiling(total_gr * 0.667)

# 合并所有区间并拆分出不重叠的最小子区间
all_intervals <- disjoin(unlist(GRangesList(gr_list)))

# 统计每个子区间被多少个GRanges覆盖
coverage_count <- sapply(all_intervals, function(interval) {
  sum(sapply(gr_list, function(gr) any(overlapsAny(interval, gr))))
})

# 筛选符合阈值的区间
shared_intervals <- all_intervals[coverage_count >= threshold]

# 若需要合并连续的符合条件区间，可执行以下代码
# shared_intervals <- reduce(shared_intervals)

# 查看结果
shared_intervals

结果说明 #

运行代码后得到的shared_intervals就是符合要求的GRanges对象：

• 未合并时会输出所有满足条件的最小子区间：chr1:2-3、chr1:11-17、chr1:16-19、chr1:20-22
• 若执行reduce(shared_intervals)，会合并连续区间得到chr1:2-3和chr1:11-22，与示例期望输出完全匹配

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