roll join con start / end window

Considera i seguenti data.table s. Il primo definisce un insieme di regioni con posizioni iniziali e finali per ciascun gruppo

 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 

Il secondo rappresenta le osservazioni per ciascun gruppo

 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 

In definitiva mi piacerebbe essere in grado di estrarre le righe in d2 che si trovano in una regione per il valore x corrispondente in d1. Il risultato desiderato è

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

Le posizioni di inizio / fine per qualsiasi gruppo x non si sovrappongono mai, ma potrebbero esserci lacune di valori non in nessuna regione.

Ora, credo che dovrei usare un rolling join. Da quello che posso dire, non posso usare la colonna “fine” nel join.

ho provato

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

e ottenuto

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

che è l’insieme giusto di righe che voglio; Tuttavia “pos” è diventato “start” e l’originale “start” è stato perso. C’è un modo per preservare tutte le colonne con il roll join in modo da poter segnalare “start”, “pos”, “end” come desiderato?

Overlap joins è stato implementato con commit 1375 in data.table v1.9.3 , ed è disponibile nella versione stabile corrente, v1.9.4 . La funzione è chiamata foverlaps . Dalle NOTIZIE :

29) Overlap joins # 528 ora è qui, finalmente !! Tranne che per gli argomenti type="equal" e maxgap e minoverlap , tutto il resto è implementato. Guarda ?foverlaps e gli esempi lì sul suo utilizzo. Questa è una caratteristica importante aggiunta a data.table .

Consideriamo x, un intervallo definito come [a, b] , dove a <= b e y, un altro intervallo definito come [c, d] , dove c <= d . Si dice che l'intervallo y si sovrapponga x, iff d >= a e c <= b 1 . E y è interamente contenuto in x, iff a <= c,d <= b 2 . Per i diversi tipi di sovrapposizioni implementate, si prega di dare un'occhiata a ?foverlaps .

La tua domanda è un caso particolare di un join di sovrapposizione: in d1 hai veri intervalli fisici con posizioni start e end . D'altra parte in d2 ci sono solo posizioni ( pos ), non intervalli. Per poter eseguire un join di sovrapposizione, è necessario creare intervalli anche in d2 . Ciò si ottiene creando una variabile addizionale pos2 , che è identica a pos ( d2[, pos2 := pos] ). Quindi, ora abbiamo un intervallo in d2 , anche se con coordinate iniziali e finali identiche. Questo 'intervallo virtuale, a larghezza zero' in d2 può quindi essere usato in foverlap per fare un join di sovrapposizione 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 di default è la key(y) , quindi l'abbiamo saltata. by.x per default prende la key(x) se esiste, e se non prende la key(y) . Ma una chiave non esiste per d2 e non possiamo impostare le colonne da y , perché non hanno lo stesso nome. Quindi, impostiamo by.x esplicitamente.

Il tipo di sovrapposizione è all'interno e vorremmo avere tutte le partite, solo se c'è una corrispondenza.

NB: foverlaps utilizza la funzione di ricerca binaria data.table (insieme al roll dove necessario) sotto il cofano, ma alcuni argomenti di funzione (tipi di sovrapposizioni, maxgap, minoverlap ecc.) Sono ispirati alla funzione findOverlaps() dal pacchetto IRanges , un pacchetto eccellente (così come GenomicRanges , che estende IRanges per Genomics).


Quindi qual è il vantaggio?

Un punto di riferimento sul codice sopra riportato sui dati risulta in foverlaps() più lento della risposta di Gabor (Tempi: Gabor's data.table solution = 0,004 vs foverlaps = 0,021 secondi). Ma importa davvero a questa granularità?

Quello che sarebbe davvero interessante è vedere come si scala - in termini di velocità e memoria . Nella risposta di Gabor, ci uniamo in base alla colonna chiave x . E quindi filtrare i risultati.

Cosa succede se d1 ha circa 40.000 righe e d2 ha 100.000 righe (o più)? Per ogni riga in d2 che corrisponde a x in d1 , tutte le righe verranno confrontate e restituite, solo per essere filtrate in seguito. Ecco un esempio del tuo Q ridimensionato solo leggermente:

Genera dati:

 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 

Ci sono voluti circa 1 GB di memoria in totale, di cui ans1 è 420MB. La maggior parte del tempo trascorso qui è sottoinsieme davvero. Puoi controllarlo impostando l'argomento verbose=TRUE .

Le soluzioni di 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 

E questo ha richiesto un totale di ~ 3,5 GB.

Ho appena notato che Gabor menziona già la memoria richiesta per risultati intermedi. Quindi, provando 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 

Ha richiesto un totale di ~ 1,4 GB. Quindi, usa decisamente meno memoria di quella mostrata sopra.

[Le risposte sono state verificate per essere identiche dopo aver rimosso pos2 da ans1 e impostato la chiave su entrambe le risposte.]

Si noti che questa unione di sovrapposizione è progettata con problemi in cui d2 non ha necessariamente coordinate iniziali e finali identiche (es: genomica, il campo da cui provengo, dove d2 è solitamente circa 30-150 milioni o più righe).


foverlaps() è stabile, ma è ancora in fase di sviluppo, il che significa che alcuni argomenti e nomi potrebbero essere cambiati.

NB: Da quando ho menzionato GenomicRanges sopra, è anche perfettamente in grado di risolvere questo problema. Usa gli alberi ad intervalli sotto il cofano, ed è anche abbastanza efficiente per la memoria. Nei miei benchmark sui dati di genomica, foverlaps() è più veloce. Ma questo è per un altro post (blog), un'altra volta.

1) sqldf Questo non è data.table ma i criteri di join complessi sono facili da specificare in modo diretto in 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 Per una risposta data.table prova questo:

 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 

Si noti che questo ha lo svantaggio di formare il risultato intermeidato possibilmente grande d1[d2] che SQL non può fare. Anche le restanti soluzioni potrebbero avere questo problema.

3) dplyr Questo suggerisce la corrispondente soluzione dplyr. Usiamo anche da 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) unione / sottoinsieme usando solo la base di 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 

Aggiunto Si noti che la soluzione della tabella di dati qui è molto più veloce di quella nell'altra risposta:

 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 

data.table v1.9.8+ ha una nuova funzionalità – join non equi . Con ciò, questa operazione diventa ancora più semplice:

 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