Возможно ли использовать функцию re.table.f.foverlaps, чтобы найти пересечение перекрывающихся диапазонов в двух таблицах?

Ответ 1

@Seth обеспечил самый быстрый способ решения проблемы перекрытий пересечения с использованием функции data.table foverlaps. Однако это решение не учитывало того факта, что файлы входного ложа могут иметь перекрывающиеся диапазоны, которые необходимо было свести к отдельным областям. @Martin Morgan решил, что с его решением, использующим пакет GenomicRanges, это и пересекало, и уменьшало диапазон. Однако решение Мартина не использовало функцию foverlaps. @Arun отметил, что перекрывающиеся диапазоны в разных строках внутри таблицы в настоящее время не возможны с использованием foverlaps. Благодаря предоставленным ответам и некоторым дополнительным исследованиям stackoverflow, я придумал это гибридное решение.

Создайте пример BED файлов без перекрывающихся областей внутри каждого файла.

chr <- c(1:22,"X","Y","MT")

#bedA contains 5 million rows
bedA <- data.table(
    CHROM = as.vector(sapply(chr, function(x) rep(x,200000))),
    START = rep(as.integer(seq(1,200000000,1000)),25),
    STOP = rep(as.integer(seq(500,200000000,1000)),25),
    key = c("CHROM","START","STOP")
    )

#bedB contains 500 thousand rows
bedB <- data.table(
  CHROM = as.vector(sapply(chr, function(x) rep(x,20000))),
  START = rep(as.integer(seq(200,200000000,10000)),25),
  STOP = rep(as.integer(seq(600,200000000,10000)),25),
  key = c("CHROM","START","STOP")
)

Теперь создайте новый файл кровати, содержащий пересекающиеся области в постели A и bedB.

#This solution uses foverlaps
system.time(tmpA <- intersectBedFiles.foverlaps(bedA,bedB))

user  system elapsed 
1.25    0.02    1.37 

#This solution uses GenomicRanges
system.time(tmpB <- intersectBedFiles.GR(bedA,bedB))

user  system elapsed 
12.95    0.06   13.04 

identical(tmpA,tmpB)
[1] TRUE

Теперь измените bedA и bedB так, чтобы они содержали перекрывающиеся области:

#Create overlapping ranges
makeOverlaps <-  as.integer(c(0,0,600,0,0,0,600,0,0,0))
bedC <- bedA[, STOP := STOP + makeOverlaps, by=CHROM]
bedD <- bedB[, STOP := STOP + makeOverlaps, by=CHROM]

Время тестирования для пересечения файлов с перекрывающимися диапазонами с использованием либо foverlaps, либо функций GenomicRanges.

#This solution uses foverlaps to find the intersection and then run GenomicRanges on the result
system.time(tmpC <- intersectBedFiles.foverlaps(bedC,bedD))

user  system elapsed 
1.83    0.05    1.89 

#This solution uses GenomicRanges
system.time(tmpD <- intersectBedFiles.GR(bedC,bedD))

user  system elapsed 
12.95    0.04   12.99 

identical(tmpC,tmpD)
[1] TRUE

Победитель: foverlaps!

ИСПОЛЬЗУЕМЫЕ ФУНКЦИИ

Это функция, основанная на foverlaps, и будет вызывать функцию GenomicRanges (сокращениеBed.GenomicRanges), если существуют перекрывающиеся диапазоны (которые проверяются для использования функции rowShift).

intersectBedFiles.foverlaps <- function(bed1,bed2) {
  require(data.table)
  bedKey <- c("CHROM","START","STOP")
  if(nrow(bed1)>nrow(bed2)) {
    bed <- foverlaps(bed1, bed2, nomatch = 0)
  } else {
    bed <- foverlaps(bed2, bed1, nomatch = 0)
  }
  bed[, START := pmax(START, i.START)]
  bed[, STOP := pmin(STOP, i.STOP)]
  bed[, `:=`(i.START = NULL, i.STOP = NULL)]
  if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey)
  if(any(bed[, STOP+1 >= rowShift(START), by=CHROM][,V1], na.rm = T)) {
    bed <- reduceBed.GenomicRanges(bed)
  }
  return(bed)
}

rowShift <- function(x, shiftLen = 1L) {
  #Note this function was described in this thread:
  #http://stackoverflow.com/questions/14689424/use-a-value-from-the-previous-row-in-an-r-data-table-calculation
  r <- (1L + shiftLen):(length(x) + shiftLen)
  r[r<1] <- NA
  return(x[r])
}

reduceBed.GenomicRanges <- function(bed) {
  setnames(bed,colnames(bed),bedKey)
  if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey)
  grBed <- makeGRangesFromDataFrame(bed,
    seqnames.field = "CHROM",start.field="START",end.field="STOP")
  grBed <- reduce(grBed)
  grBed <- data.table(
    CHROM=as.character(seqnames(grBed)),
    START=start(grBed),
    STOP=end(grBed),
    key = c("CHROM","START","STOP"))
  return(grBed)
}

Эта функция строго использовала пакет GenomicRanges, производит тот же результат, но примерно в 10 раз медленнее, чем функция foverlaps.

intersectBedFiles.GR <- function(bed1,bed2) {
  require(data.table)
  require(GenomicRanges)
  bed1 <- makeGRangesFromDataFrame(bed1,
    seqnames.field = "CHROM",start.field="START",end.field="STOP")
  bed2 <- makeGRangesFromDataFrame(bed2,
    seqnames.field = "CHROM",start.field="START",end.field="STOP")
  grMerge <- suppressWarnings(intersect(bed1,bed2))
  resultTable <- data.table(
    CHROM=as.character(seqnames(grMerge)),
    START=start(grMerge),
    STOP=end(grMerge),
    key = c("CHROM","START","STOP"))
  return(resultTable)
}

