使用R语言查找范围重叠

Question

我有两个数据框，每个数据框有三列：chrom、start和stop，我们称它们为rangesA和rangesB。对于rangesA的每一行，我需要找到哪些（如果有）rangesB的行完全包含了rangesA的行-我指的是rangesAChrom == rangesBChrom, rangesAStart >= rangesBStart and rangesAStop <= rangesBStop。

目前我正在做以下操作，但我不太喜欢这种方法。请注意，我正在循环处理rangesA的行，但除了阅读此特定解决方案更加清晰外，其他原因都不太重要。

rangesA：

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

rangesB：

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

对于rangesA中的每一行：

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

我猜肯定有一种更好的方法（通过更快的方式处理大范围的rangesA和rangesB），而不是循环这个结构。有什么想法吗？

- geoffjentry

Answer 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

- Aaron Statham

---

IRanges上的好指针，忘记了。由于各种原因，它并不适合我的情况，但仍然很好知道。我的玩具示例漏掉了一些关键部分（我的错），这使得我难以弄清楚如何使用IRanges，而提供的merge()解决方案则大大加快了速度。此外，虽然它可以大规模扩展，但我们也看到了许多非常小的情况，我认为S4的开销在这些情况下减慢了速度。 - geoffjentry

有没有办法计算部分重叠的情况？ - Cina

Answer 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

- Joshua Ulrich

Answer 3

data.table 包自 v1.9.4 版本起新增了一个函数 foverlaps()，它能够在区间范围内进行合并：

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()按提供的列（在本例中是所有列，因为我们没有提供任何列）对数据表进行排序，并标记这些列为已排序，稍后我们将使用它们执行连接操作。

• foverlaps()有效地执行重叠连接。有关详细解释和与其他方法的比较，请参见此答案。

- Arun

Answer 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

- Joe

---

想象一种情况，其中rangesA有20k行，rangesB有300k行 -> 极其巨大的连接 -> 不可能适合R RAM内存。使用IRanges的解决方案更好。 - Marcin

Answer 5

RangesA和RangesB显然是BED语法，可以在命令行中使用BEDtools执行此操作，它非常快速灵活，并有其他十几个选项与基因组区间一起使用。然后再将结果放回R中。 https://code.google.com/p/bedtools/

- mikyatope

---

从经验来看，如果数据集的大小适合在R中运行，那么在一个编程环境中保持一致性比频繁切换工具更容易复现结果。 - lmrta

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可以实现此操作，使得rangesA中的每一列成为一行，而rangesB中的每一行则对应一个范围:

> 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

- eyjo