<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>

CC BY-SA 4.0
2023-12-20T05:36:25.087Z
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
Colin D
CC BY-SA 3.0
2018-03-02T06:29:51.457Z
2018-03-02T06:29:51.457Z
user172818
<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>

CC BY-SA 4.0
2018-11-26T01:59:37.483Z
2018-11-26T01:59:37.483Z
Colin D
<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>

CC BY-SA 3.0
2018-03-09T17:58:55.840Z
Is there a standard method to grab the string from the BAM within python? It looks like this is parsed automatically with pysam
ShanZhengYang
CC BY-SA 3.0
2018-03-12T10:24:00.130Z
I guess pysam would be the standard way
heathobrien
CC BY-SA 3.0
2018-03-01T18:42:45.407Z
2018-03-01T18:42:45.407Z
heathobrien
<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>

CC BY-SA 3.0
2018-03-01T08:56:38.190Z
Could you clarify your question: You want to store the alignments in a matrix? Or the CIGAR format in a matrix?
llrs
CC BY-SA 3.0
2018-03-01T14:28:32.567Z
@Llopis A tsv/csv file with two columns: the first column has the alignments, the second column has the sequence. Does this make sense?
ShanZhengYang
CC BY-SA 3.0
2018-03-01T14:51:41.553Z
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.
conchoecia
CC BY-SA 3.0
2018-03-01T19:25:25.010Z
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.
user172818
CC BY-SA 3.0
2018-03-01T21:56:56.953Z
@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.
ShanZhengYang
<r><python><sam><phylogenetics>
Generating the reconstructed alignment from BAM
# 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**: @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.  

---

## 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 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.
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>


---

