Count the number of overlaps between groups

Question

I have two large data sets that look like this.

library(tidyverse)


dat1 <- tibble(chrom=c(rep(c("Chr1","Chr2"),each=5)),
               start=c(9885,11944, 13271,15104,19059,25793,97514,104718,118862,120950),
               end=c(11008,17644,20164,23807,25264,106001,119205, 121576,124981,138514)
)

head(dat1,n=4)
#> # A tibble: 10 × 3
#>    chrom  start    end
#>    <chr>  <dbl>  <dbl>
#>  1 Chr1    9885  11008
#>  2 Chr1   11944  17644
#>  3 Chr1   13271  20164
#>  4 Chr1   15104  23807


dat2 <- tibble(chrom=c(rep(c("Chr1","Chr2"),each=5)),
               start=c(9885,11944, 13271,15104,19059,25793,97514,104718,118862,120950),
               end=c(10203,12546,13669,15638,19283,26703,97773, 105102,119388,121331)
               )

head(dat2, n=4)
#> # A tibble: 10 × 3
#>    chrom  start    end
#>    <chr>  <dbl>  <dbl>
#>  1 Chr1    9885  10203
#>  2 Chr1   11944  12546
#>  3 Chr1   13271  13669
#>  4 Chr1   15104  15638

^{Created on 2022-12-05 with reprex v2.0.2}

I want to group my data based on the chrom and find which ranges [start-end], from the Chr1, from dat2 overlap with ranges [start-end] of Chr1 from dat1.

What I have tried

I have found a lovely package to work it around, but I feel I need to break the data sets into data frames of different chromosomes to compare them.

library(plyranges)

dat1 <- dat1 %>% 
  as_iranges()

dat2  <- dat2 %>% 
  as_iranges()

dat1 %>% 
  mutate(n_olap = count_overlaps(., dat2),
                 n_olap_within = count_overlaps_within(., dat2))

IRanges object with 10 ranges and 3 metadata columns:
           start       end     width |       chrom    n_olap n_olap_within
       <integer> <integer> <integer> | <character> <integer>     <integer>
   [1]      9885     11008      1124 |        Chr1         1             0
   [2]     11944     17644      5701 |        Chr1         3             0
   [3]     13271     20164      6894 |        Chr1         3             0
   [4]     15104     23807      8704 |        Chr1         2             0
   [5]     19059     25264      6206 |        Chr1         1             0
   [6]     25793    106001     80209 |        Chr2         3             0
   [7]     97514    119205     21692 |        Chr2         3             0
   [8]    104718    121576     16859 |        Chr2         3             0
   [9]    118862    124981      6120 |        Chr2         2             0
  [10]    120950    138514     17565 |        Chr2         1             0

To get what I want from here, I need to filter my data and then compare it. But there should a way or dplyr trick to find a solution

dat1 <- dat1 %>% 
  as_iranges() %>% 
  filter(chrom=="Chr1")

dat2  <- dat2 %>% 
  as_iranges() %>% 
  filter(chrom=="Chr1")

dat1 %>% 
  mutate(n_olap = count_overlaps(., dat2),
                 n_olap_within = count_overlaps_within(., dat2))

Is there a way to compare only chromosomes with each other?

Answer 1

I'll stick with the data.table theme set forth by RicVillalba, but I think the foverlaps function is meant for things like this (especially with larger datasets).

library(data.table)
setDT(dat1)
setDT(dat2)
setkey(dat1, chrom, start, end)
setkey(dat2, chrom, start, end)
dat1[, id := .I]
foverlaps(dat1, dat2)
#      chrom  start    end i.start  i.end    id
#     <char>  <num>  <num>   <num>  <num> <int>
#  1:   Chr1   9885  10203    9885  11008     1
#  2:   Chr1  11944  12546   11944  17644     2
#  3:   Chr1  13271  13669   11944  17644     2
#  4:   Chr1  15104  15638   11944  17644     2
#  5:   Chr1  13271  13669   13271  20164     3
#  6:   Chr1  15104  15638   13271  20164     3
#  7:   Chr1  19059  19283   13271  20164     3
#  8:   Chr1  15104  15638   15104  23807     4
#  9:   Chr1  19059  19283   15104  23807     4
# 10:   Chr1  19059  19283   19059  25264     5
# ---                                          
# 13:   Chr2 104718 105102   25793 106001     6
# 14:   Chr2  97514  97773   97514 119205     7
# 15:   Chr2 104718 105102   97514 119205     7
# 16:   Chr2 118862 119388   97514 119205     7
# 17:   Chr2 104718 105102  104718 121576     8
# 18:   Chr2 118862 119388  104718 121576     8
# 19:   Chr2 120950 121331  104718 121576     8
# 20:   Chr2 118862 119388  118862 124981     9
# 21:   Chr2 120950 121331  118862 124981     9
# 22:   Chr2 120950 121331  120950 138514    10

