genomicranges outside sequence bounds

Question

[reviewer3]

Hi,

I'm new to genomic ranges and getting a bit frustrated with it. I'm reading the reference manual and tutorials, but haven't found an answer to the following question:

I've constructed my gr object and assigned sequence lengths to it. I must have made an error, as I'm getting the warning:

In valid.GenomicRanges.seqinfo(x) :
  'ranges' contains values outside of sequence bounds. See ?trim to subset ranges.

When I use trim(gr), I simply get back a new gr that is trimmed. But it doesn't show me offending ranges. How can I do that so I can learn where I went wrong in constructing gr in the first place?

Thank you!

Answer 1

[seidel]

You didn't mention how you are constructing your GRanges object, but if you've already constructed it you can check to see if any of your GRanges elements are longer than the sequence lengths you are trying to set as follows (assuming you have a BSgenome object with the chromosome lengths, the example below uses hg19 (not required, any object/vector with sequence lengths would do)):

library(BSgenome.Hsapiens.UCSC.hg19)
which(end(mygr) > seqlengths(BSgenome.Hsapiens.UCSC.hg19)[as.character(seqnames(mygr))])

This returns the ends of each GRanges element, and checks it against the length of the corresponding chromosome.

Also, if you adjust the lengths of one gr object and want to check it against the untrimmed on to see which elements differ, you might try making a copy before trimming, and then comparing the trimmed object to the original

which(mygr.original != mygr.trimmed)

assuming they are in the same order.

[reviewer3]

Thank you, seidel @seidel

which(mygr.original != mygr.trimmed)

gave me exactly what I needed. I got the sense that these gr containers aren't very intuitive to use, but this is obviously a great example where they are. Thank you!

[reviewer3]

And,

which(end(mygr) > seqlengths(BSgenome.Hsapiens.UCSC.hg19)[as.character(seqnames(mygr))])

gives me the same exact answer. Thank you!