Pysam pileup working unpredictably

[polarise]

Hi,

I'm trying to count allelic frequencies at SNPs and Pysam seems to be the best solution to do that. However, when I make a call to get the number of reads that occur at a single site I get many more sites than expected.

Example:
len(list(f.pileup(reference='chr1',start=1756880,end=1756881)))

gives a result of 135 as opposed to the expected value of 1

and

len(list(f.pileup(reference='chr1',start=1756870,end=1756881)))

gives a result of 146 as opposed to 10.

Can anyone please explain what is going on?

Regards,

P. Korir

[ffinkernagel]

I suspect that pileup givess you pileups for every base covered by a read that was within the query range.

edit: You might want to print the pos of each PileupColumn. for p in f.pileup(...): print p.pos, p.n ...

[polarise]

Thanks for your reply ffinkernagel.

Well, I'm surprised that when I give it a query range of 1 it gives more than one result. Anyway, I have used the for loop like in your suggestion though the solution is not as clean as I wanted it to be.

f.pileup() returns an iterator. I'm not proficient in the use of iterators.

P. Korir

[ffinkernagel]

you simply do
for something in my_iterator:
do_something_with(something)

no difference to a list there.

The for loop should tell you if my hunch is correct.

[polarise]

OK. I'll give it a shot.

Thank you.