Score seeds correctly

src/algorithms/alignment_path_offsets.cpp

@@ -35,7 +35,7 @@ alignment_path_offsets(const PathPositionHandleGraph& graph,
         // Find the position of this end of this mapping
         pos_t mapping_pos = make_pos_t(mapping.position());
         // Find the positions for this end of this Mapping
-        auto pos_offs = algorithms::nearest_offsets_in_paths(&graph, mapping_pos, nearby ? search_limit : -1, path_filter);
+        auto pos_offs = algorithms::nearest_offsets_in_paths(&graph, mapping_pos, nearby ? search_limit : -1, {PathSense::REFERENCE, PathSense::GENERIC}, path_filter);
         for (auto look_at_end : end) {
             // For the start and the end of the Mapping, as needed
             for (auto& p : pos_offs) {
@@ -84,7 +84,7 @@ multipath_alignment_path_offsets(const PathPositionHandleGraph& graph,
         for (size_t j = 0; j < subpath.path().mapping_size(); ++j) {
             // get the positions on paths that this mapping touches
             pos_t mapping_pos = make_pos_t(subpath.path().mapping(j).position());
-            subpath_search_results[j] = nearest_offsets_in_paths(&graph, mapping_pos, 0, path_filter);
+            subpath_search_results[j] = nearest_offsets_in_paths(&graph, mapping_pos, 0, {PathSense::REFERENCE, PathSense::GENERIC}, path_filter);
             // make sure that offsets are stored in increasing order
             for (pair<const path_handle_t, vector<pair<size_t, bool>>>& search_record : subpath_search_results[j]) {
                 sort(search_record.second.begin(), search_record.second.end());

src/algorithms/chain_items.cpp

@@ -13,6 +13,7 @@
 
 //#define debug_chaining
 //#define debug_transition
+//#define debug_dp
 
 namespace vg {
 namespace algorithms {
@@ -564,8 +565,10 @@ TracedScore chain_items_dp(vector<TracedScore>& chain_scores,
     // We will run this over every transition in a good DP order.
     auto iteratee = [&](const transition_info& transition) {
         if (show_work) {
+#ifdef debug_dp
             cerr << "DP: " << transition.from_anchor << "->" << transition.to_anchor 
                  << " rd " << transition.read_distance << " gd " << transition.graph_distance << endl;
+#endif
         }
 
         crash_unless(chain_scores.size() > transition.to_anchor);
@@ -590,8 +593,10 @@ TracedScore chain_items_dp(vector<TracedScore>& chain_scores,
         auto& source = to_chain[transition.from_anchor];
 
         if (show_work) {
+#ifdef debug_dp
             cerr << "\t\tCome from score " << chain_scores[transition.from_anchor]
                 << " across " << source << " to " << here << endl;
+#endif
         }
 
         // How much does it pay (+) or cost (-) to make the jump from there
@@ -608,9 +613,11 @@ TracedScore chain_items_dp(vector<TracedScore>& chain_scores,
         size_t possible_match_length = std::min(transition.read_distance, transition.graph_distance);
 
         if (show_work) {
+#ifdef debug_dp
             cerr << "\t\t\tFor read distance " << transition.read_distance << " and graph distance " << transition.graph_distance
                  << " an indel of length " << indel_length
                  << ((transition.read_distance > transition.graph_distance) ? " seems plausible" : " would be required") << endl;
+#endif
         }
 
         if (indel_length > max_indel_bases) {
@@ -661,8 +668,10 @@ TracedScore chain_items_dp(vector<TracedScore>& chain_scores,
             chain_scores[transition.to_anchor] = std::max(chain_scores[transition.to_anchor], from_source_score);
 
             if (show_work) {
+#ifdef debug_dp
                 cerr << "\t\tWe can reach #" << transition.to_anchor << " with " << source_score << " + " << jump_points
                      << " from transition + " << item_points << " from item = " << from_source_score << endl;
+#endif
             }
 
             if (diagram) {
@@ -692,7 +701,9 @@ TracedScore chain_items_dp(vector<TracedScore>& chain_scores,
             }
         } else {
             if (show_work) {
+#ifdef debug_dp
                 cerr << "\t\tTransition is impossible." << endl;
+#endif
             }
         }
     };
@@ -751,7 +762,9 @@ TracedScore chain_items_dp(vector<TracedScore>& chain_scores,
         best_score.max_in(chain_scores, to_anchor);
 
         if (show_work) {
+#ifdef debug_dp
             cerr << "\tBest chain end so far: " << best_score << endl;
+#endif
         }
 
     }

src/algorithms/nearest_offsets_in_paths.cpp

@@ -18,6 +18,7 @@ using namespace std;
 path_offset_collection_t nearest_offsets_in_paths(const PathPositionHandleGraph* graph,
                                                   const pos_t& pos,
                                                   int64_t max_search,
+                                                  const std::unordered_set<PathSense>& desired_senses,
                                                   const std::function<bool(const path_handle_t&)>* path_filter) {
 
     // init the return value
@@ -48,7 +49,7 @@ path_offset_collection_t nearest_offsets_in_paths(const PathPositionHandleGraph*
         << " in " << (search_left ? "leftward" : "rightward") << " direction at distance " << dist << endl;
 #endif
 
-        graph->for_each_step_of_sense(here, {PathSense::REFERENCE, PathSense::GENERIC}, [&](const step_handle_t& step) {
+        graph->for_each_step_of_sense(here, desired_senses, [&](const step_handle_t& step) {
             // For each path visit that occurs on this node
 #ifdef debug
             cerr << "handle is on step at path offset " << graph->get_position_of_step(step) << endl;

src/algorithms/nearest_offsets_in_paths.hpp

@@ -33,9 +33,10 @@ using path_offset_collection_t = unordered_map<path_handle_t, vector<pair<size_t
 ///
 /// If path_filter is set, ignores paths for which it returns false.
 ///
-/// Doesn't consider haplotype paths.
+/// Doesn't consider haplotype paths. by default.
 path_offset_collection_t nearest_offsets_in_paths(const PathPositionHandleGraph* graph,
                                                   const pos_t& pos, int64_t max_search,
+                                                  const std::unordered_set<PathSense>& desired_senses = {PathSense::REFERENCE, PathSense::GENERIC},
                                                   const std::function<bool(const path_handle_t&)>* path_filter = nullptr);
 
 /// Wrapper for the above to support some earlier code. Only looks for paths

src/minimizer_mapper.hpp

@@ -427,6 +427,10 @@ class MinimizerMapper : public AlignerClient {
     /// Track linear reference position for placements in log output.
     static constexpr bool default_track_position = false;
     bool track_position = default_track_position;
+
+    /// Track positions along haplotypes
+    static constexpr bool default_haplotype_positions = false;
+    bool haplotype_positions = default_haplotype_positions;
 
     /// If set, log what the mapper is thinking in its mapping of each read.
     static constexpr bool default_show_work = false;
@@ -1467,7 +1471,8 @@ class MinimizerMapper : public AlignerClient {
                                    const std::vector<std::vector<size_t>>& chains,
                                    const std::vector<std::vector<bool>>& chain_rec_flags,
                                    const std::vector<size_t>& chain_source_tree,
-                                   const PathPositionHandleGraph* path_graph);
+                                   const PathPositionHandleGraph* path_graph,
+                                   bool haplotype_positions);
 
     /// Dump a graph
     static void dump_debug_graph(const HandleGraph& graph);

src/minimizer_mapper_from_chains.cpp

@@ -194,7 +194,8 @@ void MinimizerMapper::dump_debug_chains(const ZipCodeForest& zip_code_forest,
                                          const std::vector<std::vector<size_t>>& chains,
                                          const std::vector<std::vector<bool>>& chain_rec_flags,
                                          const std::vector<size_t>& chain_source_tree,
-                                         const PathPositionHandleGraph* path_graph) {
+                                         const PathPositionHandleGraph* path_graph,
+                                         bool haplotype_positions) {
     if (!path_graph) {
         // We don't have a path positional graph for this
         return;
@@ -242,6 +243,10 @@ void MinimizerMapper::dump_debug_chains(const ZipCodeForest& zip_code_forest,
 
         // Determine the positions of all the involved seeds.
         std::unordered_map<size_t, algorithms::path_offset_collection_t> seed_positions;
+        std::unordered_set<PathSense> wanted_senses {PathSense::REFERENCE, PathSense::GENERIC};
+        if (haplotype_positions) {
+            wanted_senses.insert(PathSense::HAPLOTYPE);
+        }
         for (auto& kv : seed_sets) {
             for (const std::vector<size_t> included_seeds : kv.second) {
                 for (auto& seed_num : included_seeds) {
@@ -251,7 +256,7 @@ void MinimizerMapper::dump_debug_chains(const ZipCodeForest& zip_code_forest,
                     auto found = seed_positions.find(seed_num);
                     if (found == seed_positions.end()) {
                         // If we don't know the seed's positions yet, get them
-                        found = seed_positions.emplace_hint(found, seed_num, algorithms::nearest_offsets_in_paths(path_graph, seed.pos, 100));
+                        found = seed_positions.emplace_hint(found, seed_num, algorithms::nearest_offsets_in_paths(path_graph, seed.pos, 100, wanted_senses));
                         for (auto& handle_and_positions : found->second) {
                             std::string path_name = path_graph->get_path_name(handle_and_positions.first);
                             for (auto& position : handle_and_positions.second) {
@@ -268,6 +273,19 @@ void MinimizerMapper::dump_debug_chains(const ZipCodeForest& zip_code_forest,
                                 seedpos.field(ss.str());
                             }
                         }
+                        if (found->second.empty()) {
+                            // The seed doesn't have any linear positions, but might still participate in the winning chain traceback.
+                            // Report it.
+                            seedpos.line();
+                            seedpos.field(seed_anchors.at(seed_num).read_start());
+                            seedpos.field("");
+                            seedpos.field("");
+                            seedpos.field("");
+                            seedpos.field(seed_num);
+                            std::stringstream ss;
+                            ss << seed_anchors.at(seed_num);
+                            seedpos.field(ss.str());
+                        }
                     }
                 }
             }
@@ -735,7 +753,7 @@ vector<Alignment> MinimizerMapper::map_from_chains(Alignment& aln) {
     }
 #endif
 
-    // Turn all the seeds into anchors. Either we'll fragment them directly or
+    // Turn all the seeds into anchors. Either we'll chain them directly or
     // use them to make gapless extension anchors over them.
     // TODO: Can we only use the seeds that are in trees we keep?
     vector<algorithms::Anchor> seed_anchors = this->to_anchors(aln, minimizers, seeds);
@@ -781,7 +799,7 @@ vector<Alignment> MinimizerMapper::map_from_chains(Alignment& aln) {
 
     // Dump all chains if requested (do this before alignments, while chains still exist)
     if (show_work && !chains.empty() && this->path_graph != nullptr) {
-        dump_debug_chains(zip_code_forest, seeds, minimizers, seed_anchors, chains, chain_rec_flags, chain_source_tree, this->path_graph);
+        dump_debug_chains(zip_code_forest, seeds, minimizers, seed_anchors, chains, chain_rec_flags, chain_source_tree, this->path_graph, this->haplotype_positions);
     }
 
     if (alignments.size() == 0) {
@@ -1630,9 +1648,13 @@ void MinimizerMapper::do_chaining_on_trees(Alignment& aln, const ZipCodeForest&
                     // Add position annotations for the good-looking chains.
                     // Should be much faster than full correctness tracking from every seed.
                     crash_unless(this->path_graph);
+                    std::unordered_set<PathSense> wanted_senses {PathSense::REFERENCE, PathSense::GENERIC};
+                    if (haplotype_positions) {
+                        wanted_senses.insert(PathSense::HAPLOTYPE);
+                    }
                     for (auto& boundary : {anchor_view[scored_chain.second.front()].graph_start(), anchor_view[scored_chain.second.back()].graph_end()}) {
                         // For each end of the chain
-                        auto offsets = algorithms::nearest_offsets_in_paths(this->path_graph, boundary, 100);
+                        auto offsets = algorithms::nearest_offsets_in_paths(this->path_graph, boundary, 100, wanted_senses);
                         for (auto& handle_and_positions : offsets) {
                             for (auto& position : handle_and_positions.second) {
                                 // Tell the funnel all the effective positions, ignoring orientation
@@ -2706,7 +2728,11 @@ Alignment MinimizerMapper::find_chain_alignment(
             #pragma omp critical (cerr)
             {
                 cerr << log_name() << "Need to align graph from " << (*here).graph_end() << " to " << (*next).graph_start()
-                    << " separated by ~" << graph_dist << " bp and sequence \"" << linking_bases << "\"" << endl;
+                    << " separated by ~" << graph_dist << " bp";
+                if (linking_bases.size() < 200) {
+                    cerr << " and sequence \"" << linking_bases << "\"";
+                }
+                cerr << endl;
             }
         }
 #endif
@@ -2894,7 +2920,7 @@ Alignment MinimizerMapper::find_chain_alignment(
         if (show_work) {
             #pragma omp critical (cerr)
             {
-                cerr << log_name() << "Aligned and added link of " << link_length << " bp read and " << link_graph_path_length << " bp graph via " << link_alignment_source << " with score of " << link_score << std::endl;
+                cerr << log_name() << "Aligned and added link to " << *next << " of " << link_length << " bp read and " << link_graph_path_length << " bp graph via " << link_alignment_source << " with score of " << link_score << std::endl;
             }
         }
 
@@ -3828,7 +3854,7 @@ algorithms::Anchor MinimizerMapper::to_anchor(const Alignment& aln, const Vector
     // Work out how many points the anchor is.
     // TODO: Always make sequence and quality available for scoring!
     // We're going to score the anchor as the full minimizer, and rely on the margins to stop us from taking overlapping anchors.
-    int score = aligner->score_exact_match(aln, read_start - margin_left, length + margin_right);
+    int score = aligner->score_exact_match(aln, read_start - margin_left, margin_left + length + margin_right);
     return algorithms::Anchor(read_start, graph_start, length, margin_left, margin_right, score, seed_number, &(seed.zipcode), hint_start, source.is_repetitive, paths); 
 }