5
$\begingroup$

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.

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

2 Answers 2

5
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ 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.... $\endgroup$
    – EB2127
    Apr 9, 2018 at 19:26
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 24, 2018 at 11:35
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.