Соединение с перекрытием с начальной и конечной позициями

ProgramBox

Соединение с перекрытием с начальной и конечной позициями

Post author:admin
Запись опубликована:17 февраля, 2023
Post category:Вопросы по программированию

#r #join #merge #data.table

#r #Присоединиться #слияние #данные.таблица

Вопрос:

Рассмотрим следующие data.table s. Первый определяет набор областей с начальной и конечной позициями для каждой группы ‘x’:

 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

Второй набор данных имеет ту же группирующую переменную ‘x’ и позиции ‘pos’ внутри каждой группы:

 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

В конечном счете, я хотел бы извлечь строки в ‘d2’, где ‘pos’ попадает в диапазон, определяемый ‘start’ и ‘end’, внутри каждой группы x . Желаемый результат

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

Начальные / конечные позиции для любой группы x никогда не будут перекрываться, но могут быть пробелы значений не в любом регионе.

Теперь, я полагаю, я должен использовать скользящее соединение. Из того, что я могу сказать, я не могу использовать столбец «end» в соединении.

Я пробовал

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

и получил

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

какой правильный набор строк я хочу; Однако «pos» стал «start», а исходный «start» был потерян. Есть ли способ сохранить все столбцы с объединением по рулону, чтобы я мог сообщить «начало», «pos», «конец» по желанию?

Ответ №1:

Перекрытие соединения было реализовано с фиксацией 1375 в data.table версии v1.9.3 и доступно в текущей стабильной версии v1.9.4. Вызывается функция foverlaps . Из НОВОСТЕЙ:

29) Overlap joins # 528 наконец-то здесь!! За исключением type="equal" maxgap minoverlap аргументов and и, все остальное реализовано. Ознакомьтесь ?foverlaps с примерами его использования. Это важное дополнение к функции data.table .

Давайте рассмотрим x, интервал, определенный как [a, b] , где a <= b , и y, другой интервал, определенный как [c, d] , где c <= d . Говорят, что интервал y перекрывает x вообще, если d >= a и c <= b 1. И y полностью содержится в x, если a <= c,d <= b 2. Пожалуйста, ознакомьтесь с различными типами реализованных перекрытий ?foverlaps .

Ваш вопрос является частным случаем соединения с перекрытием: у d1 вас есть истинные физические интервалы с start позициями и end позициями. С d2 другой стороны, есть только позиции ( pos ), а не интервалы. Чтобы иметь возможность выполнять соединение с перекрытием, нам также нужно создать интервалы d2 . Это достигается путем создания дополнительной переменной pos2 , которая идентична pos ( d2[, pos2 := pos] ) . Таким образом, теперь у нас есть интервал в d2 , хотя и с идентичными начальными и конечными координатами. Затем этот «виртуальный интервал нулевой ширины» d2 можно использовать foverlap для выполнения соединения с перекрытием с 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 по умолчанию это key(y) , поэтому мы пропустили его. by.x по умолчанию принимает key(x) , если оно существует, и если не принимает key(y) . Но ключа для не существует d2 , и мы не можем установить столбцы из y , потому что у них разные имена. Итак, мы устанавливаем by.x явно.

Тип перекрытия находится внутри, и мы хотели бы иметь все совпадения, только если совпадение есть.

ПРИМЕЧАНИЕ: foverlaps использует функцию двоичного поиска data.table (вместе с roll тем, где это необходимо) под капотом, Но некоторые аргументы функции (типы перекрытий, maxgap, minoverlap и т. Д.) Вдохновлены функцией findOverlaps() из пакета Bioconductor IRanges , отличного пакета (и so is GenomicRanges , который распространяется IRanges на геномику).

Так в чем же преимущество?

Тест на приведенном выше коде для ваших данных выполняется foverlaps() медленнее, чем ответ Габора (тайминги: решение Габора data.table = 0,004 против foverlaps = 0,021 секунды). Но действительно ли это имеет значение при такой степени детализации?

Что было бы действительно интересно, так это посмотреть, насколько хорошо это масштабируется — с точки зрения скорости и памяти. В ответе Габора мы объединяемся на основе ключевого столбца x . А затем отфильтруйте результаты.

Что, если d1 имеет около 40 тыс. строк и d2 имеет 100 тыс. строк (или более)? Для каждой строки в d2 , которая соответствует x in d1 , все эти строки будут сопоставлены и возвращены, только для последующей фильтрации. Вот пример вашего Q, масштабируемого незначительно:

Генерировать данные:

 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

Это заняло в общей сложности ~ 1 ГБ памяти, из которых ans1 420 МБ. Большая часть времени, потраченного здесь, на самом деле относится к подмножеству. Вы можете проверить это, установив аргумент verbose=TRUE .

Решения Габора:

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

И это заняло в общей сложности ~ 3,5 ГБ.

Я только что отметил, что Габор уже упоминает память, необходимую для промежуточных результатов. Итак, пробуем 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

Заняло в общей сложности ~ 1,4 ГБ. Таким образом, он определенно использует меньше памяти, чем показано выше.

[Ответы были проверены на идентичность после удаления pos2 из ans1 и установки ключа для обоих ответов.]

Обратите внимание, что это соединение с перекрытием разработано с учетом проблем, когда d2 не обязательно иметь одинаковые начальные и конечные координаты (например: геномика, область, откуда я родом, где d2 обычно около 30-150 миллионов или более строк).

foverlaps() является стабильным, но все еще находится в стадии разработки, что означает, что некоторые аргументы и имена могут быть изменены.

