text: user1770: I'm writing a python function from scratch to do this, but I feel like this must exist in some standard bioinformatics library already. In principle, this is a simply regex operation which many must have written previously. With the goal of having a centralized reference as opposed to everyone writing their own scripts to accomplish the same task (which is time-consuming and not necessarily reproducible), here is the "construct a CIGAR" function:
  Here's an alignment: 
  ACGT-TGC
  
  Here is a reference:
  ACGTATGC
  
  The reference will be required in order to record indels. One could in principle have a series of "queries", e.g. an R vector or a Python list. 
  This would generate the following CIGAR: 
  4=1D3=
  
  Which R packages/ python library exists to do this operation? I suspect there is something related to Bioconductor as well.
  
  user1770: Here is what I would do in Python. There are a few extra checks I've included, and my definitions of a few cases are different than Heng Li's above, @user172818. The performance is fairly similar. 
  def create_cigar(reference, alignment):
      '''Construct CIGAR from reference and alignment'''
      ## if reference and alignment have unequal lengths, stop
      if len(reference) != len(alignment):
          raise Exception('Error: reference and alignment are of unequal length')
      ## replace any '-' and '*' to '_'
      reference = reference.translate(bytes.maketrans(b"*-", b"__"))
      alignment = alignment.translate(bytes.maketrans(b"*-", b"__"))
      ### creates string with a CIGAR symbol for each position using zip()
      compare_strings = ''.join(['=', 'X'][ref == '_'] if ref == algn else 'X' if ref != '_' and algn != '_' else ['I', 'D'][algn == '_' and ref != '_'] for ref, algn in zip(reference, alignment))
      ## now, parse string into CIGAR format; faster than itertools
      total_match = 0         # total number of adjacent matches
      current_base = compare_strings[0]     # current character being processed
      result = ''     # result of function
      for i in range(len(compare_strings)):
          if compare_strings[i] == current_base:
              total_match += 1
          else:
              result += str(total_match) + current_base
              current_base = compare_strings[i]
              total_match = 1
      result += str(total_match) + current_base
      return result
  
  Following the example:
  print(create_cigar('ACGTAT-CT', 'ACGT-TGGA'))
  
  which outputs:
  4=1D1=1I2X
  
  (1) I standard the inputs, such that gaps can be either -, _, or *. (Users could add their own.)
  (2) There's no N, skipped regions, though based on the use case, I could see users writing this, e.g. a deletion >= size some number could work. 
  (3) Hard and soft padding do not make sense in this context. These will be classified as X. 
  (4) By my definition above, if there's a gap in the reference and the alignment, this is a padding P. 
  
  user37: Your request is easy to implement. I wouldn't use any libraries in this case. As you haven't showed your implementation, I will provide one:
  def gen_cigar(ref, qry):
      if len(ref) != len(qry):
          raise Exception('unequal length')
      cigar = []
      for i in range(len(ref)):
          r, q = ref[i], qry[i];
          if r is '-' and q is '-':
              raise Exception('both gaps')
          op = '=' if r is q else 'I' if r is '-' else 'D' if q is '-' else 'X';
          if len(cigar) > 0 and cigar[-1][1] is op: # add to the last operation
              cigar[-1][0] += 1
          else: cigar.append([1, op]);              # a new operation
      return "".join(map(lambda x: str(x[0]) + x[1], cigar)); # turn to string
  
  print(gen_cigar('ACGTAT-CT', 'ACGT-TGGA'));
  
  The output is:
  4=1D1=1I2X
  
  PS: It is interesting that you often approach problems from a different angle. The solutions that you think "must exist" actually rarely exist because those are not the ways we (me at least) work with data. This is why I prefer an asker to explain the bigger picture in addition to isolated intermediate problems – there could be simpler and more straightforward solutions without the need of solving the contrived problems.
  
