# Is there a straightforward way to get mismatches/indels in a BAM file using pysam?

Ideally, I'd like to do something like the following:

f = pysam.AlignmentFile(filename, index_filename=index_filename)
read = list(f.fetch('Chromosome', 1, 10000))[0]
read.get_variations()


And have that output something like:

[
(1, 'A', 'C'), # An A -> C mismatch at position 1
(15, 'I', 'CGGG'), # An insertion of the string 'CGGG' starting at position 15 of the read
(32, 'D', 4), # A deletion of 4 base pairs starting at position 32
]


I know that it's possible to figure out all of those values by parsing the CIGAR string and MD tags but I'm wondering if there's a more intuitive interface.

## 1 Answer

You're looking for pileup, which is the htslib (and thus samtools/bcftools) method for finding variants. See this section of the pysam documentation. You may want to peruse about the output as well. The gist is:

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
for pileupread in pileupcolumn.pileups:
if not pileupread.is_del and not pileupread.is_refskip:
# query position is None if is_del or is_refskip is set.
...handle mismatches...

• Thanks for the response! I think this would work but it's less than ideal since I would be iterating over each read multiple times, no? I'm wondering if it's possible to iterate over the list of reads that span a region and retrieve that same information without going over pileup columns. Dec 16 '20 at 16:12
• Also, this bit of code for pileupread in pileupcolumn.pileups crashes the Python kernel when running inside of an IPython notebook :-/ Dec 16 '20 at 16:17