Find all the bases for given reference position

Question

Hi im looking for the aligned bases in the reads for a given reference position. im using the following script from the pysam documentataion. I adjusted it to find the specified position. in this case 24793.

import pysam

samfile = pysam.AlignmentFile("generated_alignment_sorted.bam", "rb" )
for pileupcolumn in samfile.pileup("chr05_modified.copy0", 10, 52000000):
    if pileupcolumn.pos == 24793:
        print ("\ncoverage at base %s = %s" %
               (pileupcolumn.pos, pileupcolumn.n))
        for pileupread in pileupcolumn.pileups:
                print ('\tbase in read %s = %s' %
                      (pileupread.alignment.query_name.split(';')[0],
                       pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

This outputs:

coverage at base 24793 = 6
        base in read m419941/6207/CCS Read=419941 = C

For all the position I tried it doesnt print out the same number of bases as it says that the coverage is. shouldnt this print out 6 times: base in read m419941/6207/CCS Read=419941 = C for different reads?

Why is this?

Answer 1

pileupcolumn.n shows the total number of the reads that cover this basepair (in this case 24793). However, pileupcolumn.pileups iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting min_base_quality parameter when you call pileup.

To have a similar filter applied to both coverage and reads, you can use piliecolumn.get_num_aligned() to get the number of segments and pileupcolumn.pileups to get the segments themselves.