Question: How do I subset a GRanges on chromosome, region and strand?
1
gravatar for endrebak852
12 months ago by
endrebak852100
github.com/endrebak
endrebak852100 wrote:
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 bioconductor R • 2.2k views
ADD COMMENTlink modified 12 months ago by Bastien HervĂ©4.5k • written 12 months ago by endrebak852100

Read about findOverlaps?

ADD REPLYlink written 12 months ago by zx87548.7k
1

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

ADD REPLYlink written 12 months ago by ATpoint26k
4
gravatar for zx8754
12 months ago by
zx87548.7k
London
zx87548.7k wrote:

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)
ADD COMMENTlink written 12 months ago by zx87548.7k
2
gravatar for bernatgel
12 months ago by
bernatgel2.2k
Barcelona, Spain
bernatgel2.2k wrote:

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.

ADD COMMENTlink modified 12 months ago • written 12 months ago by bernatgel2.2k
1

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

ADD REPLYlink written 12 months ago by zx87548.7k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1170 users visited in the last hour