Поиск перекрытия в диапазонах с R

У меня есть два data.frames, каждый из которых содержит три столбца: chrom, start и stop, пусть называет их диапазонами A и диапазонамиB. Для каждой строки диапазонов A я ищу, чтобы найти, какая строка (если есть) в диапазонах B полностью содержит строку rangeA, под которой я подразумеваю rangesAChrom == rangesBChrom, rangesAStart >= rangesBStart and rangesAStop <= rangesBStop.

Сейчас я делаю следующее, что мне очень не нравится. Обратите внимание, что я зацикливаюсь на рядах диапазонов A по другим причинам, но ни одна из этих причин, скорее всего, не будет большой проблемой, а просто заканчивает тем, что делает вещи более читаемыми с учетом этого конкретного решения.

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

для каждой строки в диапазонах A:

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

Я полагаю, что лучше быть лучше (и, лучше, я имею в виду более быстрое использование больших экземпляров rangeA и диапазонов), чтобы сделать это, чем перебирать эту конструкцию. Любые идеи?

Ответы

Ответ 1

Это будет намного проще/быстрее, если вы можете сначала объединить два объекта.

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

Ответ 2

Используйте пакеты IRanges/GenomicRanges от Bioconductor, которые предназначены для решения этих точных проблем (и масштабируются в массовом порядке)

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

Ответ 3

Пакет 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() сортирует таблицу данных по предоставленным столбцам (в этом случае все столбцы, поскольку мы не предоставляли их), и маркирует эти столбцы как отсортированные, которые мы будем использовать позже для выполнения объединения на.

  • foverlaps() эффективно перекрывает перекрывающееся соединение. См. этот ответ для подробного объяснения и сравнения с другими подходами.

Ответ 4

Для данных примера:

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

Ответ 5

RangesA и RangesB являются явно синтаксисом BED, это можно сделать вне R в командной строке с помощью BEDtools, чрезвычайно быстро и гибко с помощью еще дюжины других опций для работы с геномными интервалами. Затем снова верните результаты в R.

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

Ответ 6

Я добавляю решение 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