Infer the length of a sequence using the CIGAR
<p>To give you a bit of context: I am trying to convert a sam file to bam</p>

<pre><code>samtools view -bT reference.fasta sequences.sam &gt; sequences.bam
</code></pre>

<p>which exits with the following error</p>

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

<p>and the offending line looks like this:</p>

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

<p>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.</p>

<pre><code>#!/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...
</code></pre>

<p>As you can see, to translate the CIGAR into an integer showing the length, I am using the module <a href="https://pypi.python.org/pypi/cigar/0.1" rel="nofollow">CIGAR</a>. 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?</p>

<p><strong>Sidenote:</strong> 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:</p>

<pre><code>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
</code></pre>

<p>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 <code>M</code>atch operation prefixed by the overall query length). In that case the fix simply requires changing <code>101M</code> to <code>98M</code>. </p>

<p>However, for more complex CIGAR strings (e.g. those that include <code>I</code>nsertions, <code>D</code>eletions, 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. </p>

<p>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. </p>

<p><code>samtools</code> calculates this using the <code>htslib</code> function <a href="https://github.com/samtools/htslib/blob/19c189438f852e6e62dbda73f854d465cebb3d9f/sam.c#L325-L332" rel="nofollow">bam_cigar2qlen</a>.</p>

<p>The other functions that <code>bam_cigar2qlen</code> calls are defined in <a href="https://github.com/samtools/htslib/blob/bf753361dab9b1640cf64f7886dbfe35357a43c5/htslib/sam.h#L76-L105" rel="nofollow">sam.h</a>, including a helpful comment showing the truth table for which operations consume query sequence vs reference sequence.</p>

<p>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 <code>M</code>, <code>I</code>, <code>S</code>, <code>=</code>, or <code>X</code> and ignore the length of CIGAR operations for any of the other operations. </p>

<p>The current version of the python cigar module seem to be using the <a href="https://github.com/brentp/cigar/blob/master/cigar.py#L68" rel="nofollow">same set of operations</a>, and the algorithm for calculating the query length (which is what <code>len(Cigar(cigar))</code> would return) looks right to me. What makes you think that it isn't giving the correct results?</p>

<p>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 <code>mask_left</code> or <code>mask_right</code> method with <code>mask="H"</code>.</p>

python module bioinformatics samtools