This is a test version of Biostars. For the public version, visit https://www.biostars.org.
reads in pysam

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!

genome snp pysam alignment

1 answer

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.

Log in to answer this question.