5

To give you a bit of context: I am trying to convert a sam file to bam

samtools view -bT reference.fasta sequences.sam > sequences.bam

which exits with the following error

[E::sam_parse1] CIGAR and query sequence are of different length
[W::sam_read1] parse error at line 102
[main_samview] truncated file

and the offending line looks like this:

SRR808297.2571281       99      gi|309056|gb|L20934.1|MSQMTCG   747     80      101M    =       790     142     TTGGTATAAAATTTAATAATCCCTTATTAATTAATAAACTTCGGCTTCCTATTCGTTCATAAGAAATATTAGCTAAACAAAATAAACCAGAAGAACAT      @@CFDDFD?HFDHIGEGGIEEJIIJJIIJIGIDGIGDCHJJCHIGIJIJIIJJGIGHIGICHIICGAHDGEGGGGACGHHGEEEFDC@=?CACC>CCC      NM:i:2  MD:Z:98A1A

My sequence is 98 characters long but a probable bug when creating the sam file reported 101 in the CIGAR. I can give myself the luxury to loss a couple of reads and I don't have access at the moment to the source code that produced the sam files, so no opportunity to hunt down the bug and re-run the alignment. In other words, I need a pragmatic solution to move on (for now). Therefore, I devised a python script that counts the length of my string of nucleotides, compares it with what is registered in the CIGAR, and saves the "sane" lines in a new file.

#!/usr/bin/python
import itertools
import cigar

with open('myfile.sam', 'r') as f:
    for line in itertools.islice(f,3,None): #Loop through the file and skip the first three lines
            cigar=line.split("\t")[5]
            cigarlength = len(Cigar(cigar)) #Use module Cigar to obtain the length reported in the CIGAR string
            seqlength = len(line.split("\t")[9])

            if (cigarlength == seqlength):
                    ...Preserve the line in a new file...

As you can see, to translate the CIGAR into an integer showing the length, I am using the module CIGAR. To be honest, I am a bit wary of its behavior. This module seems to miscalculate the length in very obvious cases. Is there another module or a more explicit strategy to translate the CIGAR into the length of the sequence?

Sidenote: Interesting, to say the least, that this problem has been widely reported but no pragmatic solution can be found in the internet. See the links below:

https://github.com/COMBINE-lab/RapMap/issues/9
http://seqanswers.com/forums/showthread.php?t=67253
http://seqanswers.com/forums/showthread.php?t=21120
https://groups.google.com/forum/#!msg/snap-user/FoDsGeNBDE0/nRFq-GhlAQAJ
1
  • "To be honest, I am a bit wary of its behavior. This module seems to miscalculate the length in very obvious cases." - such as? The underlying code is short and simple; it's hard to envisage a way in which it could fail.
    – Mark Amery
    Commented Jun 28, 2019 at 14:14

2 Answers 2

5

The SAM spec offers us this table of CIGAR operations which indicates which ones "consume" the query or the reference, complete with explicit instructions on how to calculate sequence length from a CIGAR string:

                                                             Consumes  Consumes
Op  BAM Description                                             query  reference
M   0   alignment match (can be a sequence match or mismatch)   yes   yes
I   1   insertion to the reference                              yes   no
D   2   deletion from the reference                             no    yes
N   3   skipped region from the reference                       no    yes
S   4   soft clipping (clipped sequences present in SEQ)        yes   no
H   5   hard clipping (clipped sequences NOT present in SEQ)    no    no
P   6   padding (silent deletion from padded reference)         no    no
=   7   sequence match                                          yes   yes
X   8   sequence mismatch                                       yes   yes
  • “Consumes query” and “consumes reference” indicate whether the CIGAR operation causes the alignment to step along the query sequence and the reference sequence respectively.

...

  • Sum of lengths of the M/I/S/=/X operations shall equal the length of SEQ.

This lets us trivially calculate the length of a sequence from its CIGAR by adding up the lengths of all the "consumes query" ops in the CIGAR. This is exactly what happens in the cigar module (see https://github.com/brentp/cigar/blob/754cfed348364d390ec1aa40c951362ca1041f7a/cigar.py#L88-L93), so I don't know why the OP here reckoned that module's implementation was wrong.

If we extract out the relevant code from the (already very short) cigar module, we're left with something like this as a short Python implementation of the summing operation described in the quote above:

from itertools import groupby

def query_len(cigar_string):
    """
    Given a CIGAR string, return the number of bases consumed from the
    query sequence.
    """
    read_consuming_ops = ("M", "I", "S", "=", "X")
    result = 0
    cig_iter = groupby(cigar_string, lambda chr: chr.isdigit())
    for _, length_digits in cig_iter:
        length = int(''.join(length_digits))
        op = next(next(cig_iter)[1])
        if op in read_consuming_ops:
            result += length
    return result
3

I suspect the reason there isn't a tool to fix this problem is because there is no general solution, aside from performing the alignment again using software that does not exhibit this problem. In your example, the query sequence aligns perfectly to the reference and so in that case the CIGAR string is not very interesting (just a single Match operation prefixed by the overall query length). In that case the fix simply requires changing 101M to 98M.

However, for more complex CIGAR strings (e.g. those that include Insertions, Deletions, or any other operations), you would have no way of knowing which part of the CIGAR string is too long. If you subtract from the wrong part of the CIGAR string, you'll be left with a misaligned read, which is probably worse for your downstream analysis than just leaving the whole read out.

That said, if it happens to be trivial to get it right (perhaps your broken alignment procedure always adds extra bases to the first or last CIGAR operation), then what you need to know is the correct way to calculate the query length according to the CIGAR string, so that you know what to subtract from it.

samtools calculates this using the htslib function bam_cigar2qlen.

The other functions that bam_cigar2qlen calls are defined in sam.h, including a helpful comment showing the truth table for which operations consume query sequence vs reference sequence.

In short, to calculate the query length of a CIGAR string the way that samtools (really htslib) does it, you should add the given length for CIGAR operations M, I, S, =, or X and ignore the length of CIGAR operations for any of the other operations.

The current version of the python cigar module seem to be using the same set of operations, and the algorithm for calculating the query length (which is what len(Cigar(cigar)) would return) looks right to me. What makes you think that it isn't giving the correct results?

It looks like you should be able to use the cigar python module to hard clip from either the left or right end using the mask_left or mask_right method with mask="H".

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.