Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-10-01T01:45:39.707Z
Answers.LastActivityDate: 2023-10-02T19:45:47.970Z
Answers.LastEditDate: 2023-10-02T19:45:47.970Z
Answers.UserDisplayName: Fravadona
Body: <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>
  
Tags: <python><cigar>
Title: Generating CIGAR strings for MSA with reference
ThreadText: # 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:
  
  ```py
  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:
  
  ```py
  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:
  
  ```py
  # 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)
  
  ```
  
  ```py
  [(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:
  
  ```py
  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) )
  
  ```
  
  ```py
  13M
  1I1M1I3M1I2M2D5M1I
  1M1I3M1I3M1D5M
  1I3M1D3M1D5M1I
  1D2M1D2M2D5M1I
  1M1I2M2D2M1D5M
  1I3M1D3M1D5M1I
  1I1M1I2M1D2M2D5M
  
  ```
  
  
  ---
  
  