[Bioc-devel] odd param behavior with scanBam from Rsamtools 1.14.1
Ok, that makes more sense. Thanks! Jason -----Original Message----- From: Martin Morgan [mailto:mtmorgan at fhcrc.org] Sent: Thursday, October 31, 2013 2:03 PM To: Jason Byars; bioc-devel at r-project.org Subject: Re: [Bioc-devel] odd param behavior with scanBam from Rsamtools 1.14.1
On 10/31/2013 12:54 PM, Jason Byars wrote:
Using any permutation of ScanBamParam() options seems to work fine with
countBam, filterBam, etc. But, when I use a ScanBamParam object without anything specified for "what" with scanBam, I get a length 0 list back. If all I add is what=ScanBamWhat(), the call appears to work fine. I have tried this on both Windows and Linux, and with a clean install of R and Bioconductor 2.13. Is this behavior correct? Example script follows: yes, its the way it was designed, however poorly -- you don't ask for any fields, so you don't get any fields. countBam and filterBam are about visiting records, and records exist even when none of their fields are consulted. It's like > data.frame(x=1:10)[,FALSE, drop=FALSE] data frame with 0 columns and 10 rows where there are 10 rows (records), but no columns (fields). Martin
library(Rsamtools)
# param seems to work find on filterBam, but I'm getting no results on
scanBam on Linux
sampleFile <- "F:/nessstuff/rnaseq/mybunaligned/P61_3_MBC24-6-1-M3_mybunaligned.bam"
# just picked MYB because we work on it. It's not important to show the problem.
mybStart <- 135502453
mybEnd <- 135540311
mybSize <- 37859
padding <- 10000
flag <- scanBamFlag(isUnmappedQuery=F) which <-
GRanges(seqname=Rle(c("chr6"), c(1)), ranges=IRanges(mybStart-padding,
mybEnd+padding)) # what <- c("rname", "strand", "pos", "qwidth",
"seq") param <- ScanBamParam(flag=flag,which=which)
countBam(sampleFile) #813k records
countBam(sampleFile,param=param) #2397 records
# 813k records loaded
bam <- scanBam(sampleFile)[[1]]
summary(bam)
# 0 length list returned
bam <- scanBam(sampleFile,param=param)[[1]]
summary(bam)
# specify anything for what and it works as expected #2397 records
imported param <- ScanBamParam(flag=flag,which=which,
what=scanBamWhat()) bam <- scanBam(sampleFile,param=param)[[1]]
summary(bam)
sessionInfo()
R version 3.0.2 (2013-09-25) Platform: x86_64-w64-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] Rsamtools_1.14.1 Biostrings_2.30.0 GenomicRanges_1.14.3 XVector_0.2.0 IRanges_1.20.4 [6] BiocGenerics_0.8.0 BiocInstaller_1.12.0 loaded via a namespace (and not attached): [1] bitops_1.0-6 stats4_3.0.2 tools_3.0.2 zlibbioc_1.8.0 [[alternative HTML version deleted]]
_______________________________________________ Bioc-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/bioc-devel
-- Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone: (206) 667-2793