text: user146: I have a (small) BAM file with CIGAR and MD fields. 
  Question 1: What tools exists in Python and/or R to reconstruct the alignment between the reference and the read in a BAM? Given that this is a very standard analysis, I feel that there should be R packages or Python libraries which have functions which do this task given a BAM input...
  My goal was to create a *tsv file with two columns, one with the reconstructed alignment and the other column with the corresponding sequence. This would be useful in my analysis. (I would either use R data.table, or pandas in Python.) 
  Here is read and reference:
  CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
  |||||||||||||||| ||||||||||||||||| 
  CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG
  
  Here is the desired output dataframe (with only the one example row). We could use an R data.table/data.frame, or a pandas DataFrame. Column1 are the alignments (i.e. the reconstruction of the alignment between the reference and the read based on the CIGAR and MD strings), Column 2 are the corresponding sequences. :
  ReadAligned    Reference   
  CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG    CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
  ...    ....
  
  Question 2:
  Using the CIGAR/MD tags, one could then create columns with the number of insertions and the number of deletions. This problem feels much more straightforward, e.g. the python library cigar or the cigar parsing functions of pysam
  
  user37: I am answering question 1. However, as I am unfamiliar with python, I am providing a javascript solution only to show the basic logic. To try the code at the end, you need to install node.js and run it with node this-script.js.
  
  Given that this is a very standard analysis
  
  I tend to believe this is a very rare analysis. I have only seen it done twice. Converting CIGAR+MD to padded alignment is complicated in coding, as you have to keep sequence, cigar and MD in sync all the time. The code below generates two strings in one go. It may be logically easier to generate the reference and the query strings separately.
  In general, MD is ill defined. It is too difficult to work with. If possible, avoid this tag. You can do a lot of things with the combination of CIGAR and NM.
  function md2pad(CIGAR, MD, SEQ) {
      var m, re_cigar = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g;
      var cigar = [];
      while ((m = re_cigar.exec(CIGAR)) != null) // parse CIGAR into an array
          if (m[2] != 'H') // do nothing for hard clipping, as it has no effect
              cigar.push([parseInt(m[1]), m[2]]);
      var k = 0, sx = '', sy = ''; // k: k-th cigar operator; x: ref; y: query
      // cx/cy: start coordinate of cigar[k] on ref/query; mx/my: start of current MD op
      var cx = 0, cy = 0, mx = 0, my = 0;
      while ((m = re_MD.exec(MD)) != null) {
          if (m[2] != null) { // deletion from the reference
              var len = m[2].length - 1; // length of deletion
              sx += m[2].substr(1), sy += Array(len+1).join("-");
              mx += len, cx += len, ++k; // consume a deletion D from ref
          } else { // copy or mismatch
              var ml = m[1] != null? parseInt(m[1]) : 1; // length of this MD operation
              while (k < cigar.length && cigar[k][1] != 'D') {
                  var cl = cigar[k][0], op = cigar[k][1];
                  if (op == 'M' || op == 'X' || op == '=') {
                      if (my + ml < cy + cl) { // an MD ends in the middle of cigar M
                          if (ml > 0) {
                              if (m[3] != null) sx += m[3], sy += SEQ[my];
                              else sx += SEQ.substr(my, ml), sy += SEQ.substr(my, ml);
                          }
                          mx += ml, my += ml, ml = 0;
                          break;
                      } else {
                          var dl = cy + cl - my;
                          sx += SEQ.substr(my, dl), sy += SEQ.substr(my, dl);
                          cx += cl, cy += cl, ++k; // consume an M operator
                          mx += dl, my += dl, ml -= dl;
                      }
                  } else if (op == 'I') {
                      sy += SEQ.substr(cy, cl), sx += Array(cl+1).join("-");
                      cy += cl, my += cl, ++k; // consume an insertion I to ref
                  } else if (op == 'S') {
                      cy += cl, my += cl, ++k; // consume a soft-clipping S
                  } else throw Error("inconsistent MD"); // not working with N
              }
              if (ml != 0) throw Error("inconsistent MD");
          }
      }
      if (cx != mx || cy != my) throw Error("inconsistent MD");
      return [sx, sy];
  }
  
  var cigar = '16025H3M2I55M1D108M306H';
  var MD = '7T15C24T9^C0G2T6T16T3A1A7G21T17G6C1T6T5A4';
  var seq = 'GATCACAGGCCTATCACCCTATTAATCACTCACGGGAGCTCTCCATGCATCTGGTATTTTTTCGGGGGGGGATGCACGCGATAGCATCGCGGGCCGCTGGAACCGGAGCACCCTATGTCGCAGGATCTGTCTTTGATTCCTACCTCATGCCATTATTAATCGCGCCTA';
  var s = md2pad(cigar, MD, seq);
  console.log(s[0]);
  console.log(s[1]);
  
  
  user581: If all you want are the number of insertions/deletions, you can use this to extract that info:
  from cigar import Cigar
  c = c=Cigar('20M5I10M3I5M')
  num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
  insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])
  
  Deletions would be the same thing with 'D' in place of 'I'
  
