This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How do I subset a GRanges on chromosome, region and strand?
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
• 21,994 views • link
updated by Bastien Hervé Bastien Hervé • written by endrebak852 • endrebak852

Read about findOverlaps?

• 1 views • link
written by zx8754 zx8754

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

• 1 views • link
written by ATpoint ATpoint

3 answers

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)
• 1 views • link
written by zx8754 zx8754

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 views • link
written by bernatgel bernatgel

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

• 1 views • link
written by zx8754 zx8754 
gr=GRanges(seqnames=c("chr1","chr2","chr2"), ranges=IRanges(start=c(50,150,200),end=c(100,200,300)), strand=c("+","-","-"))
gr_filter=GRanges(seqnames="chr2", ranges=IRanges(start=200,end=300), strand="-")

subsetByOverlaps(gr, gr_filter)

#GRanges object with 2 ranges and 0 metadata columns:
#seqnames    ranges strand
#       <Rle> <IRanges>  <Rle>
#[1]     chr2   150-200      -
#[2]     chr2   200-300      -
#-------
#seqinfo: 2 sequences from an unspecified genome; no seqlengths
• 0 views • link
written by Bastien Hervé Bastien Hervé

How is this different from my post?

• 1 views • link
written by zx8754 zx8754

Ah ! I wrote my answer while you were typing yours I guess :) I will delete it

• 1 views • link
written by Bastien Hervé Bastien Hervé

Log in to answer this question.

More posts like this