pysam pileup wierdness

[mbreese]

I'm trying to track down a strange behavior we recently noticed using pysam's pileup function. We are trying to find SNPs at specific positions. Depending on how we call the pileup function, we get different results.

If we do the call like this:

for p in bamfile.pileup():
 if p.pos == MYPOS: #assume only one reference/chromosome in bam file
  for read in p:
   print read.alignment.qname

we get one set of reads.

If we call it like this:

for p in bamfile.pileup(chrom,pos,pos+1):
 if p.pos == MYPOS:
  for read in p:
   print read.alignment.qname

We get a completely different set of reads (some overlap, but not much). From my reading of the source, I can't find any reason why the results should be different. Or even which set of reads is correct. Has anyone else ever seen this?

Thanks

pysam: 0.4.1
samtools: 0.1.12a (pysam embedded)

[mbreese]

We found the problem... and it's a known issue. In case anyone else has this problem: samtools pileup can only handle ~8000 reads at a position. We were getting different results depending on how the function was called, but it's probably due to a limited internal (samtools) buffer.

http://code.google.com/p/pysam/issues/detail?id=63