Personalization features for non-autosomal chromosomes

[eizengaj-roche]

Changelog Entry

To be copied to the draft changelog by merger:

• vg haplotypes adds a --wrap option for circular chromosomes, which doubles selected haplotypes after personalization
• vg gbwt adds a --wrap-contig option that performs the same operation on an existing GBWT/GBZ
• vg haplotypes adds new scoring modes for non-diploid chromosomes: --high-cov-contig/--high-cov-num-haps, --half-cov-contig/--half-cov-num-haps, and --exclude-contig

Description

These options are intended to better support non-autosome chromosomes, which need somewhat different treatment. Previously, the best practices personalization pipeline applied a uniform, diploid scoring model to all chromosomes, which is normed to the expected coverage of a unique, diploid region of the genome. This led to a number of issues in a human pangenome:

• chrM has much higher coverage than the autosomes (the factor varies by tissue, but in the range of 10s-1000s fold higher), leading the actual informative k-mers to be identified as frequent and treated as uninformative
• chrX and chrY have lower coverage than the autosomes in XY karyotypes, leading the actual informative k-mers to be identified as heterozygous and treated as uninformative during the initial phase of haplotype selection
• The second pass of diploid scoring does not make sense for chrM or for chrX/chrY in XY karyotypes

There are also further complications for mapping to chrM due to the fact that it is a circular chromosome. GBZ does not natively support circular paths, in the sense of the libhandlegraph interface for them. The paths can be made "pseudo-circular" by wrapping the haplotypes around the path a second time. However, in that case, vg haplotypes refuses to operate on the graph because it contains a top-level chain that is a cycle. You might think that the vg giraffe --supplementary option might help obtain (at least some) alignments that wrap the 0-coordinate of chrM. Unfortunately, this doesn't help, because of the prevalence of NuMTs that include the 0-coordinate: there is always a contiguous alignment available for an entire short read, so the supplementary code path never has a reason to execute.

This PR includes two main contributions to address these issues. The first is an addition of two new scoring models that can be selected on a per-contig basis:

Model	Standard	Half coverage	High coverage
Absent score	-a	-a	-a
Heterozygous score	0	1	-a
Present score	1	0	-a
Frequent score	0	0	1

The relevant contigs are specified by new command line options in vg haplotypes like this:

vg haplotypes --high-cov-contig chrM --half-cov-contig chrX --half-cov-contig chrY

In addition, there are accompanying parameters that indicate the number of haplotypes sampled under these models: --high-cov-num-haps and --half-cov-num-haps. Any contigs specified with these scoring models are automatically excluded from the second phase of diploid sampling, if diploid sampling is specified. Finally, there is also an option --exclude-contig, which can exempt a contig from personalization entirely.

The next contribution is meant to address the issues with chrM being circular. Since having truly circular contigs is unsupported in GBZ, I went the less-principled-but-sufficient route of doubling the haplotype paths. This functionality can be accessed in vg gbwt during the modification step via the option --wrap-contig <CONTIG>. I also created a second access point as a mid-execution modification in vg haplotypes with the option --wrap <CONTIG>. The motivation for the second access point is that it allows you to apply the modification after haplotype sampling (which requires a linear top-level chain) but while still working with a mutable GBWTBuilder. In contrast, the vg gbwt entrypoint has to first convert back from the static GBWT to the GBWTBuilder to wrap a contig. The wrapping step is only applied to haplotype-sense paths out of fear of messing up coordinate systems by wrapping reference- or generic-sense paths.

The PR's implementation does have some weaknesses that I'm aware of:

• The by-contig interface for scoring models doesn't account for something like the PAR. Since the HPRC pangenomes follow the convention of hard-masking the PAR, it's probably always better to have diploid sampling in the PAR, even in XY karyotypes, but the current interface wouldn't allow that. I'm open to suggestions on improving this aspect.
• The fact that the XX/XY karyotype has to be known ahead of time is a bit annoying to me, since it should be possible to tell what karyotype is fairly easily from the KFF. I couldn't think of a strategy for doing so that felt suitably general to warrant being included in VG though.
• Wrapping a haplotype can create incorrect adjacencies if a haplotype is clipped at the end. At the beginning, this can be detected by a non-zero start coordinate, but at the end we would need to know something about the "theoretical length" of the path as well. In chrM this isn't really a problem because it's so short that it's hard not to assemble T2T, but this could potentially be an issue in other circular chromosomes.

To be transparent, I've done more benchmarking on chrM than on chrX/chrY. I can pretty definitively say that you get better results with --wrap chrM --high-cov-contig chrM. If you also use --supplementary in vg giraffe or vg surject, you can get the 0-coordinate spanning reads nicely represented in the output BAM as well. For the allosomes, I think the logic is sound, but I need to do more benchmarking to be confident that it's a net positive in practice.

