[BioC] Rsamtools::countBam returning in correct counts

[Martin Morgan]

Hi Hubert --

On 07/05/2012 11:33 AM, Hubert Rehrauer wrote:
> Dear developers
>
> In the following case countBam does not give the expected result:
>
>>   param = ScanBamParam()

this creates bamWhat(param) == character(), which is to say 'don't read 
anything in', and explains why scanBam returns a 0-length named list.

>>   bamFlag(param) = scanBamFlag(isUnmappedQuery=TRUE)

It's also possible to create this as

   ScanBamParam(flag=scanBamFlag(isUnmappedQuery=TRUE))

though probably you want

   ScanBamParam(flag=scanBamFlag(isUnmappedQuery=TRUE),
                what="qname")

>>   cb = countBam(bamFile, param=param)
>>   cb
>     space start end width     file records nucleotides
> 1    NA    NA  NA    NA 4541.bam 3157955   133205755
>>   bamReads = scanBam(bamFile, param=param)[[1]]$qname
>>   scanBam(bamFile, param=param)[[1]]
> named list()
>
>
> There are no unmapped reads in my bamFile and so it should return 0, or am I wrong?

So maybe there are unmapped reads (if setting what="qname" returns 
something)?

   example(scanBamFlag)

and then exploring the bam file pointed to by 'fl' seem to suggest that 
the flag is being interpreted correctly. If you're not convince, it 
would help to have sessionInfo() and bamFile, and if possible a small 
sample (use ?filterBam?) of your bam file?

Martin

>
> regards,
> hubert
>
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> 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


-- 
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