3
$\begingroup$

I have a list containing the results of a MSA:

msa = [
    '-G-ATG-CGGATACCG-', # reference sequence
    'CGTATCTCG--TACCGA',
    '-GTATCTCGT-GACCG-',
    'CG-AT--CGG-GACCGA',
    '---AT--CG--GACCGA',
    '-GTAT---GC-GACCG-',
    'CG-AT--CGG-GACCGA',
    'CGTAT--CG--GA CG-',
]

I would like to generate the CIGAR strings of those alignments (against the reference sequence); the expected result would be something like:

cigars = [
    '13M',
    '1I1M1I3M1I2M2D5M1I',
    '1M1I3M1I3M1D5M',
    '1I3M1D3M1D5M1I',
    '1D2M1D2M2D5M1I',
    '1M1I2M2D2M1D5M',
    '1I3M1D3M1D5M1I',
    '1I1M1I2M1D2M2D5M',
]

I don't know what would be a good approach for this. I managed to get the "aligned pairs" of each alignment but I don't think that'll help for determining the cigar strings:

# mapping = position in alignment => position in reference
mapping = {y: x for x,y in enumerate(y for y,nuc in enumerate(msa[0]) if nuc != '-')}

for aln in msa:
    aligned_pairs = [(y, mapping[z]) for y,z in enumerate(z for z,nuc in enumerate(aln) if nuc != '-') if z in mapping]
    print(aligned_pairs)
[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12)]
[(1, 0), (3, 1), (4, 2), (5, 3), (7, 4), (8, 5), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(0, 0), (2, 1), (3, 2), (4, 3), (6, 4), (7, 5), (8, 6), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(0, 1), (1, 2), (2, 4), (3, 5), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12)]
[(0, 0), (2, 1), (3, 2), (4, 5), (5, 6), (6, 8), (7, 9), (8, 10), (9, 11), (10, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(1, 0), (3, 1), (4, 2), (5, 4), (6, 5), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]

How can I generate the CIGAR strings efficiently?

$\endgroup$

1 Answer 1

4
$\begingroup$

If you loop through the characters of two aligned sequences in parallel then you can easily determine the operations that compose the CIGAR. Then you just have to "collapse" the generated list of operations:

import itertools

ref = msa[0]

for aln in msa:
    cig = []

    for i,nuc in enumerate(aln):
        if nuc != '-':
            if ref[i] != '-':
                cig.append('M')
            else:
                cig.append('I')
        elif ref[i] != '-':
                cig.append('D')

    cig = [(len(list(it)), op) for op,it in itertools.groupby(cig)]

    print( ''.join(str(n) + op for n,op in cig) )
13M
1I1M1I3M1I2M2D5M1I
1M1I3M1I3M1D5M
1I3M1D3M1D5M1I
1D2M1D2M2D5M1I
1M1I2M2D2M1D5M
1I3M1D3M1D5M1I
1I1M1I2M1D2M2D5M
$\endgroup$

Your Answer

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

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