# Is there a Python/R package with the ability to convert an alignment and reference into a CIGAR?

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.

• Maybe GenomicAlignments ? Apr 4 '18 at 22:05
• @zx8754 Do you have a function in mind? Apr 4 '18 at 23:58
• @gringer Thanks. I've tried to edit the above to make this more clear. Apr 5 '18 at 0:04
• 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
Apr 5 '18 at 3:20

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.

• 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.... Apr 9 '18 at 19:26
• 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. If you really want to use is for efficiency, please intern the strings explicitly beforehand. Apr 24 '18 at 11:35

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.