Pysam Snp Calling

Question

[Biomed]

I originally asked this question on the pysam mailing list but my experience has been that the developers are not supporting pysam through that mailing list (none of my questions were answered in the same day and some were not answered at all, and the documentation is not very detailed). I still like the idea of using pysam and if biostar proves to be a better platform for support I plan to continue using it.

Here is my question: If I wanted to call snps over a region of 100KB in chr16 can I use IteratorSNPCalls and specify a startposition and iterate from there or do I need to iterate over the whole file from start until the position I am interested? Thanks

my example per below answer: SNPCaller example:

pileup_iter = samfile.pileup (stepper = 'samtools' , fastafile=fastafile)
    snpcaller= pysam.SNPCaller(pileup_iter)
    sc= snpcaller.call('chr16',73310450)

IteratorSNPCalls example:

pileup_iter = samfile.pileup ('chr16',73310450,73310460, stepper = 'samtools' , fastafile=fastafile)
snpcall_iter= pysam.IteratorSNPCalls(pileup_iter)
for call in snpcall_iter:
    print str(call)

[Sean Davis]

Thanks for the code addition.

[Biomed]

No problem, thanks for the right direction.

Answer 1

[Sean Davis]

From the documentation for IteratorSNPCalls, the argument to IteratorSNPCalls is a pileup iterator. You can create a pileup iterator over your region of interest and then use that as input to IteratorSNPCalls. Refer to the pileup engine section of the documentation for details.

[Biomed]

Of course, the region is selected through the pileup iterator. thanks a lot.

[Biomed]

Yes I was able to get an interval based pileup_iterator object and run it that way. One needs to remember that the position information is supplied through SNPCaller.call() whereas with Iterator.SNPCalls(pileup) you don't need the interval in the snpcaller because it is already coming from the pileup engine.