I have a problem working with BAM files and pysam whenever I try to examine bases at a certain position within a sequencing read. What I want to do is look at all the reads covering a region (which works fine with pysam, using fetch()) and compare the base of each read at a position, thus getting e.g. all read information about the 16,078,956. base pair in the genome. Something like Read.query_alignment_sequence gives me the first base in the aligned read, but getting the element at the position I want is a bit complicated - it now works by comparing coordinates.
Whenever a read contains an indel, comparing coordinates doesn't work anymore because "query_alignment_sequence" then is "shifted" and gets the positions wrong.
Is there any efficient way to work this out? If somebody has an idea, please let me know.
Thanks in advance!