1
$\begingroup$

I want to retrieve all overlaps between two sets of genomic intervals (GRanges) using Bioconductor. I am interested in the widths of the resulting overlapping ranges.

In cases with discontinuous ranges, intersect() works, but I am interested in the case where both sets contain contiguous ranges (no gaps).

I am using the findOverlaps function as in the example below:

gr1 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 11, 21, 31), c(10, 20, 30, 40)))
gr2 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 5, 15), c(4, 14, 30)))
hits <- findOverlaps(gr1, gr2)

This code yields the following intersections in a "Hits" object:

Hits object with 5 hits and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  [2]         1           2
  [3]         2           2
  [4]         2           3
  [5]         3           3
  -------
  queryLength: 4 / subjectLength: 3

But I could not find a way to extract the intersections widths from the resulting object.

It this the wrong approach ?

Note: Using intersect(gr1, gr2) in this case yields a single contiguous overlap range spanning the whole chromosome, but I want each overlap to be returned independently, just like findOverlaps does.

$\endgroup$
2
  • 1
    $\begingroup$ Please add a representative output example. $\endgroup$
    – user3051
    Jan 21, 2021 at 14:26
  • $\begingroup$ Thanks, I just added it $\endgroup$
    – cmdoret
    Jan 21, 2021 at 14:36

1 Answer 1

1
$\begingroup$

Whilst intersect is sort of doing what you want, it merges adjacent overlaps together which is not what you want I suppose. There doesn't seem to be dedicated Grange function to do this out of the box, but you could go through all overlapping pairs resulting from findOverlap and calculate it's widths using pintersect:

punion, pintersect, psetdiff and pgap are generic functions that compute the element-wise (aka "parallel") union, intersection, (asymmetric!) difference and gap between each element in x and its corresponding element in y.

"pintersect"(x, y, resolve.empty=c("none", "max.start", "start.x"), ...) # x, y IRanges objects.

hits <- findOverlaps(gr1, gr2)

width(pintersect(gr1[queryHits(hits)],
                 gr2[subjectHits(hits)]))
[1]  2  6  4  6 10
$\endgroup$
5
  • 1
    $\begingroup$ If I understand correctly, pintersect only computes the intersection between pairs of ranges at the same position. (e.g. 2nd vs 2nd). This implies you have the same number of ranges in gr1 and gr2. I am interested in the all-vs-all intersection. I will update my example to make that clearer, sorry. $\endgroup$
    – cmdoret
    Jan 21, 2021 at 13:42
  • $\begingroup$ Thanks that answers the question :) The only issue is that the R loop is super slow. In my case I have 468,000 intersections and it's still running after 10mins. I have found an alternative which runs in 0.1 sec, but it is imperfect as I have to truncate intervals by 1bp (slightly wrong results): mywidths <- width(intersect(narrow(gr1, start=2), narrow(gr2, start=2))); comp_df <- data.frame(feat1 = ol@from, feat2 = ol@to, size = mywidths) $\endgroup$
    – cmdoret
    Jan 21, 2021 at 15:45
  • $\begingroup$ That a good point. I've added another possibly quicker solution, hope that works! $\endgroup$ Jan 21, 2021 at 17:48
  • $\begingroup$ That last solution is a great idea, thanks ! $\endgroup$
    – cmdoret
    Jan 21, 2021 at 17:51
  • $\begingroup$ If it works for, I'd consolidate the answer since the first iterations are superflous then :) $\endgroup$ Jan 21, 2021 at 17:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.