Hi there, first of all, thanks for making bcftools! It's really amazing to have these kind of tools available.
I find the documentation around the mpileup option --adjust-MQ a bit lacking though and the option seems to sometimes give unexpected results.
-
Just from the documentation, I would expect the adjusted mapping quality to follow the given formula:
$adjust = \sqrt{\frac{option - quality}{option}}\cdot option$,
which would result in the following relationship between mapping quality and adjusted mapping quality (assuming that mapping quality values larger than the value given are left as is):
this however seems counterintuitive, as it would boost really bad aligning reads and reduce the value for better aligning ones (or is this what is intended?)
-
Looking at the code it seems like bcftools calculates the mapping quality used for the adjustment independently of the reported mapping quality from the aligner (see https://github.com/samtools/htslib/blob/develop/realn.c#L39 which is what is called in mpileup). I have not looked at this specifically in illumina data, but I have an example in 10x data where all reads at a position have a reported mapping quality of 60, but some of them seem to get adjusted to below 20 (using -C 50 -q 20), judging from the resulting DP from mpileup. It seems like there is no documentation on how this quality is calculated (and I have not managed to fully understand it from the code) and for which type of data this is applicable.
I'm happy to contribute to the documentation to make this option a bit clearer once I actually understand how it works.
Hi there, first of all, thanks for making bcftools! It's really amazing to have these kind of tools available.
I find the documentation around the mpileup option --adjust-MQ a bit lacking though and the option seems to sometimes give unexpected results.
Just from the documentation, I would expect the adjusted mapping quality to follow the given formula:
$adjust = \sqrt{\frac{option - quality}{option}}\cdot option$ ,
which would result in the following relationship between mapping quality and adjusted mapping quality (assuming that mapping quality values larger than the value given are left as is):
this however seems counterintuitive, as it would boost really bad aligning reads and reduce the value for better aligning ones (or is this what is intended?)
Looking at the code it seems like bcftools calculates the mapping quality used for the adjustment independently of the reported mapping quality from the aligner (see https://github.com/samtools/htslib/blob/develop/realn.c#L39 which is what is called in mpileup). I have not looked at this specifically in illumina data, but I have an example in 10x data where all reads at a position have a reported mapping quality of 60, but some of them seem to get adjusted to below 20 (using -C 50 -q 20), judging from the resulting DP from mpileup. It seems like there is no documentation on how this quality is calculated (and I have not managed to fully understand it from the code) and for which type of data this is applicable.
I'm happy to contribute to the documentation to make this option a bit clearer once I actually understand how it works.