original_text: 146: <p>I have a (small) BAM file with CIGAR and MD fields. </p>
  
  <p>Question 1: What tools exists in Python and/or R to reconstruct the alignment between the reference and the read in a BAM? Given that this is a very standard analysis, I feel that there should be R packages or Python libraries which have functions which do this task given a BAM input...</p>
  
  <p>My goal was to create a *tsv file with two columns, one with the reconstructed alignment and the other column with the corresponding sequence. This would be useful in my analysis. (I would either use R data.table, or pandas in Python.) </p>
  
  <p>Here is read and reference:</p>
  
  <pre><code>CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
  |||||||||||||||| ||||||||||||||||| 
  CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG
  </code></pre>
  
  <p>Here is the desired output dataframe (with only the one example row). We could use an R data.table/data.frame, or a pandas DataFrame. Column1 are the alignments (i.e. the reconstruction of the alignment between the reference and the read based on the CIGAR and MD strings), Column 2 are the corresponding sequences. :</p>
  
  <pre><code>ReadAligned    Reference   
  CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG    CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
  ...    ....
  </code></pre>
  
  <p>Question 2:
  Using the CIGAR/MD tags, one could then create columns with the number of insertions and the number of deletions. This problem feels much more straightforward, e.g. the python library <code>cigar</code> or the cigar parsing functions of <code>pysam</code></p>
   37: <p>I am answering question 1. However, as I am unfamiliar with python, I am providing a javascript solution only to show the basic logic. To try the code at the end, you need to install node.js and run it with <code>node this-script.js</code>.</p>
  
  <blockquote>
    <p>Given that this is a very standard analysis</p>
  </blockquote>
  
  <p>I tend to believe this is a very rare analysis. I have only seen it done twice. Converting CIGAR+MD to padded alignment is complicated in coding, as you have to keep sequence, cigar and MD in sync all the time. The code below generates two strings in one go. It may be logically easier to generate the reference and the query strings separately.</p>
  
  <p>In general, MD is ill defined. It is too difficult to work with. If possible, avoid this tag. You can do a lot of things with the combination of CIGAR and NM.</p>
  
  <pre class="lang-js prettyprint-override"><code>function md2pad(CIGAR, MD, SEQ) {
      var m, re_cigar = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g;
      var cigar = [];
      while ((m = re_cigar.exec(CIGAR)) != null) // parse CIGAR into an array
          if (m[2] != 'H') // do nothing for hard clipping, as it has no effect
              cigar.push([parseInt(m[1]), m[2]]);
      var k = 0, sx = '', sy = ''; // k: k-th cigar operator; x: ref; y: query
      // cx/cy: start coordinate of cigar[k] on ref/query; mx/my: start of current MD op
      var cx = 0, cy = 0, mx = 0, my = 0;
      while ((m = re_MD.exec(MD)) != null) {
          if (m[2] != null) { // deletion from the reference
              var len = m[2].length - 1; // length of deletion
              sx += m[2].substr(1), sy += Array(len+1).join("-");
              mx += len, cx += len, ++k; // consume a deletion D from ref
          } else { // copy or mismatch
              var ml = m[1] != null? parseInt(m[1]) : 1; // length of this MD operation
              while (k &lt; cigar.length &amp;&amp; cigar[k][1] != 'D') {
                  var cl = cigar[k][0], op = cigar[k][1];
                  if (op == 'M' || op == 'X' || op == '=') {
                      if (my + ml &lt; cy + cl) { // an MD ends in the middle of cigar M
                          if (ml &gt; 0) {
                              if (m[3] != null) sx += m[3], sy += SEQ[my];
                              else sx += SEQ.substr(my, ml), sy += SEQ.substr(my, ml);
                          }
                          mx += ml, my += ml, ml = 0;
                          break;
                      } else {
                          var dl = cy + cl - my;
                          sx += SEQ.substr(my, dl), sy += SEQ.substr(my, dl);
                          cx += cl, cy += cl, ++k; // consume an M operator
                          mx += dl, my += dl, ml -= dl;
                      }
                  } else if (op == 'I') {
                      sy += SEQ.substr(cy, cl), sx += Array(cl+1).join("-");
                      cy += cl, my += cl, ++k; // consume an insertion I to ref
                  } else if (op == 'S') {
                      cy += cl, my += cl, ++k; // consume a soft-clipping S
                  } else throw Error("inconsistent MD"); // not working with N
              }
              if (ml != 0) throw Error("inconsistent MD");
          }
      }
      if (cx != mx || cy != my) throw Error("inconsistent MD");
      return [sx, sy];
  }
  
  var cigar = '16025H3M2I55M1D108M306H';
  var MD = '7T15C24T9^C0G2T6T16T3A1A7G21T17G6C1T6T5A4';
  var seq = 'GATCACAGGCCTATCACCCTATTAATCACTCACGGGAGCTCTCCATGCATCTGGTATTTTTTCGGGGGGGGATGCACGCGATAGCATCGCGGGCCGCTGGAACCGGAGCACCCTATGTCGCAGGATCTGTCTTTGATTCCTACCTCATGCCATTATTAATCGCGCCTA';
  var s = md2pad(cigar, MD, seq);
  console.log(s[0]);
  console.log(s[1]);
  </code></pre>
   581: <p>If all you want are the number of insertions/deletions, you can use <a href="https://pypi.python.org/pypi/cigar/0.1" rel="noreferrer">this</a> to extract that info:</p>
  
  <pre><code>from cigar import Cigar
  c = c=Cigar('20M5I10M3I5M')
  num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
  insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])
  </code></pre>
  
  <p>Deletions would be the same thing with <code>'D'</code> in place of <code>'I'</code></p>
  