roll join con ventana de inicio/final

merge data.table

Question

Considere los siguientes data.table s. El primero define un conjunto de regiones con posiciones de inicio y final para cada grupo

library(data.table)

d1 <- data.table(x=letters[1:5], start=c(1,5,19,30, 7), end=c(3,11,22,39,25))
setkey(d1, x,start)

#    x start end
# 1: a     1   3
# 2: b     5  11
# 3: c    19  22
# 4: d    30  39
# 5: e     7  25

El segundo representa observaciones para cada grupo

d2 <- data.table(x=letters[c(1,1,2,2,3:5)], pos=c(2,3,3,12,20,52,10))
setkey(d2, x,pos)

#    x pos
# 1: a   2
# 2: a   3
# 3: b   3
# 4: b  12
# 5: c  20
# 6: d  52
# 7: e  10

En última instancia, me gustaría poder extraer las filas en d2 que están en una región para el valor x correspondiente en d1. El resultado deseado es

#    x pos start  end
# 1: a   2     1    3
# 2: a   3     1    3
# 3: c  20    19   22
# 4: e  10     7   25

Las posiciones de inicio / finalización para cualquier grupo x nunca se superpondrán, pero puede haber espacios vacíos de valores que no se encuentran en ninguna región.

Ahora, creo que debería usar una unión progresiva. Por lo que puedo decir, no puedo usar la columna "final" en la unión.

He intentado

d1[d2, roll=T, nomatch=0, mult="all"][start<=end]

y consiguió

#    x start end
# 1: a     2   3
# 2: a     3   3
# 3: c    20  22
# 4: e    10  25

cuál es el conjunto correcto de filas que quiero; Sin embargo, "pos" se ha convertido en "inicio" y el "inicio" original se ha perdido. ¿Hay alguna manera de preservar todas las columnas con la combinación de rollo para que pueda informar "inicio", "pos", "final" como desee?

Answer 1

Las combinaciones de superposición se implementaron con la confirmación 1375 en data.table v1.9.3 , y está disponible en la versión estable actual, v1.9.4 . La función se llama foverlaps . De NEWS :

29) ¡La Overlap joins # 528 ahora está aquí, finalmente! A excepción de los argumentos type="equal" y maxgap y minoverlap , todo lo demás se implementa. Echa un vistazo a ?foverlaps y los ejemplos allí en su uso. Esta es una característica principal además de data.table .

Consideremos x, un intervalo definido como [a, b] , donde a <= b , y y, otro intervalo definido como [c, d] , donde c <= d . Se dice que el intervalo y solapa a x, iff d >= a y c <= b 1 . Y y está completamente contenida dentro de x, si a <= c,d <= b 2 . Para los diferentes tipos de superposiciones implementadas, eche un vistazo a ?foverlaps .

Su pregunta es un caso especial de una combinación de solapamiento: en d1 tiene intervalos físicos verdaderos con posiciones start y end . En d2 por otro lado, solo hay posiciones ( pos ), no intervalos. Para poder hacer una unión por solapamiento, también necesitamos crear intervalos en d2 . Esto se logra creando una variable adicional pos2 , que es idéntica a pos ( d2[, pos2 := pos] ). Por lo tanto, ahora tenemos un intervalo en d2 , aunque con idénticas coordenadas de inicio y fin . Este ''intervalo virtual de ancho cero'' en d2 se puede usar en foverlap para hacer una combinación de superposición con d1 :

require(data.table) ## 1.9.3
setkey(d1)
d2[, pos2 := pos]
foverlaps(d2, d1, by.x = names(d2), type = "within", mult = "all", nomatch = 0L)
#    x start end pos pos2
# 1: a     1   3   2    2
# 2: a     1   3   3    3
# 3: c    19  22  20   20
# 4: e     7  25  10   10

by.y de forma predeterminada es la key(y) , por lo que se salteó. by.x toma por defecto la key(x) si existe, y si no toma la key(y) . Pero no existe una clave para d2 , y no podemos establecer las columnas desde y , porque no tienen los mismos nombres. Entonces, establecemos por by.x explícitamente.

El tipo de superposición está dentro , y nos gustaría tener todas las coincidencias, solo si hay una coincidencia.

NB: foverlaps usa la función de búsqueda binaria de data.table (junto con el roll donde sea necesario) bajo el capó, pero algunos argumentos de función (tipos de solapamientos, maxgap, minoverlap, etc.) están inspirados en la función findOverlaps() del paquete Bioconductor 2 , un paquete excelente (y también lo es GenomicRanges , que extiende IRanges para Genomics).

Entonces, ¿cuál es la ventaja?

Un punto de referencia en el código anterior en sus datos da como resultado foverlaps() más lento que la respuesta de Gabor (Tiempos: la solución data.table de Gabor = 0.004 vs foverlaps = 0.021 segundos). Pero, ¿realmente importa en esta granularidad?

Lo que sería realmente interesante es ver qué tan bien se escala, en términos de velocidad y memoria . En la respuesta de Gabor, nos unimos en base a la columna clave x . Y luego filtra los resultados.

¿Qué pasa si d1 tiene aproximadamente 40K filas y d2 tiene 100K filas (o más)? Para cada fila en d2 que coincida con x en d1 , todas esas filas se emparejarán y se devolverán, para luego ser filtradas. Aquí hay un ejemplo de su Q escalada solo ligeramente:

Generar datos:

require(data.table)
set.seed(1L)
n = 20e3L; k = 100e3L
idx1 = sample(100, n, TRUE)
idx2 = sample(100, n, TRUE)
d1 = data.table(x = sample(letters[1:5], n, TRUE), 
                start = pmin(idx1, idx2), 
                end = pmax(idx1, idx2))

