BandedGlobalAligner gets stuck in traceback, but only when doing patch_alignment on particular PacBio reads

[adamnovak]

After I fixed the bug in #662 I found this one.

If you grab Shilpa's graph from that issue and align a particular read (which I have extracted here), the BandedGlobalAligner will get stuck:

vg map true-yeast-graph-x.vg -x true-yeast-graph-x.xg -g true-yeast-graph-x.gcsa -s `cat seq.txt` -t 1 > test.gam

error:[BandedGlobalAligner] traceback stuck at node boundary
vg: src/banded_global_aligner.cpp:1349: void vg::BandedGlobalAligner<IntType>::BAMatrix::traceback_internal(vg::BandedGlobalAligner<IntType>::BABuilder&, vg::BandedGlobalAligner<IntType>::AltTracebackStack&, int64_t, int64_t, vg::BandedGlobalAligner<IntType>::matrix_t, bool, int8_t*, int8_t*, int8_t, int8_t, bool, IntType) [with IntType = short int; int64_t = long int; int8_t = signed char]: Assertion `0' failed.
Aborted (core dumped)

It's happening inside patch_alignment in the Mapper.

You can get the files to reproduce this here: http://hgwdev.soe.ucsc.edu/~anovak/halIssues/banded_problem.zip

@jeizenga, you might be the right person to look at this.

[jeizenga]

@adamnovak @ekg

I think I've tracked down the problem, more or less. It's falsely identifying reversed nodes as source nodes in the DP. I think we need to think out some semantics before I fix it though.

The graph it's getting is a DAG in the sense that there are no directed or bidirected cycles. Accordingly, I could relatively easily induce an orientation on all of the nodes with a traversal. However, that introduces the ambiguity that the sources and sinks could just as easily be the sinks and sources. If we don't break symmetry somehow, I think any algorithm will sometimes try to align the read in the opposite strand to what was intended.

Up until now, I've been assuming that this would all take place in the subgraph preprocessing so that I would never see a reversing edge. That's still an option, although we'll need some auxilliary structure to keep track of which nodes have been reversed so that we can label them as such in the alignment's Mappings. Another option would be to have the function call require one NodeTraversal that's in the same orientation as the read. Then I could extend out the induced orientation to the rest of the DAGified subgraph.

Anyone have preferences?

[adamnovak]

I know in most cases for short reads we're throwing it at the alignment machinery in both orientations and seeing what sticks.

For this alignment-patching-up context things are more complicated, because we need to get alignments in "the same" orientation to a bunch of disjoint subgraphs so that we can string them together.

I'm not sure that that's actually technically possible, because there's no way we can know by looking at the subgraphs what orientations you can get between them in.

I'm not really sure how to resolve that, but I can say that testing on Cactus graphs for things like this is useful, because they will run the primary path forward through some graph regions and backward through others.

It might be smartest to move removal of all reversing edges up to dagify or unfold and make worrying about consistency when pasting together alignments the problem of the code doing the pasting. If an arbitrary decision about which orientation is "forward" gets made inside the aligner, that's bad.

On Wed, Feb 22, 2017 at 3:17 PM, Jordan Eizenga notifications@github.com wrote:

@adamnovak https://github.com/adamnovak @ekg https://github.com/ekg

I think I've tracked down the problem, more or less. It's falsely identifying reversed nodes as source nodes in the DP. I think we need to think out some semantics before I fix it though.

The graph it's getting is a DAG in the sense that there are no directed or bidirected cycles. Accordingly, I could relatively easily induce an orientation on all of the nodes with a traversal. However, that introduces the ambiguity that the sources and sinks could just as easily be the sinks and sources. If we don't break symmetry somehow, I think any algorithm will sometimes try to align the read in the opposite strand to what was intended.

Up until now, I've been assuming that this would all take place in the subgraph preprocessing so that I would never see a reversing edge. That's still an option, although we'll need some auxilliary structure to keep track of which nodes have been reversed so that we can label them as such in the alignment's Mappings. Another option would be to have the function call require one NodeTraversal that's in the same orientation as the read. Then I could extend out the induced orientation to the rest of the DAGified subgraph.

Anyone have preferences?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/678#issuecomment-281837277, or mute the thread https://github.com/notifications/unsubscribe-auth/AE0_Xyl-oIi0INnmM3x0Mjk8ABIoanXOks5rfMIUgaJpZM4MIDPH .

[jeizenga]

Agree about the arbitrary decision.

We do have the "cut points" for the subgraphs at the end of the MEMs. Those could be useful. Theoretically, they should induce an orientation, although I can think of some weird edge cases. I think it would play better with the newer subgraph extraction algorithm I implemented.

In any case, I think the two main options are:

• Pass a cut point to the global aligner and let it orient all of the nodes with a traversal
• Externally, orient the graph before passing it to the global aligner so that there are no reversing edges, record the translation, and translate the alignments back after the aligner returns them.

[ekg]

Presumably this can be resolved in the external handling. It shouldn't be that the aligner is getting something with reversing edges because we should be unfolding the graph. Maybe the check for this is failing. I'll look into this. Could be behind the problems with MSGA.

[ekg]

If I dump the graph out before alignment that we're aligning against I can't then reproduce this problem using vg align. The corresponding graph doesn't have any reversing edges either.

Do either of you have the graph the banded aligner is failing on handy?

[adamnovak]

I also couldn't get this to fail on the dumped graph, which is why I submitted the bigger test case.

Does serialization ever re-order nodes?

On Thu, Feb 23, 2017 at 8:35 AM, Erik Garrison notifications@github.com wrote:

If I dump the graph out before alignment that we're aligning against I can't then reproduce this problem using vg align. The corresponding graph doesn't have any reversing edges either.

Do either of you have the graph the banded aligner is failing on handy?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/678#issuecomment-282045494, or mute the thread https://github.com/notifications/unsubscribe-auth/AE0_X09-iUjJLsI0kWKMib4SdtwQHcwsks5rfbU0gaJpZM4MIDPH .

[jeizenga]

So, modification to my previous comment. This was actually two bugs: one in my DP and one in the subgraph extraction. I've fixed the DP bug and a PR will be coming soon.

As for the subgraph extraction: I looked into it a little more, and the reversing edge I'm talking about doesn't seem to be included in the Graph that's passed to the aligner at all. Check out this viz:

The graph that the global aligner gets is missing 187054 and 187055, even though it includes 187053 and 102495. This is why I was getting 187053 falsely labelled as a tip--because it actually is one in the subgraph.

Another example:

Here I'm missing 187065 and 187066 even though 102473 and 187064 are included, which causes 187064 to be falsely labeled as another tip.

Between the two examples, I think the clear culprit is the reversing edge. One option we have is to replace the subgraph extraction algorithm with the new one I implemented for the multipath alignments, which provides more algorithmic guarantees. Thoughts?

[ekg]

@jeizenga are you getting these using the serialization in #ifdef debug stanza at the beginning of the do_align local function in vg::VG::align? Using that serialization, I can't get any reversing edges in the graph that's passed into the aligner.

I have a commit that I'm merging in that resolves this particular bug by making the alignment invocation match the one in vg align. However, I don't know if it matches the spirit of how we should be running the banded aligner, so please have a look at it.

[jeizenga]

@ekg I think I haven't made myself clear. It's not that there are reversing edges are in the subgraph. Rather, it's that there are nodes/edges in the parent VG graph that should be in the graph but aren't. The commonality seems to be that these nodes are connected by reversing edges in the parent graph (see the pictures above). When those nodes get dropped, it creates false source nodes, which this global aligner is perfectly happy to use (since it does a source-to-sink alignment). Essentially, it thinks those nodes are directly adjacent to the MEM even though they are actually quite distant in the graph. That means that, after patching, the alignment will take an "edge" that doesn't actually exist.

[ekg]

OK that does clarify things some. What aspect of this is my fix exploiting? The source node selection?

On Fri, Feb 24, 2017, 7:38 PM Jordan Eizenga notifications@github.com wrote:

@ekg https://github.com/ekg I think I haven't made myself clear. It's not that there are reversing edges are in the subgraph. Rather, it's that there are nodes/edges in the parent VG graph that should be in the graph but aren't. The commonality seems to be that these nodes are connected by reversing edges in the parent graph (see the pictures above). When those nodes get dropped, it creates false source nodes, which this global aligner is perfectly happy to use (since it does a source-to-sink alignment). Essentially, it thinks those nodes are directly adjacent to the MEM even though they are actually quite distant in the graph. That means that, after patching, the alignment will take an "edge" that doesn't actually exist.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/678#issuecomment-282369038, or mute the thread https://github.com/notifications/unsubscribe-auth/AAI4EVm3xDEIMoQ2AVfRXNCVgwQKwQWQks5rfyOKgaJpZM4MIDPH .

[jeizenga]

Are we sure that it's not taking the false "edges"? It might have just wiggled things so that the alignment doesn't run into the DP bug. The false source node bug won't make it explode, just have incorrect alignments.