1
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.