d2 = data.table(x = sample(letters[1:15], k, TRUE), 
                pos1 = sample(60:150, k, TRUE))

foverlaps:

system.time({
    setkey(d1)
    d2[, pos2 := pos1]
    ans1 = foverlaps(d2, d1, by.x=1:3, type="within", nomatch=0L)
})
# user  system elapsed 
#   3.028   0.635   3.745 

Esto tomó ~ 1GB de memoria en total, de los cuales ans1 es 420MB. La mayor parte del tiempo que pasamos aquí está realmente en un subconjunto. Puede verificarlo estableciendo el argumento verbose=TRUE .

Soluciones de Gabor:

## new session - data.table solution
system.time({
    setkey(d1, x)
    ans2 <- d1[d2, allow.cartesian=TRUE, nomatch=0L][between(pos1, start, end)]
})
#   user  system elapsed 
# 15.714   4.424  20.324 

Y esto tomó un total de ~ 3.5GB.

Solo noté que Gabor ya menciona la memoria requerida para los resultados intermedios. Entonces, probando sqldf :

# new session - sqldf solution
system.time(ans3 <- sqldf("select * from d1 join 
            d2 using (x) where pos1 between start and end"))
#   user  system elapsed 
# 73.955   1.605  77.049 

Tomó un total de ~ 1.4GB. Por lo tanto, definitivamente usa menos memoria que la que se muestra arriba.

[Se verificó que las respuestas son idénticas después de eliminar pos2 de ans1 y establecer la clave en ambas respuestas.]

Tenga en cuenta que esta combinación de superposición está diseñada con problemas donde d2 no necesariamente tiene coordenadas de inicio y final idénticas (por ejemplo, genómica, el campo de donde vengo, donde d2 suele ser de 30-150 millones o más filas).

foverlaps() es estable, pero aún está en desarrollo, lo que significa que algunos argumentos y nombres pueden cambiarse.

NB: Desde que mencioné GenomicRanges arriba, también es perfectamente capaz de resolver este problema. Utiliza árboles de intervalo debajo del capó, y también es bastante eficiente en cuanto a la memoria. En mis puntos de referencia sobre datos de genómica, foverlaps() es más rápido. Pero eso es para otra publicación (blog), en otro momento.

Answer 2

data.table v1.9.8+ tiene una nueva característica - no coincide . Con eso, esta operación se vuelve aún más directa:

require(data.table) #v1.9.8+
# no need to set keys on `d1` or `d2`
d2[d1, .(x, pos=x.pos, start, end), on=.(x, pos>=start, pos<=end), nomatch=0L]
#    x pos start end
# 1: a   2     1   3
# 2: a   3     1   3
# 3: c  20    19  22
# 4: e  10     7  25

Answer 3

1) sqldf Esto no es data.table, pero los criterios de combinación complejos son fáciles de especificar de manera directa en SQL:

library(sqldf)

sqldf("select * from d1 join d2 using (x) where pos between start and end")

dando:

  x start end pos
1 a     1   3   2
2 a     1   3   3
3 c    19  22  20
4 e     7  25  10

2) data.table Para una respuesta de data.table prueba esto:

library(data.table)

setkey(d1, x)
setkey(d2, x)
d1[d2][between(pos, start, end)]

dando:

   x start end pos
1: a     1   3   2
2: a     1   3   3
3: c    19  22  20
4: e     7  25  10

Tenga en cuenta que esto tiene la desventaja de formar el resultado intermeidate posiblemente grande d1[d2] que SQL puede no hacer. Las soluciones restantes pueden tener este problema también.

3) dplyr Esto sugiere la solución dplyr correspondiente. También usamos between data.table:

library(dplyr)
library(data.table) # between

d1 %>% 
   inner_join(d2) %>% 
   filter(between(pos, start, end))

dando:

Joining by: "x"
  x start end pos
1 a     1   3   2
2 a     1   3   3
3 c    19  22  20
4 e     7  25  10

4) fusionar / subconjunto Usar solo la base de R:

subset(merge(d1, d2), start <= pos & pos <= end)

dando:

   x start end pos
1: a     1   3   2
2: a     1   3   3
3: c    19  22  20
4: e     7  25  10

Agregado Tenga en cuenta que la solución de la tabla de datos aquí es mucho más rápida que la de la otra respuesta:

dt1 <- function() {
 d1 <- data.table(x=letters[1:5], start=c(1,5,19,30, 7), end=c(3,11,22,39,25))
 d2 <- data.table(x=letters[c(1,1,2,2,3:5)], pos=c(2,3,3,12,20,52,10))
 setkey(d1, x, start)
 idx1 = d1[d2, which=TRUE, roll=Inf] # last observation carried forwards

 setkey(d1, x, end)
 idx2 = d1[d2, which=TRUE, roll=-Inf] # next observation carried backwards

 idx = which(!is.na(idx1) & !is.na(idx2))
 ans1 <<- cbind(d1[idx1[idx]], d2[idx, list(pos)])
}

dt2 <- function() {
 d1 <- data.table(x=letters[1:5], start=c(1,5,19,30, 7), end=c(3,11,22,39,25))
 d2 <- data.table(x=letters[c(1,1,2,2,3:5)], pos=c(2,3,3,12,20,52,10))
 setkey(d1, x)
 ans2 <<- d1[d2][between(pos, start, end)]
}

all.equal(as.data.frame(ans1), as.data.frame(ans2))
## TRUE

benchmark(dt1(), dt2())[1:4]
##     test replications elapsed relative
##  1 dt1()          100    1.45    1.667  
##  2 dt2()          100    0.87    1.000  <-- from (2) above