Дополнительное сравнение с использованием IRanges

Я нашел решение свернуть перекрывающиеся области с помощью IRanges, но он более чем в 10 раз медленнее, чем GenomicRanges.

reduceBed.IRanges <- function(bed) {
  bed.tmp <- bed
  bed.tmp[,group := { 
      ir <-  IRanges(START, STOP);
      subjectHits(findOverlaps(ir, reduce(ir)))
    }, by=CHROM]
  bed.tmp <- bed.tmp[, list(CHROM=unique(CHROM), 
              START=min(START), 
              STOP=max(STOP)),
       by=list(group,CHROM)]
  setkeyv(bed.tmp,bedKey)
  bed[,group := NULL]
  return(bed.tmp[,-c(1:2),with=F])
}


system.time(bedC.reduced <- reduceBed.GenomicRanges(bedC))

user  system elapsed 
10.86    0.01   10.89 

system.time(bedD.reduced <- reduceBed.IRanges(bedC))

user  system elapsed 
137.12    0.14  137.58 

identical(bedC.reduced,bedD.reduced)
[1] TRUE

Ответ 2

foverlaps() будет хорошо.

Сначала установите ключи для обеих таблиц:

setkey(table1, CHROMOSOME, START, STOP)
setkey(table2, CHROMOSOME, START, STOP)

Теперь соедините их с помощью foverlaps() с помощью nomatch = 0, чтобы удалить несогласованные строки в table2.

resultTable <- foverlaps(table1, table2, nomatch = 0)

Затем выберите подходящие значения для START и STOP и оставьте дополнительные столбцы.

resultTable[, START := pmax(START, i.START)]
resultTable[, STOP := pmin(STOP, i.STOP)]
resultTable[, `:=`(i.START = NULL, i.STOP = NULL)]

Наложение STOP на будущий START должно быть другим вопросом. Это на самом деле тот, который у меня есть, поэтому, возможно, я попрошу его и вернусь к нему здесь, когда у меня будет хороший ответ.

Ответ 3

Если вы не застряли в решении data.table, GenomicRanges

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

дает

> library(GenomicRanges)
> intersect(makeGRangesFromDataFrame(table1), makeGRangesFromDataFrame(table2))
GRanges object with 5 ranges and 0 metadata columns:
      seqnames     ranges strand
         <Rle>  <IRanges>  <Rle>
  [1]        1 [  5,  10]      *
  [2]        1 [ 20,  50]      *
  [3]        1 [100, 110]      *
  [4]        1 [130, 130]      *
  [5]        Y [  5,  50]      *
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Ответ 4

В большинстве совпадающих диапазонов проблем в геномике мы имеем один большой набор данных x (обычно секвенированные чтения) и еще один меньший набор данных y (обычно модель гена, экзоны, интроны и т.д.). Нам поручено найти, какие интервалы в x перекрываются с интервалами в y или сколько интервалов в x перекрываются для каждого интервала y.

В foverlaps() нам не нужно setkey() на большом наборе данных x - это довольно дорогостоящая операция. Но y должен иметь ключ. Для вашего случая, из этого примера кажется, что table2 больше = x и table1= y.

require(data.table)
setkey(table1) # key columns = chr, start, end
ans = foverlaps(table2, table1, type="any", nomatch=0L)
ans[, `:=`(i.START = pmax(START, i.START), 
           i.STOP = pmin(STOP, i.STOP))]

ans = ans[, .(i.START[1L], i.STOP[.N]), by=.(CHROMOSOME, START, STOP)]
#    CHROMOSOME START STOP  V1  V2
# 1:          1     1   10   5  10
# 2:          1    20   50  20  50
# 3:          1    70  130 100 130
# 4:          Y     5  200   5  50

Но я согласен, что было бы здорово сделать это за один шаг. Пока не знаете, но, возможно, используя дополнительные значения reduce и intersect для аргумента mult=.

Ответ 5

Здесь решение полностью в data.table на основе ответа Пита. Это на самом деле медленнее, чем его решение, которое использует GenomicRanges и data.table, но все же быстрее, чем решение, которое использует только GenomicRanges.

intersectBedFiles.foverlaps2 <- function(bed1,bed2) {
  require(data.table)
  bedKey <- c("CHROM","START","STOP")
  if(nrow(bed1)>nrow(bed2)) {
    if(!identical(key(bed2),bedKey)) setkeyv(bed2,bedKey)
    bed <- foverlaps(bed1, bed2, nomatch = 0)
  } else {
    if(!identical(key(bed1),bedKey)) setkeyv(bed1,bedKey)
    bed <- foverlaps(bed2, bed1, nomatch = 0)
  }
  bed[,row_id:=1:nrow(bed)]
  bed[, START := pmax(START, i.START)]
  bed[, STOP := pmin(STOP, i.STOP)]
  bed[, ':='(i.START = NULL, i.STOP = NULL)]

  setkeyv(bed,bedKey)
  temp <- foverlaps(bed,bed)

  temp[, ':='(c("START","STOP"),list(min(START,i.START),max(STOP,i.STOP))),by=row_id]
  temp[, ':='(c("START","STOP"),list(min(START,i.START),max(STOP,i.STOP))),by=i.row_id]
  out <- unique(temp[,.(CHROM,START,STOP)])
  setkeyv(out,bedKey)
  out
}