Hi all,

I am trying to use the gaps function from the GenomicRanges package, and am experienced unexpected behavior. For example, if we have the following GRanges object:

gr <- GRanges(seqnames = "chr1", strand = c("*", "*"), ranges = IRanges(start = c(1, 5), width = 3))

seqlengths(gr) <- c("chr1" = 10)

I would expect

​gaps(gr)

To return a GRanges object equal to 

GRanges(seqnames = "chr1", strand = c("*", "*"), ranges = IRanges(start = c(4, 8), end = c(4, 10)))

However, it returns a GRanges object equal to

GRanges(seqnames = "chr1", strand = c("+", "-", "*", "*"), ranges = IRanges(start = c(1, 1, 4, 8), end = c(10, 10, 4, 10)))

It doesn't seem like you would want gaps to return an GRanges object that contains ranges that are the width of the whole sequence. Is there a way to force it to return the output I was expecting?

Jacob

The manual explains that this is expected. In this case, you can just subset to the ranges in the "*" strand. In general, if I want to do this, I stay far away from the unintuitive gaps() and do:

setdiff(as(seqinfo(gr), "GRanges"), gr)

This works perfectly, thanks so much.