Answers.Body: <p>You're looking for <code>pileup</code>, which is the htslib (and thus samtools/bcftools) method for finding variants. See <a href="https://pysam.readthedocs.io/en/latest/api.html#introduction" rel="nofollow noreferrer">this section</a> of the pysam documentation. You may want to peruse about <a href="https://pysam.readthedocs.io/en/latest/api.html#pysam.PileupColumn" rel="nofollow noreferrer">the output</a> as well. The gist is:</p>
  <pre><code>import pysam
  samfile = pysam.AlignmentFile(&quot;ex1.bam&quot;, &quot;rb&quot; )
  for pileupcolumn in samfile.pileup(&quot;chr1&quot;, 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...
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-12-16T16:12:57.067Z
Answers.Text: 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.
Answers.UserDisplayName: juniper-
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-12-16T16:17:06.587Z
Answers.Text: Also, this bit of code `for pileupread in pileupcolumn.pileups` crashes the Python kernel when running inside of an IPython notebook :-/
Answers.UserDisplayName: juniper-
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-12-16T07:49:43.570Z
Answers.LastActivityDate: 2020-12-16T07:49:43.570Z
Answers.UserDisplayName: Devon Ryan
Body: <p>Ideally, I'd like to do something like the following:</p>
  <pre><code>f = pysam.AlignmentFile(filename, index_filename=index_filename)
  read = list(f.fetch('Chromosome', 1, 10000))[0]
  read.get_variations()
  </code></pre>
  <p>And have that output something like:</p>
  <pre><code>[
    (1, 'A', 'C'), # An A -&gt; 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
  ]
  </code></pre>
  <p>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.</p>
  
Tags: <sam>
Title: Is there a straightforward way to get mismatches/indels in a BAM file using pysam?
ThreadText: # Is there a straightforward way to get mismatches/indels in a BAM file using pysam?
  
  (Asked by: username_0 on 2020-12-16)
  
  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.
  
  **Tags:** sam
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2020-12-16. Score: 1
  
  You're looking for `pileup`, which is the htslib (and thus samtools/bcftools) method for finding variants. See [this section](https://pysam.readthedocs.io/en/latest/api.html#introduction) of the pysam documentation. You may want to peruse about [the output](https://pysam.readthedocs.io/en/latest/api.html#pysam.PileupColumn) 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...
  
  ```
  
  > **Comments:**
  > **username_0**: 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.  
  > **username_0**: Also, this bit of code `for pileupread in pileupcolumn.pileups` crashes the Python kernel when running inside of an IPython notebook :-/  
  
  ---
  
  