reads in pysam

Question

[yzhuan0621]

Hey, I am a newbie in pysam, and now I have a question.

I use pysam to find the reads which overlap a given position, and then I want to get the nucleotides at this position in each read. Followed my code:

position is an given location(0 based)

reads = [read for read in samfile.fetch(chr_name, position,position+1)]
    for read in reads:
         neu1 = read.seq[read.get_reference_positions().index(position)]
         neu2 = read.seq[position-read.reference_start]
         neu3 = read.query[position-read.reference_start-read.qstart]

I want to know which is right, neu1, nue2, nue3, or others? Thanks!

Answer 1

[Devon Ryan]

None of those methods are guaranteed to give you the correct bases. Use the pileup() function, which will provide ready access to the information you want.