Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-12-20T05:36:25.087Z
Answers.Text: great reference code with the caveat that MD is indeed really annoying to work with. oftentimes easier to just recreate the mismatches from the sequence itself
Answers.UserDisplayName: Colin D
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-02T06:29:51.457Z
Answers.LastActivityDate: 2018-03-02T06:29:51.457Z
Answers.UserDisplayName: user172818
Answers.Body: <p>This program in C++ is able to perform the task <a href="https://github.com/mlafave/sam2pairwise" rel="nofollow noreferrer">https://github.com/mlafave/sam2pairwise</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-26T01:59:37.483Z
Answers.LastActivityDate: 2018-11-26T01:59:37.483Z
Answers.UserDisplayName: Colin D
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-09T17:58:55.840Z
Answers.Text: Is there a standard method to grab the string from the BAM within python? It looks like this is parsed automatically with pysam
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-12T10:24:00.130Z
Answers.Text: I guess pysam would be the standard way
Answers.UserDisplayName: heathobrien
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-01T18:42:45.407Z
Answers.LastActivityDate: 2018-03-01T18:42:45.407Z
Answers.UserDisplayName: heathobrien
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-01T08:56:38.190Z
Comments.Text: Could you clarify your question: You want to store the alignments in a matrix? Or the CIGAR format in a matrix?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-01T14:28:32.567Z
Comments.Text: @Llopis A tsv/csv file with two columns: the first column has the alignments, the second column has the sequence. Does this make sense?
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-01T14:51:41.553Z
Comments.Text: Do you want sam2pairwise to generate the alignments for you, or do you want to parse the cigar strings yourself in R/python? Writing a cigar parser is tedious, but do-able. It might be easier if you just use sam2pairwise to output the alignments for each read, then you parse that into a table.
Comments.UserDisplayName: conchoecia
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-01T19:25:25.010Z
Comments.Text: Now your question is reduced to "how to count the number of insertions and deletions" in cigar, which is much easier. With perl, it is just an one-liner `echo 78M2I10M5D30M|perl -ne 's/(\d+)([ID])/$h{$2}+=$1/eg; print "$h{I}\t$h{D}\n"'`. Someone will probably post python answers.
Comments.UserDisplayName: user172818
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-01T21:56:56.953Z
Comments.Text: @user172818 I've tried to edit this to be more specific. Is this more organized? The second problem is easier; I'm primarily motivated by the first problem, which is related.
Comments.UserDisplayName: ShanZhengYang
Tags: <r><python><sam><phylogenetics>
Title: Generating the reconstructed alignment from BAM
ThreadText: # Generating the reconstructed alignment from BAM
  
  (Asked by: username_0 on 2018-03-01)
  
  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`
  
  **Tags:** r, python, sam, phylogenetics
  
  > **Comments:**
  > **username_1**: Could you clarify your question: You want to store the alignments in a matrix? Or the CIGAR format in a matrix?  
  > **username_0**: @Llopis A tsv/csv file with two columns: the first column has the alignments, the second column has the sequence. Does this make sense?  
  > **username_2**: Do you want sam2pairwise to generate the alignments for you, or do you want to parse the cigar strings yourself in R/python? Writing a cigar parser is tedious, but do-able. It might be easier if you just use sam2pairwise to output the alignments for each read, then you parse that into a table.  
  > **username_3**: Now your question is reduced to "how to count the number of insertions and deletions" in cigar, which is much easier. With perl, it is just an one-liner `echo 78M2I10M5D30M|perl -ne 's/(\d+)([ID])/$h{$2}+=$1/eg; print "$h{I}\t$h{D}\n"'`. Someone will probably post python answers.  
  > **username_0**: @username_3 I've tried to edit this to be more specific. Is this more organized? The second problem is easier; I'm primarily motivated by the first problem, which is related.  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_5 on 2018-03-01. Score: 5 | Accepted answer
  
  If all you want are the number of insertions/deletions, you can use [this](https://pypi.python.org/pypi/cigar/0.1) 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'`
  
  > **Comments:**
  > **username_0**: Is there a standard method to grab the string from the BAM within python? It looks like this is parsed automatically with pysam  
  > **username_5**: I guess pysam would be the standard way  
  
  ---
  
  ### Answer by username_3 on 2018-03-02. Score: 5
  
  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.
  
  ```js
  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]);
  
  ```
  
  > **Comments:**
  > **username_4**: great reference code with the caveat that MD is indeed really annoying to work with. oftentimes easier to just recreate the mismatches from the sequence itself  
  
  ---
  
  ### Answer by username_4 on 2018-11-26. Score: 0
  
  This program in C++ is able to perform the task <https://github.com/mlafave/sam2pairwise>
  
  
  ---
  
  