How do I subset a GRanges on chromosome, region and strand?

Question

[endrebak852]

library(GenomicRanges)

gr=GRanges(seqnames=c("chr1","chr2","chr2"),
           ranges=IRanges(start=c(50,150,200),end=c(100,200,300)),
           strand=c("+","-","-")
)

I want to get all intervals between 200-300 on chromosome 2, the minus strand. How do I do that?

gr[seqnames(gr) == "chr2" & start(gr) > 200 & end(gr) < 300 & strand(gr) == "-"]

does not work as it does not find the overlaps with 200-300, but rather the intervals strictly contained in that range.

[zx8754]

Read about findOverlaps?

[ATpoint]

Construct a second GRanges object with the target ranges and then use either findOverlaps as zx8754 suggests, or subsetByOverlaps.

Answer 1

[zx8754]

Make a query range, then subsetByOverlaps (sorry in comments by mistake mentioned findOverlaps):

q=GRanges(seqnames="chr2",
          ranges=IRanges(start = 200, end = 300),
          strand="-")

subsetByOverlaps(gr, q)

Answer 2

[bernatgel]

If you want to use a simple selector like the one you propose, you should flip the tests for start and end:

gr[seqnames(gr) == "chr2" & start(gr) < 300 & end(gr) > 200 & strand(gr) == "-"]

and that should work.

Using subsetByOverlaps could be more readable and is your best option if you have more than one region.

[zx8754]

Good one, didn't spot there was a mistake in logical comparisons. Even if subsetByOverlaps is better for general case, this is the answer to the "problem". (I will move this to an answer).