understanding pysam pileup()

Question

[GouthamAtla]

I am trying to understand how the pysam pileup() method works. Here is my example code:

for pileupcolumn in samfile.pileup("ConsensusB", 1, 2):
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            print pileupread.query_position,pileupread.alignment.query_name,pileupread.alignment.query_sequence[pileupread.query_position]

In the above example, I would expect it to give me only reads overlapping at base 1 of reference genome. But I see that it still iterates through each and every base of the reference sequence. How do I get all the bases overlapping a particular position (let's say a SNP from VCF)

Answer 1

[GouthamAtla]

Thanks. I could use https://github.com/alimanfoo/pysamstats for reference.

Great examples here

[Ram]

Why delete the question?

[Istvan Albert]

Indeed, the right etiquette is to answer your own question, accept the answer.