Question: genomicranges outside sequence bounds
0
gravatar for reviewer3
5.2 years ago by
reviewer320
United States
reviewer320 wrote:

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!

 

ADD COMMENTlink modified 5.2 years ago by seidel6.8k • written 5.2 years ago by reviewer320
3
gravatar for seidel
5.2 years ago by
seidel6.8k
United States
seidel6.8k wrote:

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.

 

 

 

 

ADD COMMENTlink modified 5.2 years ago • written 5.2 years ago by seidel6.8k

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!

 

ADD REPLYlink modified 5.2 years ago • written 5.2 years ago by reviewer320

And,

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

gives me the same exact answer. Thank you!

ADD REPLYlink written 5.2 years ago by reviewer320
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: 2115 users visited in the last hour