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
.