Extracting percent identity of mapped reads and the coordinates on the genome the read mapped from a BAM file?

Question

[Eric]

Hello,

I am working with bam files that have environmental reads mapped back to an isolate genome. I am looking for a way to extract both the percent identity of all mapped reads in the BAM file and the coordinates on the genome where they mapped. I was looking at some combination of maybe bedtools and pysam or parsing the CIGAR string but nothing I have tried has worked so far.

I was trying to find a pythonic answer, but if people have other solutions feel free to share.

Thanks!

Answer 1

[simon.burgermeister]

If you use pysam:

samfile = pysam.AlignmentFile("your_bam", "rb")
for read in samfile.fetch():
     cig=read.cigarstring #for cigarstring
     aligned=read.get_aligned_pairs() # will give you the aligned pairs

Be careful 'M' in cigar also include SNPs

[Eric]

Thanks for this example, its very helpful being new to pysam!