original_text: 1770: <p>I'm writing a python function from scratch to do this, but I feel like this must exist in some standard bioinformatics library already. In principle, this is a simply regex operation which many must have written previously. With the goal of having a centralized reference as opposed to everyone writing their own scripts to accomplish the same task (which is time-consuming and not necessarily reproducible), here is the "construct a CIGAR" function:</p>
  
  <p>Here's an alignment: </p>
  
  <pre><code>ACGT-TGC
  </code></pre>
  
  <p>Here is a reference:</p>
  
  <pre><code>ACGTATGC
  </code></pre>
  
  <p>The reference will be required in order to record indels. One could in principle have a series of "queries", e.g. an R vector or a Python list. </p>
  
  <p>This would generate the following CIGAR: </p>
  
  <pre><code>4=1D3=
  </code></pre>
  
  <p>Which R packages/ python library exists to do this operation? I suspect there is something related to Bioconductor as well.</p>
   1770: <p>Here is what I would do in Python. There are a few extra checks I've included, and my definitions of a few cases are different than Heng Li's above, @user172818. The performance is fairly similar. </p>
  
  <pre><code>def create_cigar(reference, alignment):
      '''Construct CIGAR from reference and alignment'''
      ## if reference and alignment have unequal lengths, stop
      if len(reference) != len(alignment):
          raise Exception('Error: reference and alignment are of unequal length')
      ## replace any '-' and '*' to '_'
      reference = reference.translate(bytes.maketrans(b"*-", b"__"))
      alignment = alignment.translate(bytes.maketrans(b"*-", b"__"))
      ### creates string with a CIGAR symbol for each position using zip()
      compare_strings = ''.join(['=', 'X'][ref == '_'] if ref == algn else 'X' if ref != '_' and algn != '_' else ['I', 'D'][algn == '_' and ref != '_'] for ref, algn in zip(reference, alignment))
      ## now, parse string into CIGAR format; faster than itertools
      total_match = 0         # total number of adjacent matches
      current_base = compare_strings[0]     # current character being processed
      result = ''     # result of function
      for i in range(len(compare_strings)):
          if compare_strings[i] == current_base:
              total_match += 1
          else:
              result += str(total_match) + current_base
              current_base = compare_strings[i]
              total_match = 1
      result += str(total_match) + current_base
      return result
  </code></pre>
  
  <p>Following the example:</p>
  
  <pre><code>print(create_cigar('ACGTAT-CT', 'ACGT-TGGA'))
  </code></pre>
  
  <p>which outputs:</p>
  
  <pre><code>4=1D1=1I2X
  </code></pre>
  
  <p>(1) I standard the inputs, such that gaps can be either <code>-</code>, <code>_</code>, or <code>*</code>. (Users could add their own.)</p>
  
  <p>(2) There's no <code>N</code>, skipped regions, though based on the use case, I could see users writing this, e.g. a deletion >= size some number could work. </p>
  
  <p>(3) Hard and soft padding do not make sense in this context. These will be classified as <code>X</code>. </p>
  
  <p>(4) By my definition above, if there's a gap in the reference and the alignment, this is a padding <code>P</code>. </p>
   37: <p>Your request is easy to implement. I wouldn't use any libraries in this case. As you haven't showed your implementation, I will provide one:</p>
  
  <pre class="lang-python prettyprint-override"><code>def gen_cigar(ref, qry):
      if len(ref) != len(qry):
          raise Exception('unequal length')
      cigar = []
      for i in range(len(ref)):
          r, q = ref[i], qry[i];
          if r is '-' and q is '-':
              raise Exception('both gaps')
          op = '=' if r is q else 'I' if r is '-' else 'D' if q is '-' else 'X';
          if len(cigar) &gt; 0 and cigar[-1][1] is op: # add to the last operation
              cigar[-1][0] += 1
          else: cigar.append([1, op]);              # a new operation
      return "".join(map(lambda x: str(x[0]) + x[1], cigar)); # turn to string
  
  print(gen_cigar('ACGTAT-CT', 'ACGT-TGGA'));
  </code></pre>
  
  <p>The output is:</p>
  
  <pre><code>4=1D1=1I2X
  </code></pre>
  
  <p>PS: It is interesting that you often approach problems from a different angle. The solutions that you think "must exist" actually rarely exist because those are not the ways we (me at least) work with data. This is why I prefer an asker to explain the bigger picture in addition to isolated intermediate problems – there could be simpler and more straightforward solutions without the need of solving the contrived problems.</p>
  