(Note that in addition to requiring keys, the order matters: the last two keys must be the "start" and "end" of the ranges to check for overlapping. I added chrom to ensure we do this by-chromosome.)

That's the start. I added the id column into dat1 so that we could efficiently return back to the original columns. If you inspect the columns closely, notice that the i.* columns are from dat1, so those are the ones we want to keep.

Extending that to do the aggregation you're hoping for,

overlaps <- foverlaps(dat1, dat2)[, .(n_olaps = .N, n_within = sum(between(i.start, start, end) & between(i.end, start, end))), by = .(id)]
overlaps
#        id n_olaps n_within
#     <int>   <int>    <int>
#  1:     1       1        0
#  2:     2       3        0
#  3:     3       3        0
#  4:     4       2        0
#  5:     5       1        0
#  6:     6       3        0
#  7:     7       3        0
#  8:     8       3        0
#  9:     9       2        0
# 10:    10       1        0

dat1 <- overlaps[dat1, on = .(id)]
dat1
#        id n_olaps n_within  chrom  start    end
#     <int>   <int>    <int> <char>  <num>  <num>
#  1:     1       1        0   Chr1   9885  11008
#  2:     2       3        0   Chr1  11944  17644
#  3:     3       3        0   Chr1  13271  20164
#  4:     4       2        0   Chr1  15104  23807
#  5:     5       1        0   Chr1  19059  25264
#  6:     6       3        0   Chr2  25793 106001
#  7:     7       3        0   Chr2  97514 119205
#  8:     8       3        0   Chr2 104718 121576
#  9:     9       2        0   Chr2 118862 124981
# 10:    10       1        0   Chr2 120950 138514

I chose to create overlaps and join it back onto dat1 in case there were other columns that needed to be retained. It's certainly feasible to do this more in-place in dat1, especially helpful if your dataset is too big for a temporary copy.

Answer 2

If you have huge datasets, you can consider using data.table this way:

library(data.table)

dat1 <- data.frame(chrom=c(rep(c("Chr1","Chr2"),each=5)),
               start=c(9885,11944, 13271,15104,19059,25793,97514,104718,118862,120950),
               end=c(11008,17644,20164,23807,25264,106001,119205, 121576,124981,138514)
)

dat2 <- data.frame(chrom=c(rep(c("Chr1","Chr2"),each=5)),
               start=c(9885,11944, 13271,15104,19059,25793,97514,104718,118862,120950),
               end=c(10203,12546,13669,15638,19283,26703,97773, 105102,119388,121331)
               )

setDT(dat1)
setDT(dat2)

dat1[dat2, cbind(
  .SD[i.end >= start & end >= i.start],
  start2 = i.start,
  end2 = i.end), on="chrom", by=.EACHI][!is.na(start)]
#>     chrom  start    end start2   end2
#>  1:  Chr1   9885  11008   9885  10203
#>  2:  Chr1  11944  17644  11944  12546
#>  3:  Chr1  13271  20164  13271  13669
#>  4:  Chr1  15104  23807  15104  15638
#>  5:  Chr1  19059  25264  19059  19283
#>  6:  Chr2  25793 106001  25793  26703
#>  7:  Chr2  97514 119205  97514  97773
#>  8:  Chr2 104718 121576 104718 105102
#>  9:  Chr2  97514 119205 118862 119388
#> 10:  Chr2 118862 124981 118862 119388
#> 11:  Chr2 120950 138514 120950 121331

EDIT i made a correction on overlap matching. Thanks to r2evans.