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

Question

У меня есть два фрейма данных, один называется strain_1, а другой - strain_2. Каждый фрейм данных имеет 4 столбца (st_A, ed_A, st_B, ed_B: для позиций " start " и " end "), но a разное количество рядов . st_A, ed_A и st_B, ed_B - это " start "и позиции « end » block_A и block_B соответственно (см. изображение 1 и приведенный ниже пример).

Я хочу определить общие перекрывающиеся блоки между strain_1 и strain_2.

Пример из изображения 1:

strain_1 <- data.frame(st_A=c(7,25,35,48,89), ed_A=c(9,28,38,51,91),
                       st_B=c(123,97,140,73, 13), ed_B=c(127,98,145,76,16))

strain_2 <- data.frame(st_A=c(5,20,36,49) ,  ed_A=c(8,25,39,50),    
                       st_B=c(124,95,141,105) ,  ed_B=c(129,100,147,110))

Из этого примера мы видим три перекрывающихся области (изображение 1):

Область перекрытия определяется следующим образом: мин значение st_A (или st_B) и максимальное значение ed_A (или ed_B) для block_A и block_B соответственно (см. image 2 : зеленое поле = общая область).

Цель - создать новый data frame с этими общими областями (пара блоков)

## result_desired 
result_desired <- data.frame(st_A=c(5,20,35), ed_A=c(9,28,39),
                             st_B=c(123,95,140), ed_B=c(129,100,147))

Существует 16 возможных комбинаций (см. изображение 3 ), в зависимости от размера каждого блока.

Есть ли быстрый способ сделать это? зная, что у меня есть данные с несколькими тысячами строк.

Я попробовал некоторый код, основанный на комментариях @Gregor, но не могу получить желаемый результат:

## require(dplyr)
require(dplyr)

## data 
strain_1 <- data.frame(st_A=c(7,25,35,48,89), ed_A=c(9,28,38,51,91),
                       st_B=c(123,97,140,73, 13), ed_B=c(127,98,145,76,16))

strain_2 <- data.frame(st_A=c(5,20,36,49) ,  ed_A=c(8,25,39,50),    
                       st_B=c(124,95,141,105) ,  ed_B=c(129,100,147,110))

# merge data to get cross join 
cj_data <-merge(strain_1,strain_2, by = NULL)

# Check block1 and block2
cj_filtered <- cj_data %>% mutate(c_block1= case_when(st_A.x <= st_A.y & ed_A.x <= ed_A.y | 
                                                      st_A.x >= st_A.y & ed_A.x >= ed_A.y |
                                                      st_A.x <= st_A.y & ed_A.x >= ed_A.y | 
                                                      st_A.x >= st_A.y & ed_A.x <= ed_A.y ~ "overlap_OK",
                                                      TRUE ~ "NO"),

                                  c_block2= case_when(st_B.x <= st_B.y & ed_B.x <= ed_B.y | 
                                                      st_B.x >= st_B.y & ed_B.x >= ed_B.y |
                                                      st_B.x <= st_B.y & ed_B.x >= ed_B.y | 
                                                      st_B.x >= st_B.y & ed_B.x <= ed_B.y ~ "overlap_OK",
                                                      TRUE ~ "NO"))

## cj_filtered:
st_A.x  ed_A.x  st_B.x  ed_B.x  st_A.y  ed_A.y  st_B.y  ed_B.y  c_block1    c_block2
7       9       123     127     5       8       124     129     overlap_OK  overlap_OK
25      28      97      98      5       8       124     129     overlap_OK  overlap_OK
35      38      140     145     5       8       124     129     overlap_OK  overlap_OK
48      51      73      76      5       8       124     129     overlap_OK  overlap_OK
89      91      13      16      5       8       124     129     overlap_OK  overlap_OK
7       9       123     127     20      25      95      100     overlap_OK  overlap_OK
25      28      97      98      20      25      95      100     overlap_OK  overlap_OK
35      38      140     145     20      25      95      100     overlap_OK  overlap_OK
48      51      73      76      20      25      95      100     overlap_OK  overlap_OK
89      91      13      16      20      25      95      100     overlap_OK  overlap_OK
7       9       123     127     36      39      141     147     overlap_OK  overlap_OK
25      28      97      98      36      39      141     147     overlap_OK  overlap_OK
35      38      140     145     36      39      141     147     overlap_OK  overlap_OK
48      51      73      76      36      39      141     147     overlap_OK  overlap_OK
89      91      13      16      36      39      141     147     overlap_OK  overlap_OK
7       9       123     127     49      50      105     110     overlap_OK  overlap_OK
25      28      97      98      49      50      105     110     overlap_OK  overlap_OK
35      38      140     145     49      50      105     110     overlap_OK  overlap_OK
48      51      73      76      49      50      105     110     overlap_OK  overlap_OK
89      91      13      16      49      50      105     110     overlap_OK  overlap_OK

Спасибо за вашу помощь.

Answer 1

Вот 2 варианта использования data.table:

1a) Использование неэквивалентных объединений:

cols <- c(paste0("x.", names(strain_1)), paste0("i.", names(strain_2)))
DT <- rbindlist(list(
    strain_1[strain_2, on=.(st_A>=st_A, st_A<=ed_A), nomatch=0L, mget(cols)],
    strain_1[strain_2, on=.(st_A<=st_A, ed_A>=st_A), nomatch=0L, mget(cols)]
))

1b) Использование foverlaps:

setkey(strain_1, st_A, ed_A)
setkey(strain_2, st_A, ed_A)
foverlaps(strain_1, strain_2, nomatch=0L)

И затем еще один шаг 2, чтобы получить желаемый вывод:

DT[between(x.st_B, i.st_B, i.ed_B) | between(i.st_B, x.st_B, x.ed_B), 
    .(st_A=pmin(x.st_A, i.st_A), 
        ed_A=pmax(x.ed_A, i.ed_A), 
        st_B=pmin(x.st_B, i.st_B), 
        ed_B=pmax(x.ed_B, i.ed_B))]

вывод:

   st_A ed_A st_B ed_B
1:    5    9  123  129
2:   20   28   95  100
3:   35   39  140  147

data:

library(data.table)
strain_1 <- data.frame(st_A=c(7,25,35,48,89), ed_A=c(9,28,38,51,91),
    st_B=c(123,97,140,73, 13), ed_B=c(127,98,145,76,16))

strain_2 <- data.frame(st_A=c(5,20,36,49) ,  ed_A=c(8,25,39,50),    
    st_B=c(124,95,141,105) ,  ed_B=c(129,100,147,110))

result_desired <- data.frame(st_A=c(5,20,35), ed_A=c(9,28,39),
    st_B=c(123,95,140), ed_B=c(129,100,147))

setDT(strain_1)
setDT(strain_2)
setDT(result_desired)

ps: в Bioconducter должно быть что-тос IRanges.