mpileup with & without -q: odd result?

[afkoeppel]

Hello everyone,

I'm running mpileup twice on the same BAM alignment of RNASeq reads against the same reference genome, with and without filtering for mapping quality (-q 5 or not), and getting, (I think) an odd result.

I checked the same position on the same chromosome between the filtered and unfiltered output.
In both cases, 7997 reads overlapped that position, but the counts of reads that matched the reference (,.) versus those that were mismatches (Aa) was drastically different between the two. Specifically, there were far more mismatches at that position in the filtered output than in the unfiltered output.

If I'm understanding what -q does correctly, it should simply be disregarding alignments with a mapq score of lower than 5 (if I use -q 5), right? So I would understand if the number of reads overlapping the position was reduced, but I don't understand why the count of overlapping reads would remain the same, but the proportion of matches to mismatches would change.

Can anyone help me make sense of this?

Thanks,
Alex

P.S.
samtools 0.1.18
Linux OS

[BAMseek]

Samtools has a maximum depth cut-off in order to limit the amount of memory required, which leads to some odd behavior when you hit that max. I think the max is 8,000 and it looks like you are close to that. Try supplying the -d flag with a really large number (like 10 million), and see if the results are more what you expect.

Justin

[afkoeppel]

That seems to have done the trick. Thanks.