Meaning of Negative Alignment Score in BOWTIE2 output

[foolishbrat]

I have a mIRNA seq data in which I aligned it to whole genome UCSC hg19.
For that I use bowtie2 with the following parameter:

bowtie2 -x myindex -U my_mirna_seq.fastq --n-ceil L,0,0 -k 1 -N 1  -S mydata.sam

However in the SAM output I notice the alignment score is negative:

$cut -f12  mydata.sam | grep AS | sort -u | tail 
AS:i:-29
AS:i:-3
AS:i:-30
AS:i:-4
AS:i:-5
AS:i:-6
AS:i:-7
AS:i:-8
AS:i:-9
AS:i:0

What does it mean to have -ve AS value?
How should I deal with it?

And also the MAPQ score is 255 for all the alignments, no others.
Is that normal?

[dpryan]

Regarding the negative alignment scores, the highest possible alignment score in end-to-end alignment is 0. For each mismatch or gap open/extension, that score is decremented by some amount (for mismatches, this amount depends on the phred score of the base).

Regarding the MAPQ value, this is due to you using "-k 1". What you're telling bowtie2 to do is to look for at most one valid hit and to stop looking once it's found one. That hit will often not be the best hit that bowtie2 could find and since it's being told to stop looking right away it can't possibly calculate a MAPQ, since it would need to search more thoroughly to have any clue how good the returned alignment is. So, 255 would be the appropriate value.

Just don't use "-k 1" and you'll get the single best alignment for each read and an accompanying MAPQ score.

[foolishbrat]

Thanks, two more questions:

So if I remove -k1 that means for each read it will output only 1 best position, which doesn't mean it's uniquely mapped, am I right?

If so how can I determine the uniquely mapped reads from the parameter or the SAM output?
By looking at MAPQ=255?

[dpryan]

Correct, there's no option to tell bowtie2 to only report unique alignments. For bowtie2 (I should note that this is ONLY true for bowtie2), MAPQ scores of 0 and 1 indicate non-unique mapping, so you can filter by that (though not if you use "-k 1", since this only works when you have actual MAPQ values). Alternatively (also not when "-k 1" is used, but a bit more portable to other aligners), you can often just look for reads with "NH:i:1" auxiliary tags, though that's more work.

[foolishbrat]

Thanks a million.