Marco, Now I understand the problem better. I see you are right that 'gaps' computes right and left flanking regions with respect to the ranges provided in start and end Careful reading of the docs reveals that gaps is working as advertised. It is not a bug. ? gaps,GenomicRanges states " 'gaps' returns an object of the same type as 'x' containing complemented ranges for each distinct (seqname, strand) pairing." where strand of '*' is distinct from '+' and '-' It alos refers us to gaps,IRanges for details about treatment of optional start,end parameters, where we learn: " 'gaps' returns the "normal" Ranges object representing the set of integers that remain after the set of integers represented by 'x' has been removed from the interval specified by the 'start' and 'end' arguments." and "'start=NA' and 'end=NA' are interpreted as 'start=1' and 'end=length(subject(x))'" (it is stated clearly for Views, less so for IRanges, but, we get the message). So, you just need to remove them, like this: intergenic<- intergenic[strand(intergenic)!='*'] Or use the approach I provided using setdiff. Cheers, Malcolm >-----Original Message----- >From: Blanchette, Marco >Sent: Thursday, September 05, 2013 8:08 PM >To: Cook, Malcolm >Cc: BioConductor >Subject: Re: bug in GenomicRanges gaps? > >though... sorry about the genes.noStrand... was working on a different example before posting... bad pasting... > >And yes, I was trying to find the gaps between genic regions regardless of there strand. I was doing the same operations as yours until >I fumble onto gaps() which seemed to encapsulate nicely in one function the concept of finding gaps between ranges (right??) and if >you look at the operation returned, it uses the seqinfo to compute the right and left gaps in relation to the entire chromosme. >However, it returns extra ranges that conceptually are not gaps (ie span the full chromosome)... > >But thanks for the tidbits Malcolm. >-- >Marco Blanchette, Ph.D. >Assistant Investigator >Stowers Institute for Medical Research >1000 East 50th St. > >Kansas City, MO 64110 > >Tel: 816-926-4071 >Cell: 816-726-8419 >Fax: 816-926-2018 > >On Sep 5, 2013, at 7:53 PM, "Cook, Malcolm" <mec at="" stowers.org=""> > wrote: > >> Marco, >> >> genes.noStrand is unset when you first reference it, so it is hard to know what you are trying to do. >> >> Does your definition of intergenic want to ignore strand (i.e. can the same locus be genic on one strand and integenic on another?) >> >> I'm guessing the problem you are witnessing is that the strand _is_ not '*' when you thought you had set it to '*'. Check this! >> >> Also, does your definition of intergenic want to include non-genic segments that are the very beginning/end of a chromosome (and >therefore only have a gene on one-side and thus arguably are not truly inter-genic). If so, your approach is not going to include these >chromosome-terminal regions. In which case, you might try: >> >> genic.gr<-unlist(transcriptsBy(tr.db,'gene'),use.names=FALSE) >> strandgenic.gr)<-'*' >> genic.gr<-reducegenic.gr) >> chrom.gr<-as(seqinfo(tr.db),'GRanges') >> intergenic.gr<-setdiffchrom.gr,genic.gr) >> >> stopifnot(0==length(findOverlapsintergenic.gr,ignoreSelf=TRUE))) >> >> ~ malcolm_cook at stowers.org >> >> ________________________________________ >> From: bioconductor-bounces at r-project.org [bioconductor-bounces at r-project.org] on behalf of Blanchette, Marco >> Sent: Thursday, September 05, 2013 7:14 PM >> To: BioConductor >> Subject: [BioC] bug in GenomicRanges gaps? >> >> Hi, >> >> Is this a bug of normal behavior... (I suspect it's a bug). I am trying to compute the intergenic regions using GenomicFeatures txdb >and GenomicRanges gaps() and the GRanges return contains entry spanning all chromosomes on both strand... On a mock GRanges, I >don't get this weird behavior (See below). Any suggestions? >> >> thanks >> >>> library(GenomicFeatures) >> ### Downloading the txdb from UCSC >>> txdb <- makeTranscriptDbFromUCSC(genome="dm3", >> tablename="ensGene") >> >> ### Getting the transcripts for each genes >>> allTx <- transcriptsBy(txdb, by='gene') >> >> ### Geting the GRanges for each genes >>> genes <- unlist(range(allTx)) >> >> ### Removing the strand to perform the reduce operation >>> strand(genes) <- '*' >> >> ### Now, let's find the gaps >>> genic <- reduce(genes.noStrand) >> >> ### Because of teh reduce, there should be no overlaping genic regions >>> length(findOverlaps(genic,ignoreSelf=TRUE)) >> [1] 0 >> >> ### Finding the intergenic regions >>> intergenic <- gaps(genic) >> >> ### The difference of non-overlaping ranges should return non- overlaping gaps, but... >>> length(findOverlaps(intergenic,ignoreSelf=TRUE)) >> [1] 46300 >> >> ### What's funky is that I am getting gaps spanning the full lenght of each chromosomes >> ### Why??? >>> chrs.width <- seqlengths(seqinfo(intergenic)) >>> intergenic[width(intergenic) == gene2chrsLength] >> GRanges with 30 ranges and 0 metadata columns: >> seqnames ranges strand >> <rle> <iranges> <rle> >> [1] chr2L [1, 23011544] + >> [2] chr2L [1, 23011544] - >> [3] chr2R [1, 21146708] + >> [4] chr2R [1, 21146708] - >> [5] chr3L [1, 24543557] + >> ... ... ... ... >> [26] chrXHet [1, 204112] - >> [27] chrYHet [1, 347038] + >> [28] chrYHet [1, 347038] - >> [29] chrUextra [1, 29004656] + >> [30] chrUextra [1, 29004656] - >> --- >> seqlengths: >> chr2L chr2R chr3L chr3R ... chrXHet chrYHet chrUextra >> 23011544 21146708 24543557 27905053 ... 204112 347038 29004656 >> >> >> ### On Fake data >>> genes <- GRanges(seqnames=c(rep("2L",3),rep("3L",3)), >> ranges=IRanges(start=rep(c(20,200,350),2),end=rep(c(80,275,450),2))) >>> length(findOverlaps(intergenic,ignoreSelf=TRUE)) >> [1] 0 >>> intergenic >> GRanges with 6 ranges and 0 metadata columns: >> seqnames ranges strand >> <rle> <iranges> <rle> >> [1] 2L [ 1, 19] * >> [2] 2L [ 81, 199] * >> [3] 2L [276, 349] * >> [4] 3L [ 1, 19] * >> [5] 3L [ 81, 199] * >> [6] 3L [276, 349] * >> --- >> seqlengths: >> 2L 3L >> NA NA >> >> -- >> Marco Blanchette, Ph.D. >> Assistant Investigator >> Stowers Institute for Medical Research >> 1000 East 50th St. >> >> Kansas City, MO 64110 >> >> Tel: 816-926-4071 >> Cell: 816-726-8419 >> Fax: 816-926-2018 >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor >>