I am trying to find intersection of two genomic ranges (gr1 and gr2) and keep metadata from one of themgr1chrI [1, 100] * | 0.1chrI [101, 200] * | 0.2gr2chrI [50, 150] + | intersect(gr1, gr2) will remove metadata, wheresubsetByOverlaps(gr1, gr2) return:chrI [1, 100] * | 0.1chrI [101, 200] * | 0.2where 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.1chrI [101, 150] * | 0.2
Intersection of two genomic ranges to keep metadata
r;bioconductor;granges
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 seqlengthsIt'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).