Promote secondary when primary alignment fails to surject

[gaoj66-roche]

Changelog Entry

To be copied to the draft changelog by merger:

• Add --promote-secondary to vg surject, giraffe, map, and mpmap for promoting a read's secondary alignment to primary alignment when the primary alignment fails to surject (becomes unmapped)

Description

When a read's primary alignment fails to surject (becomes unmapped) but a secondary surjects, the best-scoring mapped, non-supplementary secondary is promoted to primary instead of emitting an unmapped primary alongside mapped secondaries, which can be considered a violation of the BAM/SAM spec and cause issues with variant calling in regions where a duplication event is present in the pangenome but not the linear reference. This PR adds a CLI option to handle this situation. When enabled, we will "promote" the surjectable secondary alignment and replace the unsurjectable (unmapped) primary alignment with an auxiliary tag ps (promoted secondary) to indicate that it is promoted from a secondary alignment.

Only alignments already present are considered with --max-multimaps > 1 (and the user will be warned if --promote-secondary is enabled without --max-multimaps > 1).

For standalone vg surject, input must be collated by read name. A grouped parallel reader keeps each read's primary and secondaries in the same group within a batch, never splitting a group across worker threads. In order to support this new type of parallelism, a group-aware parallel iterator is added to libvgio (See vgteam/libvgio#92). vg surject with GAF input is not supported because GAF does not handle secondary alignments and it is not possible to distinguish a secondary alignment from a primary alignment in a GAF file.

One important caveat for the vg surject command is that the input has to be collated (alignments with the same query name must be adjacent to each other). I admit that this is a bit awkward since we don't have a tool for sorting GAM by query name so the only way to get a collated GAM right now is directly from output of the aligner. Once a GAM file is sorted by graph position, it would not possible to recover the collated GAM if the user deletes the original aligner output. I might add the ability to sort alignments by query name or collate alignments to gamsort in the future to keep the scope of this PR focused on secondary promotion during surjection.

[gaoj66-roche]

A few questions and thoughts

• This initially came up when we are developing our copy number caller where we noticed BAMs aligned with vg giraffe would miss common duplication calls when the duplication is present in the pangenome reference but not in the linear reference and we realized that the root cause is because reads supporting the duplication would map to the duplication node on the graph but become unmapped after surjection. We developed our own tool for rescuing these reads but it has a pretty siginificant IO overhead by having to wait for giraffe to finish, read the file, write to separate file, and sort the final BAM. I believe that it should be handled by vg instead of a separate cleanup/rescue pass.
• BAM spec says "When one segment is present in multiple lines to represent a multiple mapping of the segment, only one of these records should have the secondary alignment flag bit (0x100) unset. RNEXT and PNEXT point to the primary line of the next read in the template." Although it does not explicitly prohibit unmapped primary with mapped secondary, it does say that secondary alignments should be used in the case of "multiple mapping" so an unmapped primary does not really make sense because it is not a "mapping". Moreover, the convention seems to be choosing the best alignment as primary and set everything else as secondary. Even though in our case, the original primary alignment is indeed the best alignment in the graph space, it is obviously not the best alignment in the linear space and we should be basing our decisions on what happens in the linear space during surjection (another example: alignment score of surjected record should be calculated based on alignment to the linear reference, not graph alignment).
• Currently, secondary alignments do not carry meaningful mapping qualities. But I think if a secondary alignment is promoted to a primary, it should have a meaningful mapping quality. There are a few ways one can implement this and I would like to know what you think is the best option.
  • Calculate mapping quality on the fly during surjection when a secondary is promoted. The issue is we don't have all the information at the surjection stage in order to calculate the mapping quality using the same algorithm. In particular, we keep only max_multimaps many alignments but mapq is calculated using all alignments of a read.
  • Calculate and store mapping quality for secondary alignment records; only output pre-calculated mapping quality when a secondary is promoted; otherwise, set mapq back to 0. This would work fine with BAM/SAM/CRAM output but it would also mean secondary alignments in GAM/GAF gets a mapq. This is a deviation from the current convention.
  • Let downstream tools use the alignment score as a proxy for mapping quality. This is what we are currently doing in our tools but requires prior knowledge of how the secondary promotion process works and would mean that even though a secondary is "promoted" to a primary, it still would not behave exactly like a primary.

