I am trying to write a python script that finds bases at specific sites in each read (aligning to the site). I am using pysam pileup, but cannot figure out how to use it properly. Could someone help me out with this ?
the input data is a list of coordinates and a samfile = pysam.AlignmentFile("bam", "rb")
I just want to iterate through each read in a bam aligning to a single base site, and extract bases at each of the sites.
I tried doing this, but it's outputting each base multiple times.
for pileupcolumn in samfile.pileup('1', 14693, 14694): for pileupread in pileupcolumn.pileups: print pileupread.alignment.query_sequence[pileupread.query_position]