Use MD/NM on CRAM for samtools consensus#2156
Conversation
This is awkward as there's no right solution for BAM vs CRAM, but the previous method of not computing MD for CRAM was both illogical and also removing user control. Now it basically relies on --input-fmt-option decode_md=1 (or =0) which is adjustable, defaulting to 1 (decode). So BAM vs CRAM may differ now if the original BAM has no MD. Previously it differed if the original BAM had an MD. The latter is much more common place, hence it was a bad choice. TODO: Sadly CRAM doesn't have a mechanism to record whether the original data has MD/NM tags, so it always generates or always doesn't generate. Consider adding a private htslib specific tag to hold this CRAM meta-data, or later extending the specification with additional CF (cram flag) bits. We coould copy the CRAMv4 idea of a tag indicating the aux tag ordering, so if it contains MD and NM then we know to decode them, as well as which position in the tag string (along with RG too). That's essentially how CRAMv4 works but via the TL dictionary.
|
For the curious, the reason the MD/NM tags affect consensus generation is by default, consensus uses the mapping quality to adjust the base qualities when computing the pileup. Mapping quality however is woefully inaccurate for long read data, therefore we also have a localised mapping quality adjustment based on the number of mismatches within the proximity of this base. A large number of nearby mismatches reduces the local mapping quality and in turn the base quality. Particularly useful on platforms with varying levels of "lucidity"! A similar effect is applied for mismatches near the end of an alignment or close to a soft-clip position. It's all heuristics, but ones demonstrated to work. Obviously this requires the MD and NM tags to be present, so adding them when they weren't there originally or removing them when they were there will change consensus unless the |
This is awkward as there's no right solution for BAM vs CRAM, but the previous method of not computing MD for CRAM was both illogical and also removing user control.
Now it basically relies on --input-fmt-option decode_md=1 (or =0) which is adjustable, defaulting to 1 (decode).
So BAM vs CRAM may differ now if the original BAM has no MD. Previously it differed if the original BAM had an MD. The latter is much more common place, hence it was a bad choice.
TODO: Sadly CRAM doesn't have a mechanism to record whether the original data has MD/NM tags, so it always generates or always doesn't generate. Consider adding a private htslib specific tag to hold this CRAM meta-data, or later extending the specification with additional CF (cram flag) bits. We could copy the CRAMv4 idea of a tag indicating the aux tag ordering, so if it contains MD and NM then we know to decode them, as well as which position in the tag string (along with RG too). That's essentially how CRAMv4 works but via the TL dictionary.