[adamnovak]

I think promoting the secondary alignment to the primary surjection is the right idea, at least in some workflows like copy number calling where you don't actually want the read to be at its true location when the true location is not in the linear reference.

@Sagorikanag in our lab has been thinking about this in the context of surjecting to a diploid reference, and has some ideas on handling mapping quality in a way that can encompass both uncertainty around which graph placement is the true source and uncertainty around which surjection is the true source. Maybe she has ideas here?

[Sagorikanag]

I have been working on a related approach for diploid surjection (Mapping to a pangenome graph containing a donor specific diploid assembly (DSA) and linear references, then surjecting back to the diploid DSA) that may also help with the promoted-secondary MAPQ problem.

I group all consecutive GAM/GAF records with the same read name so that the primary and retained secondary graph placements are processed together. Each graph placement is then surjected to all overlapping target paths (Anything that is set as REF in the graph, for me that is the DSA).

For every individual graph placement, I calculate a haplotype quality (HQ) by comparing the scores of that placement's haplotype-path surjections. I also add haplotype tags identifying the preferred and alternative haplotype surjections, currently represented as HP:Z:pri_hap and HP:Z:sec_hap.

After that, I pool the surjections produced from all graph placements for the read and calculate a separate global quality across the complete retained candidate set. This global step selects exactly one overall primary surjection and marks the remaining non-supplementary candidates as secondary. The original graph mapping quality is also retained in a separate tag before the output MAPQ is replaced. (Right now we cap the Global Quality at the reads' Mapping Quality, but that is something we're working on) So in this step we see many secondary graph placements yielding us the global primary surjection including when the primary graph placement remains unsurjected.

For paired-end reads, the implementation also incorporates the estimated fragment-length distribution when evaluating the paired candidates.

This would not recreate the mapper's original MAPQ exactly, because at the surjection stage we still only have the retained max_multimaps candidates rather than every candidate originally considered by the mapper. Instead, it produces a new surjection-level quality from the candidate set that is actually available.

Would you be interested in this?

[adamnovak]

@gaoj66-roche I also checked with @Sagorikanag and as part of the work she's been doing here at UCSC on her thesis project, involving diploid-reference surjection and these new notions of quality, she also has her own implementation of the grouped IO stuff that you're proposing in vgteam/libvgio#92, to support her in-progress implementation of a diploid mode for vg surject.

I haven't seen her code recently, but apparently her grouped input implementation in libvgio handles both GAF and GAM and shares a little more code with the existing codepaths, so I think it might be better to take her implementation instead, depending on what it looks like, and then get your feature here to sit on top of that.

It might be possible that, if we can land her diploid surject changes in vg soon, the whole --promote-secondary feature could sit on top of or automatically fall out of the logic she has for juggling the set of all surjections from a group of graph alignments, instead of needing to be its own pass.

I talked to her and she's going to get that work up as PRs as soon as she can. Then we can read through everything and figure out how to merge everything together in the most parsimonious way. I want to get you all at Roche the feature you need, without creating a jillion merge conflicts between our grad student and her thesis.

[gaoj66-roche]

Thanks both! I am happy to put this PR on hold till @Sagorikanag 's grouping iterator implementation and everything else is merged and rebase my PR with that. I also agree with @adamnovak that the MAPQ should ideally capture the confidence of both the mapping in the graph space and the linear space but it could potentially be a larger change than what is necessary to address the current issue of secondary with unmapped primary.

In my current PR, I also added a tag to the primary alignments promoted from a secondary due to the original primary failing surjection so that downstream tools can decide if they want to treat them differently or even use them at all. I think it would be useful (at least until we can figure out how to assign a consistent and meaningful MAPQ to promoted secondaries) to have some way of distinguishing promoted primaries from a real primary.

Another possiblity would be to not strip the secondary flag at all and only append a tag to let downstream tools decide whether it wants to treat it as a secondary or primary. This is probably the lowest-risk approach since many variant callers already discard secondary alignments by default so the secondary promotion feature would be less likely to cause any regression in other tools that use vg aligned BAMs as input, unless they deliberately decide to use the promoted alignments.

Please let me know what you all think and keep me posted on any discussion/decision related to this. Thanks!