使用 R 查找范围内的重叠答案

【问题标题】：Finding overlap in ranges with R使用 R 查找范围内的重叠
【发布时间】：2010-10-12 15:18:18
【问题描述】：

我有两个 data.frame，每个包含三列：chrom、start 和 stop，我们称它们为 rangeA 和 rangeB。对于rangeA 的每一行，我正在寻找rangeB 中的哪一行（如果有）完全包含rangeA 行——我的意思是rangesAChrom == rangesBChrom, rangesAStart >= rangesBStart and rangesAStop <= rangesBStop。

现在我正在做以下事情，我只是不太喜欢。请注意，由于其他原因，我正在遍历 rangeA 的行，但这些原因都不是什么大问题，考虑到这个特定的解决方案，它最终只会使事情更具可读性。

范围A：

chrom   start   stop
 5       100     105
 1       200     250
 9       275     300

范围B：

chrom    start    stop
  1       200      265
  5       99       106
  9       275      290

对于 rangeA 中的每一行：

matches <- which((rangesB[,'chrom']  == rangesA[row,'chrom']) &&
                 (rangesB[,'start'] <= rangesA[row, 'start']) &&
                 (rangesB[,'stop'] >= rangesA[row, 'stop']))

我认为必须有一种更好的（更好，我的意思是在 rangeA 和 rangeB 的大型实例上更快）的方法来执行此操作，而不是循环遍历此构造。有什么想法吗？

【问题讨论】：

标签： r bioinformatics

【解决方案1】：

使用来自 Bioconductor 的 IRanges/GenomicRanges 包，它专门用于处理这些确切的问题（并且可以大规模扩展）

source("http://bioconductor.org/biocLite.R")
biocLite("IRanges")

不同染色体上的范围有几个合适的容器，一个是 RangesList

library(IRanges)
rangesA <- split(IRanges(rangesA$start, rangesA$stop), rangesA$chrom)
rangesB <- split(IRanges(rangesB$start, rangesB$stop), rangesB$chrom)
#which rangesB wholly contain at least one rangesA?
ov <- countOverlaps(rangesB, rangesA, type="within")>0

【讨论】：

IRanges 上的好指针，忘记了。我最终没有接受这个，因为由于各种原因它不适合我自己的情况，但仍然很高兴知道。我的玩具示例遗漏了几个关键位（我自己的错），这让我很难弄清楚使用 IRange 工作，而 merge() 解决方案提供了巨大的加速。此外，虽然它可以大规模扩展，但我们也看到了很多非常小的案例，我假设 S4 的开销在这些案例中会减慢它的速度。
是否也可以计算部分重叠？

【解决方案2】：

如果您可以先合并这两个对象，这将更容易/更快。

ranges <- merge(rangesA,rangesB,by="chrom",suffixes=c("A","B"))
ranges[with(ranges, startB <= startA & stopB >= stopA),]
#  chrom startA stopA startB stopB
#1     1    200   250    200   265
#2     5    100   105     99   106

【讨论】：

【解决方案3】：

RangesA 和 RangesB 显然是 BED 语法，这可以在 R 之外的命令行中使用 BEDtools 完成，非常快速和灵活，还有很多其他选项可以处理基因组区间。然后将结果再次放回 R 中。

https://code.google.com/p/bedtools/

【讨论】：

【解决方案4】：

我添加了dplyr 解决方案。

library(dplyr)
inner_join(rangesA, rangesB, by="chrom") %>% 
  filter(start.y < start.x | stop.y > stop.x)

输出：

  chrom start.x stop.x start.y stop.y
1     5     100    105      99    106
2     1     200    250     200    265

【讨论】：

想象这样一种情况，rangeA 有 20k 行，rangeB 有 300k 行 -> 疯狂的巨大连接 -> 无法放入 R RAM 内存。使用 IRange 的解决方案更好

【解决方案5】：

data.table 包有一个函数foverlaps()，它能够从 v1.9.4 开始合并跨区间范围：

require(data.table)
setDT(rangesA)
setDT(rangesB)

setkey(rangesB)
foverlaps(rangesA, rangesB, type="within", nomatch=0L)
#    chrom start stop i.start i.stop
# 1:     5    99  106     100    105
# 2:     1   200  265     200    250

setDT()通过引用将data.frame转换为data.table
setkey() 按提供的列（在本例中为所有列，因为我们没有提供任何列）对 data.table 进行排序，并将这些列标记为已排序，稍后我们将使用它来执行连接开。
foverlaps() 有效地进行重叠连接。有关详细说明以及与其他方法的比较，请参阅 this answer。

【讨论】：

【解决方案6】：

对于您的示例数据：

rangesA <- data.frame(
    chrom = c(5, 1, 9),
    start = c(100, 200, 275),
    stop = c(105, 250, 300)
)
rangesB <- data.frame(
    chrom = c(1, 5, 9),
    start = c(200, 99, 275),
    stop = c(265, 106, 290)
)

这将使用sapply 来实现，这样每一列是范围A 中的一行，每一行是范围B 中的对应行：

> sapply(rangesA$stop, '>=', rangesB$start) & sapply(rangesA$start, '<=', rangesB$stop)
      [,1]  [,2]  [,3]
[1,] FALSE  TRUE FALSE
[2,]  TRUE FALSE FALSE
[3,] FALSE FALSE  TRUE

【讨论】：