!! Please pardon the line-breaks MS Outlook auto-inserted into my code in a previous version of this post.
Martin, Hmmm, Thanks .... I think I'm getting it.... Following your lead, I can directly re-order cnt in the OriginalOrder as: > cnt[sort(unlist(split(values(which2), > seqnames(which2)))$OriginalOrder,index.return=TRUE)$ix,] space start end width file records nucleotides 2 seq2 100 1000 901 ex1.bam 1169 41235 1 seq1 1000 2000 1001 ex1.bam 612 21549 3 seq2 1000 2000 1001 ex1.bam 642 22640 ... which can usefully(?) be abstracted as countBamWhich <- function (file,which,index=file,...) { ### wrapper for countBam that reorders its results to aggree with ### 'which', a required ScanBamParam. ... params are additional ### ScanBamParam options. ### ### ASSUMES: internal implementation detail of ScanBamParam. c.f. ### http://permalink.gmane.org/gmane.science.biology.informatics.conductor/34208) param <- ScanBamParam(which=which,...) values(which)[['OriginalOrder']] <- 1:length(which) CBW = countBam(file,index,param=ScanBamParam(which=which)) CBW[sort(unlist(split(values(which), seqnames(which)))$OriginalOrder,index.return=TRUE)$ix,] } allowing me to write > countBamWhich(fl, which2) space start end width file records nucleotides 2 seq2 100 1000 901 ex1.bam 1169 41235 1 seq1 1000 2000 1001 ex1.bam 612 21549 3 seq2 1000 2000 1001 ex1.bam 642 22640 All in favor? ~ Malcolm Malcolm Cook Stowers Institute for Medical Research - Bioinformatics Kansas City, Missouri USA > -----Original Message----- > From: Martin Morgan [mailto:mtmor...@fhcrc.org] > Sent: Thursday, March 31, 2011 11:41 AM > To: Cook, Malcolm > Cc: 'myo...@wehi.edu.au'; 'bioconduc...@r-project.org'; > 'Bioc-sig-sequencing@r-project.org' > Subject: Re: [BioC] Re(surrecting): [Bioc-sig-seq] Rsamtools > countBam labeling > > On 03/31/2011 08:32 AM, Cook, Malcolm wrote: > > Matt et. al., > > > > I wonder if a satisfactory resolution to the issue of "the order > > changes between the GRanges object and the countBam data.frame > > > > > http://www.mail-archive.com/bioc-sig-sequencing@r-project.org/ > msg01144.html > > > > I am presented with the same issue and poised to tackle it > but wondr > > if a generic solution emerged from you inquiries& efforts. > > Hi Malcolm -- > > For a reproducible example, > > library(Rsamtools) > example(countBam) > which1 <- as(which, "GRanges") > ## which2 might be where your data actually starts > which2 <- which1[c(2,1,3)] > values(which2)[["OriginalOrder"]] <- 1:3 > param <- ScanBamParam(which=which2) > cnt <- countBam(fl, param=param) > > What happens is that ScanBamParam converts its argument to an > IRangesList, using split(ranges(which2), seqnames(which2)). So do the > same for the values and unlist > > cntVals <- unlist(split(values(which2), seqnames(which2))) > > then cbind coerced values > > cbind(cnt, as.data.frame(cntVals)) > > with > > > which2 > GRanges with 3 ranges and 1 elementMetadata value > seqnames ranges strand | OriginalOrder > <Rle> <IRanges> <Rle> | <integer> > [1] seq2 [ 100, 1000] * | 1 > [2] seq1 [1000, 2000] * | 2 > [3] seq2 [1000, 2000] * | 3 > > seqlengths > seq1 seq2 > NA NA > > cbind(cnt, as.data.frame(cntVals)) > space start end width file records nucleotides OriginalOrder > 1 seq1 1000 2000 1001 ex1.bam 612 21549 2 > 2 seq2 100 1000 901 ex1.bam 1169 41235 1 > 3 seq2 1000 2000 1001 ex1.bam 642 22640 3 > > Martin > > > > > Thanks, > > > > Malcolm Cook Stowers Institute for Medical Research - > > Bioinformatics Kansas City, Missouri USA > > > > > > _______________________________________________ Bioconductor mailing > > list bioconduc...@r-project.org > > https://stat.ethz.ch/mailman/listinfo/bioconductor Search the > > archives: > > http://news.gmane.org/gmane.science.biology.informatics.conductor > > > -- > Computational Biology > Fred Hutchinson Cancer Research Center > 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 > > Location: M1-B861 > Telephone: 206 667-2793 > _______________________________________________ Bioc-sig-sequencing mailing list Bioc-sig-sequencing@r-project.org https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing