9
$\begingroup$

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
$\endgroup$
4
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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