Repository navigation
surject changes: detect target-path repeats during anchor pruning - #5052
Sagorikanag wants to merge 4 commits into
Conversation
e03619f to
9fc5cd8
Compare
adamnovak
left a comment
adamnovak
left a comment
There was a problem hiding this comment.
This looks like it's going to work.
The tests are too interdependent for me to be able to look at any of them and tell whether they are asserting the right things or whether they are constructing the situation they claim to be constructing.
Also, I think the implementation shouldn't be based around pulling substring copies out of the path for each slide value; we should pull one string window and do all our string checks inside that, kind of like we do them inside the read already.
| // Preserve the read-side search radius, including for anchors with indels. | ||
| const size_t read_span = chunk.first.second - chunk.first.first; | ||
| const size_t read_slide_limit = min<size_t>(max_slide, read_span * 2); | ||
| const size_t read_start = chunk.first.first - sequence.begin(); | ||
| const size_t read_remaining = sequence.end() - chunk.first.second; | ||
| for (int64_t distance = -static_cast<int64_t>(read_slide_limit); | ||
| distance <= static_cast<int64_t>(read_slide_limit); ++distance) { |
There was a problem hiding this comment.
Isn't this implementing a search radius in read space, rather than "preserve"ing "the" search in read space? This is itself the search in read space.
And where's the indel handling? The comment suggests we need to be doing something special here to handle indels in anchors, but I can't see that any such special thing is being done, which makes the code look wrong.
| const auto first_step = step_ranges[i].first; | ||
| const auto first_handle = graph->get_handle_of_step(first_step); | ||
| const auto path_handle = graph->get_path_handle_of_step(first_step); | ||
| const auto& first_pos = chunk.second.mapping(0).position(); | ||
| const bool reverse_on_path = | ||
| first_pos.is_reverse() != graph->get_is_reverse(first_handle); | ||
|
|
||
| const size_t step_start = graph->get_position_of_step(first_step); | ||
| const size_t node_length = graph->get_length(first_handle); | ||
| const size_t path_length = graph->get_path_length(path_handle); | ||
| const size_t mapping_offset = first_pos.offset(); |
There was a problem hiding this comment.
This is starting to read like Javascript, where every project has a linter that will harass you to make anything that never changes const.
But in JS you can just say const to declare a variable, and in C++ it's both more typing to make all the locals cosnt and not, as far as I know, usual practice. See for example https://quuxplusone.github.io/blog/2022/01/23/dont-const-all-the-things/#local-variables-rarely-const (but see dissenting opinions at https://stackoverflow.com/a/38926475).
@Sagorikanag Can you explain why you think marking locals as const when possible is in general better? Is this helping us keep track of what's going on in this giant function? Or are you doing this for no particular reason (in which case I would drop all the consts on locals that aren't important)?
| if (ref_span > path_length - slid_start) { | ||
| continue; | ||
| } | ||
| const string slid_ref_sequence = path_window(slid_start); |
There was a problem hiding this comment.
This is going to pull almost the same sequence on each successive call, but it's going to build it from scratch every time, by scanning almost the same nodes of the path. I think that's unnecessarily inefficient.
I guess you can say the runtime is already O(path_slide_limit * ref_span) in the worst case, because we need equality tests on all those bases. But when mostly we can stop a lot of those equality tests early, the expected runtime is much faster, whereas here we're definitely spending O(path_slide_limit * ref_span) in even the best case making all those copies of substrings of the path. So I feel like it would be a lot better to traverse the graph path once, and extract the entire relevant window, and then do a bunch of std::equal comparisons within it. We could even maybe find a way to pull out a check-for-slide function and use it in both the read and reference cases, once we've pulled the strings.
(Actually, it would probably be even better to pull the window, and use std::find and std::rfind with appropriate bounds, and thus figure out whether there's an instance of the string other than the center one, because I bet they are faster than repeated std::equal. But we might not need to do that here.)
| SECTION("Read insertion preserves the read-based search radius") { | ||
| region = "ACGATTCGCCCC"; | ||
| insertion = true; | ||
| read_between = "TTACGCTAGA"; | ||
| max_slide = 12; | ||
| } |
There was a problem hiding this comment.
I think the whole set of test cases here is definitely too much of a test construction kit. If there's a good reason to have a test case construction kit, we could do it, but we'd need to explain it a lot better. We have an insertion flag that isn't documented at the top of the function, which we set here, which has the behavior of inserting "GCTA" at offset 2 in the anchor sequence, which we pulled from part of the path sequence at offset 0, which we built by concatenating a sentinel sequence onto the end of a region sequence. If this test later fails at REQUIRE(chunks.size() == (prune ? 1 : 2)); it's going to be a project to reconstruct the actual test case and what's actually supposed to be true about the result that isn't.
| // A distinct longer anchor prevents the keep-one-anchor fallback from | ||
| // hiding removal of the anchor being tested. | ||
| const string sentinel = "GATCCGTAGTCACTGACCTAGGTC"; | ||
| const string path_sequence = sentinel_first | ||
| ? sentinel + region : region + sentinel; |
There was a problem hiding this comment.
Would it make more sense to add a flag to control that behavior and turn it off, so we can more easily get at the part we're trying to test?
We could also refactor things so the part that identifies which anchors ought to be kept according to the shift test is its own piece that we can call and test separately. Then those results would later be processed to actually remove the anchors, subject to the constraint that at least one anchor remains.
There was a problem hiding this comment.
I went with your refactoring suggestion. anchor_has_nearby_repeat() now checks an anchor without modifying it, so the tests can exercise the sliding logic directly without a flag.
prune_and_trim_anchors() uses that result to mark anchors for removal, while retaining the existing keep-at-least-one-anchor behavior.
9fc5cd8 to
a0c9e33
Compare
a0c9e33 to
54be5a8
Compare
This branch has not been deployed
2 participants
Changelog Entry
vg surjectnow checks for nearby repeated anchor sequence on the target path, in addition to the existing read-side repeat check.Description
This change extends suspicious-anchor pruning so an anchor can be marked for pruning when its reference sequence occurs again nearby on the target path.
The existing read-side repeat search is preserved with its original read-based search radius. The new target-path search uses the anchor's reference span and computes orientation relative to the target path, including reverse-oriented path steps and multi-node anchors.
Tests cover read-only and path-only repeats, slide boundaries, path boundaries, forward/reverse orientation combinations, multi-node anchors, reverse-oriented nodes, and differing read/reference spans caused by insertions.
Stack created with GitHub Stacks CLI • Give Feedback 💬