From 8ceb4cd678b36cc383675673e62cbd6bfb182733 Mon Sep 17 00:00:00 2001 From: Jordan Eizenga Date: Thu, 1 Aug 2024 13:37:48 -0700 Subject: [PATCH] bugfix invalid surjected alignments --- src/multipath_alignment.cpp | 8 ++++++++ src/multipath_alignment_graph.cpp | 11 ++++++----- src/surjector.cpp | 6 ++++++ 3 files changed, 20 insertions(+), 5 deletions(-) diff --git a/src/multipath_alignment.cpp b/src/multipath_alignment.cpp index d91ae71a03a..7d7679c15d2 100644 --- a/src/multipath_alignment.cpp +++ b/src/multipath_alignment.cpp @@ -3857,6 +3857,14 @@ namespace vg { #endif return false; } + for (const auto& connection : subpath.connection()) { + if (connection.next() <= i) { +#ifdef debug_verbose_validation + cerr << "validation failure on connection topological order" << endl; +#endif + return false; + } + } } } diff --git a/src/multipath_alignment_graph.cpp b/src/multipath_alignment_graph.cpp index ce36787a12c..a43a7a2722c 100644 --- a/src/multipath_alignment_graph.cpp +++ b/src/multipath_alignment_graph.cpp @@ -6170,17 +6170,17 @@ void MultipathAlignmentGraph::align(const Alignment& alignment, const HandleGrap if (aligning_tail_length < tail_length) { // Tail is too long. Just make a softclip directly in the base graph ID space. // TODO: What if we just don't produce this? Do we get softclips for free? + + size_t cut_point = (path_node.end - alignment.sequence().begin()) + aligning_tail_length; + size_t tail_remaining = tail_length - aligning_tail_length; #ifdef debug_multipath_alignment cerr << "softclip long right" << endl; #endif - size_t cut_point = (path_node.end - alignment.sequence().begin()) + aligning_tail_length; - size_t tail_remaining = tail_length - aligning_tail_length; for (auto& alt_aln : alt_alignments) { - alt_aln.mutable_sequence()->append(alignment.sequence(), cut_point, tail_remaining); if (!alt_aln.quality().empty()) { - alt_aln.mutable_sequence()->append(alignment.quality(), cut_point, tail_remaining); + alt_aln.mutable_quality()->append(alignment.quality(), cut_point, tail_remaining); } Mapping* mapping = alt_aln.mutable_path()->mutable_mapping(alt_aln.path().mapping_size() - 1); @@ -6188,6 +6188,7 @@ void MultipathAlignmentGraph::align(const Alignment& alignment, const HandleGrap if (final_edit->from_length() == 0 && !final_edit->sequence().empty()) { // extend the softclip final_edit->set_sequence(final_edit->sequence() + alignment.sequence().substr(cut_point, tail_remaining)); + final_edit->set_to_length(final_edit->sequence().size()); } else { final_edit = mapping->add_edit(); @@ -6342,7 +6343,7 @@ void MultipathAlignmentGraph::align(const Alignment& alignment, const HandleGrap alt_aln.mutable_sequence()->append(alignment.sequence(), 0, tail_remaining); if (!alt_aln.quality().empty()) { - alt_aln.mutable_sequence()->append(alignment.quality(), 0, tail_remaining); + alt_aln.mutable_quality()->append(alignment.quality(), 0, tail_remaining); } Mapping* mapping = alt_aln.mutable_path()->mutable_mapping(0); diff --git a/src/surjector.cpp b/src/surjector.cpp index b428397223a..e7a59e2615e 100644 --- a/src/surjector.cpp +++ b/src/surjector.cpp @@ -4240,7 +4240,13 @@ using namespace std; auto mappings = path_chunk.second.mutable_mapping(); mappings->erase(mappings->begin(), mappings->begin() + mappings_to_delete); auto edits = mappings->front().mutable_edit(); + size_t deleting_from_length = 0; + for (size_t e = 0; e < edits_to_delete; ++e) { + deleting_from_length += (*edits)[e].from_length(); + } edits->erase(edits->begin(), edits->begin() + edits_to_delete); + auto position = mappings->front().mutable_position(); + position->set_offset(position->offset() + deleting_from_length); // trim ref interval for (size_t m = 0; m < mappings_to_delete; ++m) {