13

我有两个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 的大型实例上更快)的方法来执行此操作,而不是循环遍历此构造。有任何想法吗?

4

6 回答 6

21

使用来自 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
于 2010-10-13T00:54:08.020 回答
14

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

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
于 2010-10-12T15:45:41.743 回答
12

自 v1.9.4 以来,该data.table包具有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()按提供的列(在本例中为所有列,因为我们没有提供任何列)对 data.table 进行排序,并将这些列标记为已排序,稍后我们将使用它来执行连接。

  • foverlaps()有效地进行重叠连接。有关详细说明以及与其他方法的比较,请参阅此答案。

于 2014-12-27T22:49:52.553 回答
3

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

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

于 2013-07-17T15:05:37.593 回答
3

我添加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
于 2016-12-01T16:43:17.997 回答
2

对于您的示例数据:

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,使得每一列是 rangeA 中的一行,每一行是 rangeB 中的对应行:

> 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
于 2010-10-12T17:53:25.983 回答