Handle alignments without quality.#2083
Conversation
| if (bam_get_qual(p->b)[p->qpos] < opt_min_baseQ) { | ||
| --depth_at_pos; // low base quality | ||
| } else { | ||
| stats[tid].summed_baseQ += bam_get_qual(p->b)[p->qpos]; |
There was a problem hiding this comment.
While not new logic, this means that the "Mean BaseQ" reported by the tool is the average quality of bases which meet the minimum base quality criteria, rather than the average quality of all bases. This isn't what I imagined it to be and the man page is rather ambiguous given it doesn't explicitly lean one way or the other.
It's probably too late to change the program functionality, unless it's a clear bug (but depth_at_pos is also filtered depth). The man page problem just needs some tweakage to explain that the filtering takes place before computing averaged statistics.
There was a problem hiding this comment.
Alignments that don't have the minimum base quality don't get included in the coverage calculations. So it is the average quality of used alignments. I think the man page covers that.
This change does mess with the calculations, hence the warning. It is a compromise so that the program doesn't have to exit when it comes across a missing quality value.
|
From discussions in the office I was of the opinion this was a missing bounds check in the mpileup algorithm, but this is purely samtools coverage code. If I misunderstood, then great as this is clearly the less bad scenario of the two. Can you confirm that the same issues don't appear upstream anywhere in htslib? |
This is a purely samtools change. The errors were caused by bam_get_qual being called in the samtools coverage main loop on non-existent quality data. |
| } | ||
|
|
||
| if (print_value_warning) { | ||
| print_error("coverage", "Warning: Missing quality values in alignments. Mean base quality calculated only on available values.\n"); |
There was a problem hiding this comment.
I'm not sure if it's deliberate or not to add an extra line, but print_error already prints a newline. It's hardly unique in this front though. It looks like we have 18% of calls adding their own explicit newline.
There was a problem hiding this comment.
You know, I even checked what print_error did and didn't notice that.
There was a problem hiding this comment.
I did notice that nearly all other uses of print_error with a newline come from bam_markdup ;-) [Guilty too of a couple in amplicon_stats]
|
Looks good to me, barring the question about whether we want a double or single newline in the error message. It's kind of irrelevant, but thought I'd ask. |
Secondary alignments with CIGAR strings but without quality (valid in the sam spce) could cause segfaults and/or uninitialised memory access. This commit fixes that and changes the quality calculations to just use existing qualities. This changes existing behaviour but as secondary values are excluded by default the effects should be minor.
c1bf945 to
aa745f4
Compare
Secondary alignments with CIGAR strings but without quality (valid in the sam spec) could cause segfaults and/or uninitialised memory access in samtools coverage.
The presence of a CIGAR string means they have length (for the purposes of coverage) but they have no sequence or quality data. The guards on trying to access non-existent quality data were imperfect.
This commit fixes the faulty memory access and changes the quality calculations to just use existing qualities. This changes existing behaviour but as secondary values are excluded by default the effects should be minor. A warning is printed if there are missing qualities.
This should fix #2076.