[Bioc-devel] odd param behavior with scanBam from Rsamtools 1.14.1

Jason Byars jbyars at unm.edu
Thu Oct 31 21:06:54 CET 2013


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



More information about the Bioc-devel mailing list