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

CC BY-SA 3.0
2018-04-09T19:24:03.817Z
2018-04-09T19:24:03.817Z
EB2127
<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>

CC BY-SA 3.0
2018-04-09T19:26:14.857Z
Hi Heng Li, thanks for the help. I've shared my own answer below. You're right that this isn't a difficult function to write; my motivation was largely to standardize these sorts of tasks. (e.g. apparently there's a standard MATLAB function for the above, `align2cigar`). Anyways, RE: PS...thanks for the compliment!....I think....
EB2127
CC BY-SA 3.0
2018-04-24T11:35:38.157Z
Your use of `is` instead of `==` gives me the heebie-jeebies. It’s absolutely not guaranteed that this will succeed. You got lucky due to string interning (that said, I tested it on CPython 2.6.9, 2.7.14, 3.6.4 and Pypy 5.10.0 and it works everywhere). But pylint correctly complains, and [it’s strongly recommended not to rely on this for correctness](https://stackoverflow.com/a/19899966/1968). **If** you really want to use `is` for efficiency, please intern the strings explicitly beforehand.
Konrad Rudolph
CC BY-SA 3.0
2018-04-05T14:50:04.103Z
2018-04-05T14:50:04.103Z
user172818
<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>

CC BY-SA 3.0
2018-04-04T22:05:06.757Z
Maybe [GenomicAlignments](https://bioconductor.org/packages/release/bioc/html/GenomicAlignments.html) ?
zx8754
CC BY-SA 3.0
2018-04-04T23:58:08.427Z
@zx8754 Do you have a function in mind?
EB2127
CC BY-SA 3.0
2018-04-05T00:04:36.170Z
@gringer Thanks. I've tried to edit the above to make this more clear.
EB2127
CC BY-SA 3.0
2018-04-05T03:20:53.650Z
Sorry, I realise that my previous statement was wrong; only two of the three (query, reference, match) are needed to make a CIGAR string, bearing in mind that gaps in the reference should be allowed. In any case, it would be useful having a demonstrative CIGAR string that included at least I,D,X,=.
gringer
<r><python><phylogenetics><bioconductor><cigar>
Is there a Python/R package with the ability to convert an alignment and reference into a CIGAR?
# Is there a Python/R package with the ability to convert an alignment and reference into a CIGAR?

(Asked by: username_0 on 2018-04-04)

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.

**Tags:** r, python, phylogenetics, bioconductor, cigar
> **Comments:**
> **username_1**: Maybe [GenomicAlignments](https://bioconductor.org/packages/release/bioc/html/GenomicAlignments.html) ?  
> **username_0**: @zx8754 Do you have a function in mind?  
> **username_0**: @gringer Thanks. I've tried to edit the above to make this more clear.  
> **username_2**: Sorry, I realise that my previous statement was wrong; only two of the three (query, reference, match) are needed to make a CIGAR string, bearing in mind that gaps in the reference should be allowed. In any case, it would be useful having a demonstrative CIGAR string that included at least I,D,X,=.  

---

## 2 answers

### Answer by username_3 on 2018-04-05. Score: 5 | Accepted answer

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.

> **Comments:**
> **username_0**: Hi Heng Li, thanks for the help. I've shared my own answer below. You're right that this isn't a difficult function to write; my motivation was largely to standardize these sorts of tasks. (e.g. apparently there's a standard MATLAB function for the above, `align2cigar`). Anyways, RE: PS...thanks for the compliment!....I think....  
> **username_4**: Your use of `is` instead of `==` gives me the heebie-jeebies. It’s absolutely not guaranteed that this will succeed. You got lucky due to string interning (that said, I tested it on CPython 2.6.9, 2.7.14, 3.6.4 and Pypy 5.10.0 and it works everywhere). But pylint correctly complains, and [it’s strongly recommended not to rely on this for correctness](https://stackoverflow.com/a/19899966/1968). **If** you really want to use `is` for efficiency, please intern the strings explicitly beforehand.  

---

### Answer by username_0 on 2018-04-09. Score: 3

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.


---

