.fetch and .count different from samtools view

[GoogleCodeExporter]

What steps will reproduce the problem?

1. create a genome with sequence "chr1:AAACCCTTTGGG", save to genome.fasta
1.1. bowtie-build genome.fasta genome
2. bowtie genome -c AAAGG -a -v 0 --sam > test.sam
3. samtools view -bS test.sam > test.bam
4. samtools sort test.bam test-sorted
5. samtools index test-sorted.bam

With this:

samtools view -c test-sorted.bam chr1:1-5

you get the aligned read (aligned to - strand, position 5), so the count is 1.

But if you use pysam and do:

bam = pysam.Samfile("test-sorted.bam", "rb")
bam.count("chr1", 0, 4)

the result is 0.

Is this a bug or am i missing something?

Gregor

Original issue reported on code.google.com by gregor....@gmail.com on 15 Jun 2012 at 2:10

[GoogleCodeExporter]

Python (and pysam) is right exclusive - so a read starting at 4 (or 5 in 
SAM-coordinates) would not be found by your query, I believe.

Original comment by finkerna...@mathematik.uni-marburg.de on 20 Jun 2012 at 9:49

[GoogleCodeExporter]

But samtools finds it and counts it (so it's not right exclusive it seems?), 
Tnx.

Original comment by gregor....@gmail.com on 20 Jun 2012 at 9:51

[GoogleCodeExporter]

Indeed - samtools differs from pysam in its ranges. Pysam follows python 
convention.

Best wishes,
Andreas

Original comment by andreas....@gmail.com on 7 Jul 2012 at 9:37

• Changed state: WontFix