NM tags

[wwood]

Hi,

With the bioconda package in hand, I was able to put something together that passes at least cursory tests in CoverM.

However, I had to add a call to samtools callmd along these lines:

x-mapper{} --num-threads {} --reference '{}' {} --out-sam - | samtools calmd -S - '{}'

because it seems that x-mapper does not report NM tags in its sam output, which is what coverm uses to post-process the alignments and e.g. remove low identity matches according to user specifications.

This seems like double handling, because I guess x-mapper already calculates it internally - is it straightforward to add that to x-mapper output?

Thanks, ben

[mathjeff]

Hi

Thanks for the update!

Yeah, theoretically we could do that if it's helpful. Some questions, though:

• What would be the advantage of using NM, which I believe is the number of differences between the query and reference, instead of the alignment score, AS? The alignment score takes into account the user's specifications for the cost of an indel, mutation, and ambiguity in a sequence. Mapper also provides a custom Sam field, cs, which is the combined alignment score of the mates in a paired-end-read and incorporates the spacing between the alignments (including potential overlap).

• Mapper also provides its own filter for alignment penalty via --max-penalty <fraction>, and it uses its knowledge of that filter to allow it to run more quickly by skipping checking alignments that can easily be proven to be worse than that threshold - would it make more sense to specify that instead?

Thanks

[wwood]

Hi,

Good question - there are 2 reasons I can think of.

1. % ID and % read mapped are very straightforward to interpret. For instance, aligning reads at 95% is related to how well reads from one species maps to another in the same genus. People seem to like that simplicity. A percentage also removes the bias of longer reads having higher alignment scores than shorter ones.
2. Having NM would slot mapper into the coverm codebase in the same way that other mapping software fits.

CoverM has MAPQ-based filtering now as well, maybe would make sense to add AS/CS-based filtering as well for those interested.

WDYT?

[mathjeff]

Oh interesting. Is the idea that a user might be interested in knowing both what fraction of an individual read was aligned to the reference, and also how similar the aligned fraction of the read was to the reference? For example: 90% of this read was aligned, and that section of the read was 95% similar to the reference?

Also interesting: Mapper's alignment score filter is proportional to the length of the read: the default is --max-penalty 0.1, which corresponds to alignment score >= -0.1 * read length, which represents 90% similarity (Mapper's alignment scores are always <= 0). As a result, if a user is only interested in alignments having >= 95% similarity to the reference, it would usually be easiest to specify --max-penalty 0.05 (which would run more quickly and generate less output than filtering afterward).

It does make sense to me for Mapper to output what CoverM needs if we can figure out exactly what that is

[wwood]

Oh interesting. Is the idea that a user might be interested in knowing both what fraction of an individual read was aligned to the reference, and also how similar the aligned fraction of the read was to the reference? For example: 90% of this read was aligned, and that section of the read was 95% similar to the reference?

Yes, but also there is more possibilities for filtering. For instance, it is possible to specify % identity / % of read that is aligned across a alignment pair, or a hard number of bases the read needs to be aligned, etc. There is probably too many possbilities to translate them all. There is a full list here
https://wwood.github.io/CoverM/coverm-genome.html#alignment-thresholding

I'm still a bit confused with how your suggestion would completely stand-in for NM tag interpretation though. If the user wants at least 90% of the read to be aligned, so we set --max-penalty 0.1, then couldn't there be reads which align for >90%, but with imperfect identity, which would also fail that criteria? So then mapper wouldn't report them, even they would not be filtered out by using mapper with default params and reading the NM tags. Hopefully that makes sense.

Thanks, ben

[mathjeff]

Ok, I think you're saying that it's inconvenient to pass some of these filters directly to Mapper because you might have several different types of filters, including the absolute aligned length of the read, the absolute aligned length of the read plus its mate, and the fraction of the aligned length of the read.

I think this brings us back to the question of what the NM tag should be:

What is the desired behavior in case of an insertion, deletion, ambiguity, mismatch, soft clip, or paired-end read? It seems that there is a spec at https://samtools.github.io/hts-specs/SAMtags.pdf that appears to be not very confident about its own authority, saying "Note that historically this has been ill-defined and both data and tools exist that disagree with this definition."

[wwood]

I think if we cannot rely on even the Sam spec, then we are all lost. It is the best we have as a standard output of different tools. NM is output by several e.g. minimap2, BWA, strobealign, etc. and so to me is a reasonably well established convention. I'm not sure what tools the spec is referring to, or how old that statement is.

Is the SAM spec unclear in some way? You've probably given these kinds of things more thought than I have, so maybe I'm missing something?

[mathjeff]

Ok thanks. Here's the summary of my current understanding:

CoverM would like to know the number of aligned basepairs in each alignment in Mapper's .sam files so that this information can be compared to a user-specified filter, which might involve the absolute number of aligned basepairs, the fraction of aligned basepairs, or the total number of aligned basepairs in the query pair.

There are a few options for how this could be calculated, but currently the way that CoverM calculates this is by reading the CIGAR string and the NM field and doing some arithmetic to determine the number of unmutated basepairs.

CoverM would like for Mapper to output a value for the NM field as specified by the Sam spec so that it can apply its existing calculation to Mapper's output too.

Does that seem like a good description?

I think the Sam spec seems clear and specific but surprising and maybe not convenient for users.
It says each mismatch should count as 1 (reasonable), each insertion or deletion should count as 1 (slightly surprising), soft clip should be ignored (slightly surprising), and ambiguous code such as "R" and "N" should count as 1 (surprising) unless the corresponding CIGAR character is "=" (also surprising).

I think this means that in the future, if we upgrade Mapper's CIGAR flags to specify "="/"X" rather than "M", then any ambiguities in the reference will suddenly no longer be counted towards the NM value.

We could add a command line argument such as --nm to Mapper which adds the NM Sam flag and also forces the Sam CIGAR fields to contain "M" (like they do now) instead of "="/"X".

What do you think?