Merge pull request #3376 from vgteam/zip-mp-aln-tails

Merge common portions of multipath alignment tail alignments

src/aligner.cpp

@@ -1559,17 +1559,17 @@ int32_t Aligner::score_full_length_bonus(bool left_side, const Alignment& alignm
     return full_length_bonus;
 }
 
-int32_t Aligner::score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const Path& path,
+int32_t Aligner::score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const path_t& path,
                                          string::const_iterator seq_begin, bool no_read_end_scoring) const {
 
     int32_t score = 0;
     string::const_iterator read_pos = seq_begin;
     bool in_deletion = false;
     for (size_t i = 0; i < path.mapping_size(); i++) {
-        const Mapping& mapping = path.mapping(i);
+        const auto& mapping = path.mapping(i);
 
         for (size_t j = 0; j < mapping.edit_size(); j++) {
-            const Edit& edit = mapping.edit(j);
+            const auto& edit = mapping.edit(j);
 
             if (edit.from_length() > 0) {
                 if (edit.to_length() > 0) {
@@ -2255,7 +2255,7 @@ int32_t QualAdjAligner::score_full_length_bonus(bool left_side, const Alignment&
                                    alignment.quality().begin());
 }
 
-int32_t QualAdjAligner::score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const Path& path,
+int32_t QualAdjAligner::score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const path_t& path,
                                                 string::const_iterator seq_begin, bool no_read_end_scoring) const {
 
     int32_t score = 0;
@@ -2264,7 +2264,7 @@ int32_t QualAdjAligner::score_partial_alignment(const Alignment& alignment, cons
 
     bool in_deletion = false;
     for (size_t i = 0; i < path.mapping_size(); i++) {
-        const Mapping& mapping = path.mapping(i);
+        const auto& mapping = path.mapping(i);
 
         // get the sequence of this node on the proper strand
         string node_seq = graph.get_sequence(graph.get_handle(mapping.position().node_id(),
@@ -2273,7 +2273,7 @@ int32_t QualAdjAligner::score_partial_alignment(const Alignment& alignment, cons
         string::const_iterator ref_pos = node_seq.begin() + mapping.position().offset();
 
         for (size_t j = 0; j < mapping.edit_size(); j++) {
-            const Edit& edit = mapping.edit(j);
+            const auto& edit = mapping.edit(j);
 
             if (edit.from_length() > 0) {
                 if (edit.to_length() > 0) {

src/aligner.hpp

@@ -14,6 +14,7 @@
 #include "Variant.h"
 #include "Fasta.h"
 #include "handle.hpp"
+#include "path.hpp"
 #include "dozeu_interface.hpp"
 #include "deletion_aligner.hpp"
 
@@ -182,7 +183,7 @@ namespace vg {
         virtual int32_t score_full_length_bonus(bool left_side, const Alignment& alignment) const = 0;
 
         /// Compute the score of a path against the given range of subsequence with the given qualities.
-        virtual int32_t score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const Path& path,
+        virtual int32_t score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const path_t& path,
                                                 string::const_iterator seq_begin, bool no_read_end_scoring = false) const = 0;
 
         /// Returns the score of an insert or deletion of the given length
@@ -370,7 +371,7 @@ namespace vg {
 
         int32_t score_full_length_bonus(bool left_side, const Alignment& alignment) const;
 
-        int32_t score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const Path& path,
+        int32_t score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const path_t& path,
                                         string::const_iterator seq_begin, bool no_read_end_scoring = false) const;
 
     private:
@@ -429,7 +430,7 @@ namespace vg {
 
         int32_t score_full_length_bonus(bool left_side, const Alignment& alignment) const;
 
-        int32_t score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const Path& path,
+        int32_t score_partial_alignment(const Alignment& alignment, const HandleGraph& graph, const path_t& path,
                                         string::const_iterator seq_begin, bool no_read_end_scoring = false) const;
 
 

src/multipath_alignment_graph.hpp

@@ -305,6 +305,12 @@ namespace vg {
                     size_t max_alt_alns, bool dynamic_alt_alns, size_t max_gap, double pessimistic_tail_gap_multiplier,
                     size_t min_paths, unordered_set<size_t>* sources = nullptr);
 
+        /// If a list of aligned subsequences are identifical in a prefix/suffix, remove that
+        /// prefix/suffix from all of the alignments and return it as a separate alignment.
+        /// If there is no shared prefix/suffix, returns an empty path with 0 score.
+        static pair<path_t, int32_t> zip_alignments(vector<pair<path_t, int32_t>>& alt_alns, bool from_left,
+                                                    const Alignment& alignment, const HandleGraph& align_graph,
+                                                    string::const_iterator begin, const GSSWAligner* aligner);
 
         /// Memo for the transcendental pessimistic tail gap function (thread local to maintain thread-safety)
         static thread_local unordered_map<double, vector<int64_t>> pessimistic_tail_gap_memo;

src/multipath_mapper.cpp

@@ -54,7 +54,7 @@ namespace vg {
     void MultipathMapper::multipath_map_internal(const Alignment& alignment,
                                                  MappingQualityMethod mapq_method,
                                                  vector<multipath_alignment_t>& multipath_alns_out) {
-
+        
 #ifdef debug_multipath_mapper
         cerr << "multipath mapping read " << pb2json(alignment) << endl;
         cerr << "querying MEMs..." << endl;

src/position.cpp

@@ -75,4 +75,10 @@ string debug_string(const position_t& pos) {
     return to_return;
 }
 
+void from_proto_position(const Position& from, position_t& to) {
+    to.set_node_id(from.node_id());
+    to.set_offset(from.offset());
+    to.set_is_reverse(from.is_reverse());
+}
+
 }

src/position.hpp

@@ -61,7 +61,7 @@ pair<int64_t, int64_t> min_oriented_distances(const unordered_map<path_handle_t,
                                               const unordered_map<path_handle_t, vector<pair<size_t, bool> > >& path_offsets2);
 
 string debug_string(const position_t& pos);
-
+void from_proto_position(const Position& from, position_t& to);
 
 /*
  * position_t