bowtie CIGAR and MapQ errors

[skingan]

Hello,

I am mapping illumina reads of small RNAs to the Drosophila genome using bowtie1. My command:

bowtie -S -n 0 -l 15 -k 200 --best ref.fa reads.fq > file.sam

The map qualities and cigars are uniformly too good for the mapped reads. The cigars are all LM where L is the length of the read, despite in some cases there being mismatches as indicate in the NM field. MapQ is always 255 if the read is mapped.

Has anyone ever seen this? Is bowtie just outputting meaningless cigars and mapQs?

I am thinking of using bwa instead of bowtie.

Thanks in advance.

[skingan]

This thread answers the mapQ question:

http://seqanswers.com/forums/showthread.php?t=7815

But I'm still confused about the cigar.

[dpryan]

A CIGAR of LM (for L==some number) is still correct in the presence of mismatches. M in this case simply means "mapped to a base in the reference", not "matches a base in the reference". Compare "M", "=" and "X" and see the SAM spec.

[skingan]

Thanks for your reply, dpryan. But if I have a non-zero edit distance, shouldn't there be some X's in the cigar? Here are some fields from my file:

MapQ: 255
CIGAR: 30M
NM:i:2

And an alignment between the read and the reference sequence

TATTACTTGTTTTTCTTTGCGACTTCGGGT
TATTACTTGTTTTTCTTTGCGACTTCGTTT

That cigar is not right but the NM field is.

[dpryan]

No, there shouldn't be. Again, reread what I wrote and the specification. An M CIGAR operation does not imply a sequence match between the query and reference. Very few aligners use X CIGAR operations (BBMap does and STAR are the only ones that I know of off-hand), so you have to parse the MD tag.

[skingan]

Okay, I see it now. Thanks!