GMAP SAM output has an error when converting to BAM

[rkizen]

Hello,

I am using GMAP to map cDNA sequences but when I try to convert the SAM files to BAM I get this error:

Line 25706, sequence length 519 vs 499 from CIGAR
Parse error at line 25706: CIGAR and sequence length are inconsistent

I checked this sequence in the SAM file and there is indeed a discrepancy between the length of the sequence and the CIGAR.

This is the line from the SAM file that is causing the problem:

m140206_015639_42149_c100624452550000001823118308061485_s1_p0/61037/33_552_CCS 16 AADN03010840.1 1984 40 5S4M1I1M1I1M1I1M2D1M1I1M1I12M10D27M7I6M6D2M8I3M2D2M2D5M2I5M9D7M2D23M1D8M1I1M1107N8M1I12M1I48M1I56M1I41M1D28M1D6M728N78M5334N1I6M1I43M101N28M1S * 0 0 CTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCTTTTTTTTTTTTTTGTTTTTTTTCTCTTTCTTTTTTTTTTTTCGTTTTTTTCACAGACCAGCGACGTGTTATTACATCTACCGCAGTACAGGTTGTATTTGCCGAGCGCCCTTGGGGTGGGGGCTGCTGTCGGCTGCCCTTGTAGGAGCCAGAATTGCTGCTGCTGCTCAGCGAGGAGGTGGCCGTAGGTGTCCCTCCCAGGGGGCGAGGGGCGCCTGGAGGCAGGGAGAAGAGGCTGGTGTCTCAGGCCGAGTGGCTGAGGGCGCTGACGACTTAAAGAGGAGGCAAAACCACTCCGTCCATAAAGCCGCGACTTCGTGATCCGGTGGATGAGCGTCTGCTCCCCGGCAGAGACATGGCCACGGTGCCGCGGGTCCTGGGCACAGGGGACACTGCGCGGGGGCCCGCGGCTCTCAGGGGCGGCAGCGCCGGGCTCAGCGACCTGCGCCGGCCTCCATGGGTGAGGCGGCC * MD:Z:1GCA0G0C0C^GC1C0C1GCG1G2C2^CAGGCGGGCC2CAGG1GCCG1A1GCACAC3C4GC1G1^AGGGGG5^GT2^TG1TC7^GGTGGTTTG2GC3^GG0T1GT1A1T1C13^G95G78^A28^C83A70C6 NH:i:1 HI:i:1 NM:i:110 SM:i:40 XQ:i:40 X2:i:0 XO:Z:UU XS:A:- PG:Z:M

Has anyone encountered this issue before? What is causing GMAP to produce an incorrect CIGAR for that sequence?

[Brian Bushnell]

Bugs happen; I suggest you report it to the authors and remove that read.

[rkizen]

Right it seems the error has been brought up before and the newest release (2014-05-15.v2) was supposed to fix it.

I installed it and ran it and got the same errors on different reads. Actually there were more errors this time. The release I had used earlier was gmap-2014-04-06.

I wrote a script to remove the erroneous lines from the sam files but I am still curious as to what is causing this error and if there may be some hidden errors associated with it.

I will send a report to the developer.

[Brian Bushnell]

If you are encountering multiple problematic reads, you could try BBMap as an alternative. It's one of the few aligners designed for DNA and RNA sequencing of both Illumina and PacBio reads, and it does not require any reference annotation. A typical command (for PacBio RNA-seq) would be:

mapPacBio8k in=reads.fa out=mapped.sam ref=organism.fa maxindel=200000

Please let me know if you encounter any problems!