Rsamtools::countBam returning in correct counts
1
0
Entering edit mode
@hubert-rehrauer-3823
Last seen 9.6 years ago
Dear developers In the following case countBam does not give the expected result: > param = ScanBamParam() > bamFlag(param) = scanBamFlag(isUnmappedQuery=TRUE) > 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? regards, hubert [[alternative HTML version deleted]]
• 1.7k views
ADD COMMENT
0
Entering edit mode
@martin-morgan-1513
Last seen 6 weeks ago
United States
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
ADD COMMENT
0
Entering edit mode
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 > >
ADD REPLY

Login before adding your answer.

Traffic: 786 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6