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

Question

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.

Answer 1

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.

Answer 2

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.