[BioC] Rsamtools::countBam returning in correct counts

[Hubert Rehrauer]

Hi Martin

thanks for the clarification, yes it works. Sorry my fault.

And there are unmapped reads in the file. The bam holds paired-end 
alignments from the LifeScope software. If only one of the reads in a 
pair gets mapped Lifescope still reports both, and obviously one of them 
has the unmapped flag.

regards,
hubert


On 05.07.12 19:31, Martin Morgan wrote:
> 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
>
>