Counting uniquely mapped reads - ambiguous results

Question

[mickey_95]

Hi,

I have bam files aligned with Star, which according to the manual assigns a MAPQ score of 255 to all uniquely mapped reads. Using

samtools view -c -q 255 myfile.bam

I should obtain the number of uniquely mapped reads in my file. The same should be obtained using:

samtools view -c -F 260 myfile.bam

This is, however, not the case and the read numbers are quite different: ~95 mio for the first command and ~350 mio for the latter. How can this be? Am I overlooking something?

I would really appreciate your help!

Answer 1

[ATpoint]

The first command filters for mapping quality, the second one for reads being Not Unmapped and Not Primary, therefore being mapped and primary. Still, a multimapper can be a primary alignment. Therefore you are comparing apples with peers.

[mickey_95]

Thanks for the quick response! Just to be sure I am understanding this correctly: If a read is mapping to multiple loci, the read with the best alignment score would be defined as primary alignment and all others would be secondary, correct?

[ATpoint]

From what I know: If a read maps to multiple locations equally well it gets a MAPQ of 0. Typically from all possible locations one is randomly chosen and reported as the primary alignment, if you instruct the aligner to output more than one location then all other locations would be secondary alignments. Still, your flag Not Primary includes also supplementary and chimeric alignments, so read that could have something to do with structural variants etc. THerefore it is not really a fair comparison.

[mickey_95]

Ah, okay! Thanks for clarifying!