Ambiguous path production

[bricoletc]

I produced nested PRGs (--min_match_length=7, --max_nesting=5) from multiple sequence alignments of plasmodium genes and found that make_prg can produce ambiguous PRGs.

Here's an example I drew out by hand to understand

_o and _c mean site opening and closing nodes, with a numbered ID before. 7T means 7 consecutive Ts.

The path CA9T can be produced going through either sites 33/34/35 or 33/36. This can be horrible for genotyping as if the dataset has that sequence then you have to either say you can't genotype because of ambiguity, or randomly choose one of the two possibilities. The latter could cause comparison issues across datasets. For now in gramtools I've decided to bail out in these cases and null genotype, with a FILTER flag to AMBIG to signal this occurred.

[iqbal-lab]

Is it possible to define a total order on alternate paths that spell the same sequence? God help us

[bricoletc]

Here's another example which occurs when building PRG from m.tuberculosis, sequences from 1,017 samples, variation called with cortex in genomic positions: NC_000962.3:336600-337630 and this region starts at POS 337368

157.pdf

The following sequence: GCTCCGCCGGTCCCGCCGGTCC can be spelled by following two paths:

• 157/1, 160/1, 161/0
• 158/0, 159/0

[bricoletc]

At some level the problem seems to boil to excessive clustering, ie putting sequences that are very similar in different clusters.

Another example where this has clearly occurred:
209.pdf
From Pfalciparum gene DBLMSP, 2,500 input sequences

This should potentially be a single variant site enumerating all the similar sequences.

[iqbal-lab]

These must always be horizontally offset, right?
Is it possible for this to happen with just one leve of nesting?
I guess the answer is yes, you can delete an A and then add two As, or add an A and then add an A

[iqbal-lab]

is it possible without indels? i guess not?

[bricoletc]

The kmer counting underlying the clustering ignores alignment gaps- if that's what you're referring to?

Ie, gaps don't change anything: AAA-AAA and AA-AAAA are the same for clustering (at the moment)

[iqbal-lab]

Doh, of course

[iqbal-lab]

Rachel says couyld you please put the MSAs that generated these examples in?

[bricoletc]

Sure, i will need to extract them, and reproduce the graphs- the images are subgraphs from larger PRGs.

I will do this later this week, and also think we can add tests and some work to start alleviating this, if we can clearly define what's desirable and undesirable in these situations

[iqbal-lab]

I will propose that the only thing you need to avoid is having two different source/sink paths which spell the same sequence.

[leoisl]

I took a closer look at the first example, and I think this is hard to solve for all cases. We might make some improvements here and there, based on some edge cases we find by hand, and it might work for these 3 examples, but this is just a very small number of cases, we won't be sure if it improved overall. I think this is due to a misclustering of sequences, but this is unavoidable as the k-means clustering process is very heuristic, and it might undercluster or overcluster (it is really hard to set the correct inertia parameter for every cluster, if there is one). So, I think that PRGs having ambiguous paths will unavoidably happen. Having this heuristic component also makes it hard to have any algorithm that will work on all situations (I think that whatever algorithm we think about, we can't have a proof of correctness, and we will always be able to develop an edge case in such a way that a misclustering of sequence happens, and this induces ambiguous paths in the PRG).

So, we could think of post-processing algorithms to remove some edges/nodes from the PRG in order to keep a sort of path-spelling uniqueness property. That is not possible neither, I am afraid. In the first graph, it is clear we can remove the CA9T node from the PRG, as it is a redundant node: if we remove it, we can still generate all sequences that we could generate before, and without duplication. So we lose no information, and now we have path-spelling uniqueness. However, it is easy to modify this example in such a way that the CA9T node can't be removed without changing the set of paths that the PRG spells (i.e. add a sequence that can only be generated if going through the CA9T node), and this is a counter-example showing that if a PRG has ambiguous paths, we can't remove nodes or edges in order to break this ambiguity (I don't think breaking ambiguity on the cost of losing paths is worth it).

It seems to me that in some edge cases, we can have path ambiguity in the PRG, and it will be responsibility of the downstream algorithms to be aware and take care of this.

Of course I could be wrong on any of the statements above, and as in bioinformatics, maybe we don't need a correct algorithm, but one that works on our cases (and hope they work on other)...

[leoisl]

PS: I am not advocating we replace k-means clustering. I think it makes total sense in make_prg, and we should use it as we are using now

[bricoletc]

Thanks for pitching in @leoisl , I have two thoughts:
i) whatever the upstream construction algorithm, downstream application needs some kind of awareness of path ambiguity if we know it can happen. in gramtools if i detect it at genotyping, I null all the calls in the path and set a "AMBIG" filter
ii) i'm not sure we can rule out being able to guarantee that all paths in produced graph spell different sequences. even if we cannot make that guarantee, i think the current clustering can be improved so as to reduce the average similarity of sequences in different branches

[iqbal-lab]

i think it's an interaction of

• the MSA (can we rely on the MSA algorithm to align the same sequence the same way on two different lines?)
• the kmer-size for collapse k1
• the kmer size for clustering k2
• indels and local small scale indels shorter than k2.
• nesting
  So i wonder if
  a) it might be possible to detect that there is >=1 indel and local sequence identity and decide to stop nesting further.
  b) i think i need to see a bunch of examples
  for your initial example @bricoletc i think there must be some error. you say

The path CA9T can be produced going through either sites 33/34/35 or 33/36.

I don't see how you can get CAT from 33/34/35, you can get CATTTC