How do I subset a GRanges on chromosome, region and strand?
3
2
Entering edit mode
2.4 years ago
endrebak852 ▴ 110
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.

GenomicRanges R bioconductor • 7.3k views
0
Entering edit mode

Read about findOverlaps?

1
Entering edit mode

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

7
Entering edit mode
2.4 years ago
zx8754 10k

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)

3
Entering edit mode
2.4 years ago
bernatgel ★ 2.9k

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.

1
Entering edit mode

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