[jltsiren]

A few comments:

• You currently wrap paths if they are haplotypes and start at 0. If the haplotype is fragmented, the first fragment then wraps back to the beginning, which probably doesn't make sense. You could use gbwt::FragmentMap to determine if a haplotype is fragmented, and then either refuse to wrap it or append the first fragment to the last.

• resolve_chains_by_name() uses nested loops for mapping contig names to top-level chain ids. You could just check against the contig names stored in TopLevelChain objects, as those are supposed to be the canonical contig names for the chains.

• wrap_haplotype_paths() checks the given contig names against the actual contig names for each path. According to HPRC conventions, the contig names for haplotype paths are accession numbers. If a haplotype is fragmented into several assembly contigs, each fragment has a distinct contig name and starts at 0.

• Additionally, wrap_haplotype_paths() takes bidirectionality of the GBWT into account when determining sequence_id but always inserts the path in both orientations.

• Recombinator::copy_chain() considers each copied path a reference path in the statistics. Is this intended behavior?

In the long term, the .hapl file should know what kind of a sequence each top-level chain corresponds to. Sampling could then choose the appropriate behavior automatically, instead of having the user specify it by contig.

[eizengaj-roche]

Thanks for the close reading! I think I have a better understanding of some of the data structures now. My responses:

You currently wrap paths if they are haplotypes and start at 0. If the haplotype is fragmented, the first fragment then wraps back to the beginning, which probably doesn't make sense. You could use gbwt::FragmentMap to determine if a haplotype is fragmented, and then either refuse to wrap it or append the first fragment to the last.

I'm now using the next/prev functions in the FragmentMap to locate the first and last fragment. The check for 0 is still there as a guard.

resolve_chains_by_name() uses nested loops for mapping contig names to top-level chain ids. You could just check against the contig names stored in TopLevelChain objects, as those are supposed to be the canonical contig names for the chains.

I am doing it this way now.

wrap_haplotype_paths() checks the given contig names against the actual contig names for each path. According to HPRC conventions, the contig names for haplotype paths are accession numbers. If a haplotype is fragmented into several assembly contigs, each fragment has a distinct contig name and starts at 0.

Am I understanding correctly that the concern is that each fragment has a separate contig name, and could therefore be wrapped onto itself, despite actually being a mid-chromosome contig? If so, that does seem like a sticky issue. I guess we would need some strategy to associate a first and last contig for each sample on a given chain. Does machinery for that exist? I don't think I fully grok the suffix completion code in the Recombinator, but maybe that would be useful? That wouldn't automatically help with the vg gbwt entrypoint though.

Additionally, wrap_haplotype_paths() takes bidirectionality of the GBWT into account when determining sequence_id but always inserts the path in both orientations.

Good catch. The inserts now respect the bidirectionality of the GBWT.

Recombinator::copy_chain() considers each copied path a reference path in the statistics. Is this intended behavior?

No, it was not. I'm now tracking excluded contigs separately.

Regarding a change in the .hapl format, I like the idea, but I think it's probably better for me to leave that as future work for someone who understands the implications of changing the format better.

Anyway, this commit may or may not be mergeable depending on what we think the appropriate strategy is for your 3rd comment.

[jltsiren]

Am I understanding correctly that the concern is that each fragment has a separate contig name, and could therefore be wrapped onto itself, despite actually being a mid-chromosome contig? If so, that does seem like a sticky issue. I guess we would need some strategy to associate a first and last contig for each sample on a given chain. Does machinery for that exist? I don't think I fully grok the suffix completion code in the Recombinator, but maybe that would be useful? That wouldn't automatically help with the vg gbwt entrypoint though.

The way you are now handling this in Recombinator looks correct.

In the general case, there is a risk of wrapping a mid-chromosome contig onto itself. We had a proposal for haplotype-level metadata in GFA, which could be used for mapping assembly contigs to chromosomes and listing them in order. But the proposal has not been merged, because nobody has had the time to implement it.

[eizengaj-roche]

I wonder if the vg gbwt interface is premature then. For my own purposes, the personalization entry point is more important, but if that were the only one, it would preclude using the feature on, e.g., the clipped graphs. Another option would be to leave the vg gbwt interface intact, but to try to warn users if there's a mid-chromosome loop. I think the way I would check that would be to check if the haplotype walk's endpoints near the endpoints of reference-sense path.

Do you have any leaning one way or the other?

[jltsiren]

I don't have any preferences on that. The vg gbwt subcommand is kind of dangerous anyway, as it lets you do many things unconditionally without checking if it would make sense.

[eizengaj-roche]

Makes sense. My read is that this is mergeable then. I'll do that in a day or two if nobody objects.