Intersection of two genomic ranges to keep metadata

Question

I am trying to find intersection of two genomic ranges (gr1 and gr2) and keep metadata from one of them

gr1
chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

gr2
chrI [50, 150] + |

intersect(gr1, gr2) will remove metadata, where subsetByOverlaps(gr1, gr2) return:

chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

where metadata is saved, but intersection is bigger than real (the same as in initial gr1 start/end)

How can I get an intersection range and keep metadata? to get the range:

chrI [50, 100] * | 0.1
chrI [101, 150] * | 0.2

Devon Ryan · Accepted Answer · 2017-06-22T11:32:27.683

Presuming you don't have cases of multiple ovelaps

> library(GenomicRanges)
> gr1 = GRanges(c("chrI", "chrI"), IRanges(start=c(1, 101), width=c(100, 100)), mcols=data.frame(foo=c(0.1, 0.2)))
> gr2 = GRanges(c("chrI"), IRanges(start=c(50), width=c(101)))
> o = findOverlaps(gr1, gr2)
> grl1 = split(gr1[queryHits(o)], 1:length(o)) # You can't mendoapply on a GRanges object
> grl2 = split(gr2[subjectHits(o)], 1:length(o))
> foo = function(x, y) {
+     rv = x
+     start(rv) = max(start(x), start(y))
+     end(rv) = min(end(x), end(y))
+     return(rv)
+ }
> unlist(mendoapply(foo, grl1, y=grl2))

The output is:

GRanges object with 2 ranges and 1 metadata column:
    seqnames     ranges strand | mcols.foo
       <Rle>  <IRanges>  <Rle> | <numeric>
  1     chrI [ 50, 100]      * |       0.1
  2     chrI [101, 150]      * |       0.2
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

It'd be more memory efficient to iterate over o, but this gives you the idea.

In cases of multiple overlaps, one could use intersect(), then findOverlaps() with that and then continue on with the procedure demonstrated above (or iterate over the overlaps).

Intersection of two genomic ranges to keep metadata

1 Answers1