TransWikia.com

Find all the bases for given reference position

Bioinformatics Asked by Diesel__100 on August 8, 2020

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?

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP