<p>If you loop through the characters of two aligned sequences <strong>in parallel</strong> then you can easily determine the operations that compose the CIGAR. Then you just have to &quot;collapse&quot; the generated list of operations:</p>
<pre class="lang-py prettyprint-override"><code>import itertools

ref = msa[0]

for aln in msa:
    cig = []

    for i,nuc in enumerate(aln):
        if nuc != '-':
            if ref[i] != '-':
                cig.append('M')
            else:
                cig.append('I')
        elif ref[i] != '-':
                cig.append('D')

    cig = [(len(list(it)), op) for op,it in itertools.groupby(cig)]

    print( ''.join(str(n) + op for n,op in cig) )
</code></pre>
<pre class="lang-py prettyprint-override"><code>13M
1I1M1I3M1I2M2D5M1I
1M1I3M1I3M1D5M
1I3M1D3M1D5M1I
1D2M1D2M2D5M1I
1M1I2M2D2M1D5M
1I3M1D3M1D5M1I
1I1M1I2M1D2M2D5M
</code></pre>

CC BY-SA 4.0
2023-10-01T01:45:39.707Z
2023-10-02T19:45:47.970Z
2023-10-02T19:45:47.970Z
Fravadona
<p>I have a list containing the results of a MSA:</p>
<pre class="lang-py prettyprint-override"><code>msa = [
    '-G-ATG-CGGATACCG-', # reference sequence
    'CGTATCTCG--TACCGA',
    '-GTATCTCGT-GACCG-',
    'CG-AT--CGG-GACCGA',
    '---AT--CG--GACCGA',
    '-GTAT---GC-GACCG-',
    'CG-AT--CGG-GACCGA',
    'CGTAT--CG--GA CG-',
]
</code></pre>
<p>I would like to generate the CIGAR strings of those alignments (against the reference sequence); the expected result would be something like:</p>
<pre class="lang-py prettyprint-override"><code>cigars = [
    '13M',
    '1I1M1I3M1I2M2D5M1I',
    '1M1I3M1I3M1D5M',
    '1I3M1D3M1D5M1I',
    '1D2M1D2M2D5M1I',
    '1M1I2M2D2M1D5M',
    '1I3M1D3M1D5M1I',
    '1I1M1I2M1D2M2D5M',
]
</code></pre>
<p>I don't know what would be a good approach for this. I managed to get the &quot;aligned pairs&quot; of each alignment but I don't think that'll help for determining the cigar strings:</p>
<pre class="lang-py prettyprint-override"><code># mapping = position in alignment =&gt; position in reference
mapping = {y: x for x,y in enumerate(y for y,nuc in enumerate(msa[0]) if nuc != '-')}

for aln in msa:
    aligned_pairs = [(y, mapping[z]) for y,z in enumerate(z for z,nuc in enumerate(aln) if nuc != '-') if z in mapping]
    print(aligned_pairs)
</code></pre>
<pre class="lang-py prettyprint-override"><code>[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12)]
[(1, 0), (3, 1), (4, 2), (5, 3), (7, 4), (8, 5), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(0, 0), (2, 1), (3, 2), (4, 3), (6, 4), (7, 5), (8, 6), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(0, 1), (1, 2), (2, 4), (3, 5), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12)]
[(0, 0), (2, 1), (3, 2), (4, 5), (5, 6), (6, 8), (7, 9), (8, 10), (9, 11), (10, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(1, 0), (3, 1), (4, 2), (5, 4), (6, 5), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
</code></pre>
<p>How can I generate the CIGAR strings efficiently?</p>

<python><cigar>
Generating CIGAR strings for MSA with reference
# Generating CIGAR strings for MSA with reference

(Asked by: username_0 on 2023-09-30)

I have a list containing the results of a MSA:
msa = [
    '-G-ATG-CGGATACCG-', # reference sequence
    'CGTATCTCG--TACCGA',
    '-GTATCTCGT-GACCG-',
    'CG-AT--CGG-GACCGA',
    '---AT--CG--GACCGA',
    '-GTAT---GC-GACCG-',
    'CG-AT--CGG-GACCGA',
    'CGTAT--CG--GA CG-',
]

I would like to generate the CIGAR strings of those alignments (against the reference sequence); the expected result would be something like:
cigars = [
    '13M',
    '1I1M1I3M1I2M2D5M1I',
    '1M1I3M1I3M1D5M',
    '1I3M1D3M1D5M1I',
    '1D2M1D2M2D5M1I',
    '1M1I2M2D2M1D5M',
    '1I3M1D3M1D5M1I',
    '1I1M1I2M1D2M2D5M',
]

I don't know what would be a good approach for this. I managed to get the "aligned pairs" of each alignment but I don't think that'll help for determining the cigar strings:
# mapping = position in alignment => position in reference
mapping = {y: x for x,y in enumerate(y for y,nuc in enumerate(msa[0]) if nuc != '-')}

for aln in msa:
    aligned_pairs = [(y, mapping[z]) for y,z in enumerate(z for z,nuc in enumerate(aln) if nuc != '-') if z in mapping]
    print(aligned_pairs)

[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12)]
[(1, 0), (3, 1), (4, 2), (5, 3), (7, 4), (8, 5), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(0, 0), (2, 1), (3, 2), (4, 3), (6, 4), (7, 5), (8, 6), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(0, 1), (1, 2), (2, 4), (3, 5), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12)]
[(0, 0), (2, 1), (3, 2), (4, 5), (5, 6), (6, 8), (7, 9), (8, 10), (9, 11), (10, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(1, 0), (3, 1), (4, 2), (5, 4), (6, 5), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]

How can I generate the CIGAR strings efficiently?

**Tags:** python, cigar

---

## 1 answers

### Answer by username_0 on 2023-10-01. Score: 4 | Accepted answer

If you loop through the characters of two aligned sequences in parallel then you can easily determine the operations that compose the CIGAR. Then you just have to "collapse" the generated list of operations:
import itertools

ref = msa[0]

for aln in msa:
    cig = []

    for i,nuc in enumerate(aln):
        if nuc != '-':
            if ref[i] != '-':
                cig.append('M')
            else:
                cig.append('I')
        elif ref[i] != '-':
                cig.append('D')

    cig = [(len(list(it)), op) for op,it in itertools.groupby(cig)]

    print( ''.join(str(n) + op for n,op in cig) )

13M
1I1M1I3M1I2M2D5M1I
1M1I3M1I3M1D5M
1I3M1D3M1D5M1I
1D2M1D2M2D5M1I
1M1I2M2D2M1D5M
1I3M1D3M1D5M1I
1I1M1I2M1D2M2D5M


---

