vg surject : --diploid-map and other code changes

[Sagorikanag]

Purpose

This draft combines the complete diploid map surjection work into one branch for review.

This has not been rebased over #5025 yet. The paired-surjection implementation overlaps with that PR, and will be changed after reviewing by @adamnovak.

Main changes

• Add vg surject --diploid-map.
• Compare alternative graph placements jointly.
• Select haplotype-local and global primary surjections.
• Add hp, hq, GlobalQ, and aq annotations.
• Score paired surjections using alignment and fragment-length evidence.
• Learn the fragment-length distribution before enabling parallel processing.
• Handle paired reads with no compatible path/strand placement.
• Prune mapper-declared tail-region anchors.
• Classify overlapping alternatives as secondary and disjoint alternatives as
  supplementary.
• Rank paths by their best alignment.
• Realign anchors containing substitutions.
• Validate anchor sliding on target paths.
• Preserve tail annotations through GAF.
• Add implementation and evaluation documentation.

---

• Dependencies

This PR depends on two open libvgio PRs:

• Add grouped parallel I/O for graph placements of unpaired reads libvgio#93 adds grouped parallel input for alternative placements of unpaired reads.
• Add grouped parallel I/O for graph placements of interleaved reads libvgio#95 adds grouped parallel input for alternative placements of interleaved paired-end reads.
• Existing smaller surjection PRs may be superseded by this consolidated PR.

[adamnovak]

I've looked at the overall design here and the backing abstractions, and made some comments about them. I found some things about how this is written out that I would want to change (I don't think we want to merge the big Markdown reports, for example, and some things might want renaming).

I haven't really dug much into the core algorithms; I think they really need to be constrained by tests. If they give the right answers for weird edge cases I will be more confident in the core algorithms.

[adamnovak]

This report and diploid-map-surjection.md are both very thorough, and could help people understand this PR, but I don't think they make a lot of sense to carry in the vg repository, at least without organizing things so it's clear that these are reports about a particular project that was done at a particular time.

Some of the material in these would be useful on the vg wiki, because it covers how vg suject's diploid mode now works, and what the new tags mean.

But a lot of this is about how the changeset in particular was tested, and what particular changes to surjection are part of this project, and there are a lot of particular accuracy statistics and references to particular experiments you or I did at some point. That is all of limited utility to future readers of the vg repository, who will only care about how everything works at that future point, and would be better placed in a manuscript about your project.

[adamnovak]

--into-ref is described as just surject into this reference assembly, so we don't save a lot of space by referring to it, whereas it makes the user need to find another option later in the list to understand this one.

It's not quite clear from this if "haplotype-aware diploid surjection" is the name of what this turns on, or some downstream process that consumes the tags, or some other thing that can be flagged on separately and, if turned on, will emit those tags when this flag is used.

[adamnovak]

This then gets passed to e.g. vg::io::gam_paired_grouped_for_each_parallel_after_wait, which doesn't know about/depend on PairedSurjector::fragment_group_t and always wants to give you a vector of pairs of Alignments. That makes it seem like there might be some confusion around exactly what piece of code owns the idea of a fragment_group_t, in case we wanted to someday change it to a real struct or something.

I would describe the things that need to attach to libvgio in the same terms that libvgio does, and if that's too complicated and we need some abstractions for the types, the definitions should live in libvgio if we can get away with it.

[adamnovak]

What if someone provides --diploid-map with GAMP input? It doesn't look like we document that it won't work, or like we produce an error when asked to do that, but it also doesn't look like it actually activates diploid mode. So a user can think that they've run diploid surjection when they haven't.

[adamnovak]

I think most of this (maybe except for the ensure_alignment_is_for_graph() and/or set_metadata() calls) might want to be a method on Surjector (maybe an overload or variant of surject()?) because you need all these pieces to accomplish the operation of taking a collection of graph placements and producing the corresponding collection of linear placements. This is a piece we might want to be able to write a unit test for.

[adamnovak]

I'm wondering if there needs to be more to FragmentLengthDistribution so that this doesn't look so much like

vg/src/minimizer_mapper.hpp

Lines 538 to 549 in 3a37aa3

	bool fragment_distr_is_finalized () {return fragment_length_distr.is_finalized();}
	void finalize_fragment_length_distr() {
	if (!fragment_length_distr.is_finalized()) {
	fragment_length_distr.force_parameters(fragment_length_distr.mean(), fragment_length_distr.std_dev());
	}
	}
	void force_fragment_length_distr(double mean, double stdev) {
	fragment_length_distr.force_parameters(mean, stdev);
	}
	double get_fragment_length_mean() const { return fragment_length_distr.mean(); }
	double get_fragment_length_stdev() const {return fragment_length_distr.std_dev(); }
	size_t get_fragment_length_sample_size() const { return fragment_length_distr.curr_sample_size(); }

Why do we need to both have the FragmentLengthDistribution around and also wrap and expose all these methods? Does the abandonment feature need to be pushed down into FragmentLengthDistribution?

Maybe we should be inheriting from PairedEndMapper (which might need to be renamed to encompass other kinds of paired-end read processing)?

vg/src/mapper.hpp

Lines 144 to 168 in 3a37aa3

	/**
	* A class trait/mixin which defines a mapper's paired end distribution support.
	*
	* Doesn't actually define the paired-end mapping methods.
	*/
	class PairedEndMapper {
	public:
	/// Set parameters for estimating fragment length distribution.
	/// TODO: setting alignment threads after this could mess up the internal memory for how many threads to reset to
	void set_fragment_length_distr_params(size_t maximum_sample_size = 1000, size_t reestimation_frequency = 1000,
	double robust_estimation_fraction = 0.95);
	/// Returns true if fragment length distribution has been fixed
	bool has_fixed_fragment_length_distr();
	/// Use the given fragment length distribution parameters instead of
	/// estimating them.
	void force_fragment_length_distr(double mean, double stddev);
	protected:
	/// Holds the actual fragment length distribution and estimation information
	FragmentLengthDistribution fragment_length_distr;
	};

[adamnovak]

These comments maybe want to be in the header?

[adamnovak]

Should we accept input from other mappers with MAPQs higher than 60?

[adamnovak]

If this is some kind of scratch that we use through the whole function, it might make sense to explain how we make sure that it has the values we need it to have at any given point, or what it's supposed to look like after every pass through it.

[adamnovak]

I don't really understand this. I think it was like this before, but don't we need to figure out which alignment of this read the supplementary alignment supplements, and then figure out which alignment of the pair partner is actually paired with it?

Otherwise we might have a pair of primary non-supplementary alignments at the start of a chromosome, and a pair of secondary non-supplementary alignments at the end, and then all the supplementary alignments saying they belong to the pair at the start of the chromosome even if they themselves appear at the end, right?