¿Es posible usar la función R data.table foverlaps para encontrar la intersección de los rangos superpuestos en dos tablas?

Me gustaría utilizar foverlaps para encontrar los rangos de intersección de dos archivos de cama, y ​​contraer las filas que contienen rangos superpuestos en una sola fila. En el siguiente ejemplo, tengo dos tablas con rangos genómicos. Las tablas se denominan archivos de “cama” que tienen coordenadas de inicio basadas en cero y posiciones finales basadas en uno de las características de los cromosomas. Por ejemplo, START = 9, STOP = 20 se interpreta para abarcar las bases 10 a 20, inclusive. Estos archivos de cama pueden contener millones de filas. La solución debería dar el mismo resultado, independientemente del orden en que se proporcionen los dos archivos que se intersectarán.

Primera mesa

> table1 CHROMOSOME START STOP 1: 1 1 10 2: 1 20 50 3: 1 70 130 4: X 1 20 5: Y 5 200 

Segunda mesa

 > table2 CHROMOSOME START STOP 1: 1 5 12 2: 1 15 55 3: 1 60 65 4: 1 100 110 5: 1 130 131 6: X 60 80 7: Y 1 15 8: Y 10 50 

Estaba pensando que la nueva función foverlaps podría ser una manera muy rápida de encontrar los intervalos de intersección en estas dos tablas para producir una tabla que se vería así:

Tabla de resultados:

 > resultTable CHROMOSOME START STOP 1: 1 5 10 2: 1 20 50 3: 1 100 110 4: Y 5 50 

¿Es eso posible, o hay una mejor manera de hacerlo en data.table?

También me gustaría primero confirmar que en una tabla, para cualquier CHROMOSOME dado, la coordenada STOP no se superpone con la coordenada inicial de la siguiente fila. Por ejemplo, CROMOSOMA Y: 1-15 y CROMOSOMO Y: 10-50 necesitarían colapsarse para CROMOSOMO Y: 1-50 (ver la Segunda Tabla Filas 7 y 8). Este no debería ser el caso, pero la función probablemente debería verificar eso. Un ejemplo de la vida real de cómo las superposiciones potenciales deberían colapsarse es a continuación:

CHROM START STOP 1: 1 721281 721619 2: 1 721430 721906 3: 1 721751 722042

Salida deseada:

CHROM START STOP 1: 1 721281 722042

Las funciones para crear tablas de ejemplo son las siguientes:

 table1 <- data.table( CHROMOSOME = as.character(c("1","1","1","X","Y")) , START = c(1,20,70,1,5) , STOP = c(10,50,130,20,200) ) table2 <- data.table( CHROMOSOME = as.character(c("1","1","1","1","1","X","Y","Y")) , START = c(5,15,60,100,130,60,1,10) , STOP = c(12,55,65,110,131,80,15,50) ) 

@Seth brindó la forma más rápida de resolver el problema de las superposiciones de intersecciones utilizando la función data.table foverlaps. Sin embargo, esta solución no tuvo en cuenta el hecho de que los archivos de lecho de entrada pueden tener rangos superpuestos que deben reducirse a regiones individuales. @Martin Morgan resolvió eso con su solución usando el paquete GenomicRanges, que hizo tanto la intersección como la reducción de rango. Sin embargo, la solución de Martin no usó la función foverlaps. @Arun señaló que los rangos superpuestos en diferentes filas dentro de una tabla no eran posibles en la actualidad usando foverlaps. Gracias a las respuestas proporcionadas y algunas investigaciones adicionales sobre stackoverflow, se me ocurrió esta solución híbrida.

Cree archivos BED de ejemplo sin superponer regiones dentro de cada archivo.

 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") ) 

Ahora crea un nuevo archivo de cama que contenga las regiones de intersección en bedA y 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 

Ahora, modifique bedA y bedB de manera que contengan regiones superpuestas:

 #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] 

El tiempo de prueba para intersectar los archivos de la cama con rangos superpuestos usando las funciones foverlaps o 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 

El ganador: foverlaps !

FUNCIONES UTILIZADAS

Esta es la función basada en foverlaps, y solo llamará a la función GenomicRanges (reduceBed.GenomicRanges) si hay rangos superpuestos (que están marcados para usar la función 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) } 

Esta función utiliza estrictamente el paquete GenomicRanges, produce el mismo resultado, pero es aproximadamente 10 veces más lenta que la función 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) } 

Una comparación adicional usando IRanges

Encontré una solución para colapsar regiones superpuestas usando IRanges, pero es más de 10 veces más lenta que 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 

foverlaps() hará muy bien.

Primero establece las claves para ambas tablas:

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

Ahora únase a ellos usando foverlaps() con nomatch = 0 para eliminar las filas no table2 en table2 .

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

A continuación, elija los valores adecuados para START y STOP, y suelte las columnas adicionales.

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

La PARADA superpuesta a un futuro START debería ser una pregunta diferente. En realidad, es uno que tengo, así que tal vez lo pregunte y vuelvo aquí cuando tenga una buena respuesta.

En caso de que no estés atascado en una solución de data.table , GenomicRanges

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

da

 > library(GenomicRanges) > intersect(makeGRangesFromDataFrame(table1), makeGRangesFromDataFrame(table2)) GRanges object with 5 ranges and 0 metadata columns: seqnames ranges strand    [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 

En la mayoría de los rangos de superposición de problemas en genómica, tenemos un gran conjunto de datos x (generalmente lecturas secuenciadas) y otro conjunto de datos más pequeño (por lo general, el modelo de gen, exones, intrones, etc.). Nos corresponde buscar los intervalos en x superposición con qué intervalos en y o cuántos intervalos en x superponen para cada intervalo y .

En foverlaps() , no tenemos que setkey() en el conjunto de datos más grande x – es una operación bastante costosa. Pero necesita tener su llave establecida. Para su caso, de este ejemplo, parece que table2 es más grande = x , y 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 

Pero acepto que sería genial poder hacer esto en un solo paso. Aún no estoy seguro de cómo, pero tal vez usar valores adicionales reduce e intersect para mult= argumento.