I created a much simpler version of MethylDackel to extract read counts and found some differences.
After careful debugging, I discovered that this statement was the culprit:
//Is the phred score even high enough?
if(bam_get_qual(plp->b)[plp->qpos] < config->minPhred) return 0;
That is line 126 of extract.c
Basically, bam_get_qual(plp->b)[plp->qpos) was always zero.
As it happens, the optimizing compiler was generating incorrect code.
On the other hand, the call to this function was:
rv = updateMetrics(config, plp[0]+i);
The quality check could happen here, rather than inside updateMetrics, where it worked perfectly.
I wonder if the same will happen with a different compiler or with optimizations disabled.
I created a much simpler version of MethylDackel to extract read counts and found some differences.
After careful debugging, I discovered that this statement was the culprit:
That is line 126 of extract.c
Basically, bam_get_qual(plp->b)[plp->qpos) was always zero.
As it happens, the optimizing compiler was generating incorrect code.
On the other hand, the call to this function was:
The quality check could happen here, rather than inside updateMetrics, where it worked perfectly.
I wonder if the same will happen with a different compiler or with optimizations disabled.