ПРИМЕЧАНИЕ: Поскольку я упоминал GenomicRanges выше, он также вполне способен решить эту проблему. Он использует интервальные деревья под капотом и также довольно эффективен с точки зрения памяти. В моих тестах по данным геномики foverlaps() это быстрее. Но это для другого поста (в блоге), в другой раз.

1. Вы когда-нибудь делали это сообщение в блоге? Я пытаюсь перекрыть почти каждый нуклеотид в геноме человека (оцененный по соответствующей метрике) с набором диапазонов, и я убил findOverlaps GRanges примерно через четыре часа, чтобы начать с касательной foverlaps. Наличие этого готового сообщения в блоге с инструкциями сэкономило бы мне сегодня много времени. : D

Ответ №2:

data.table v1.9.8 имеет новую функцию — неравные соединения. При этом эта операция становится еще более простой:

 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

1. не могли бы вы добавить использование памяти и время для этого подхода

2. Здесь полезно отметить, что x имеет в виду две вещи. Первый x — это столбец «x», а второй x в x.pos относится к d2.

3. @RyanWardValverde Почему это необходимо включить pos=x.pos ? Допустим, в d2 таких столбцах, как ‘a’ и ‘b’, больше столбцов. При выборе столбцов; a и b не обязательно должны иметь префикс x, однако pos имеет. Кроме того, если вы не ставите префикс pos с x тогда pos==start . Вы знаете, почему это ожидается?

4. Похоже на побочный продукт порядка операций. Удаление всего, что находится между скобками в j позиции, т. Е. .(...) Дает нам d2[d1, on = .(x, pos>=start, pos<=end), nomatch = NULL] , и дает столбцы с именами pos и pos.1 , которые соответствуют start и end , что для меня является неожиданным поведением. Это, по-видимому, указывает на то, что соединение выполняется первым, и для доступа к data.table позиции в x позиции его необходимо вызвать явно.

Ответ №3:

1) sqldf Это не data.table, но сложные критерии соединения легко указать прямым способом в SQL:

 library(sqldf)

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

предоставление:

   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 Для ответа data.table попробуйте следующее:

 library(data.table)

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

предоставление:

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

Обратите внимание, что это имеет недостаток в формировании возможно большого промежуточного результата d1[d2] , который SQL может не выполнить. У остальных решений тоже может быть эта проблема.

3) dplyr Это предполагает соответствующее решение dplyr. Мы также используем between from data.table:

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

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

предоставление:

 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) объединение / подмножество с использованием только основания R:

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

предоставление:

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

Добавлено, обратите внимание, что решение таблицы данных здесь намного быстрее, чем в другом ответе:

 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) amp; !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

1. Я ценю вашу всесторонность. Я просто очень хотел, чтобы что-то использовало преимущества оптимизированного поиска по дереву, реализованного в объединении данных.таблицы. Спасибо за sqldf предложение. Я привык писать объединения, как sqldf("select * from d1 join d2 on d1.x==d2.x and d2.pos>=d1.start and d2.pos<=d1.end") в SQL Server, которые были хороши, когда вы могли добавлять индексы в таблицу. Я думаю, что это может избежать создания полного внешнего соединения в памяти (но не тестировалось)

2. Однако обратите внимание, что код data.table здесь короче и быстрее.

3. Вы правы, between это быстрее. Я также тестировал d1[d2, roll=T, nomatch=0, mult="all"][start<=end] с вашим жгутом и был близок к промежуточному. Это просто показывает, что вам нужно все протестировать, потому что вы никогда не знаете, что будет быстрее. Спасибо, что нашли время, чтобы проверить это. Очень интересно. Может быть, когда @Arun реализует «объединение диапазона» в data.table, он сможет сделать это быстрее.

4. Не могли бы вы объяснить, как работает «d1 [d2] [между (pos, start, end)]», пожалуйста?

5. d1[d2] соединяет d1 и d2 вдоль x и [between(...)] выбирает те строки, для которых значение between(...) TRUE .

Ответ №4:

Соединения с перекрытием доступны dplyr 1.1.0 через функцию join_by .

С join_by помощью вы можете выполнить соединение с перекрытием between с помощью или вручную с >= помощью и <= :

 library(dplyr)
inner_join(d2, d1, by = join_by(x, between(pos, start, end)))
#  x pos start end
#1 a   2     1   3
#2 a   3     1   3
#3 c  20    19  22
#4 e  10     7  25

 inner_join(d2, d1, by = join_by(x, pos >= start, pos <= end))
#  x pos start end
#1 a   2     1   3
#2 a   3     1   3
#3 c  20    19  22
#4 e  10     7  25

Ответ №5:

Использование fuzzyjoin :

 result <- fuzzyjoin::fuzzy_inner_join(d1, d2, 
                           by = c('x', 'pos' = 'start', 'pos' = 'end'),
                           match_fun = list(`==`, `>=`, `<=`))
result

#  x.x     pos x.y   start   end
#  <chr> <dbl> <chr> <dbl> <dbl>
#1 a         2 a         1     3
#2 a         3 a         1     3
#3 c        20 c        19    22
#4 e        10 e         7    25

Поскольку fuzzyjoin возвращает все столбцы, нам может потребоваться выполнить некоторую очистку, чтобы сохранить нужные столбцы.

 library(dplyr)
result %>% select(x = x.x, pos, start, end)

# A tibble: 4 x 4
#  x       pos start   end
#  <chr> <dbl> <dbl> <dbl>
#1 a         2     1     3
#2 a         3     1     3
#3 c        20    19    22
#4 e        10     7    25