Merge pull request #3065 from vgteam/distance-one-node

Distance one node

src/min_distance.cpp

@@ -115,6 +115,38 @@ MinimumDistanceIndex::MinimumDistanceIndex(const HandleGraph* graph,
         }
     }
 
+    auto add_single_nodes = [&](const handle_t& h)-> bool {
+        id_t id = graph->get_id(h); 
+        if (primary_snarl_assignments[id - min_node_id] == 0) {
+            //If this node hasn't already been added to the distance index, make a fake snarl for it
+            handle_t handle = graph->get_handle(id, false);
+            int64_t node_len = graph->get_length(handle);
+
+            //TODO: Do components properly
+
+            //Make a new connected component for this one node
+            size_t component_num = component_to_chain_index.size()+1;
+            component_to_chain_index.resize(curr_component);
+            //Assign it to a chain that doesn't exist?
+            component_to_chain_index[curr_component-1] = chain_indexes.size();
+            component_to_chain_length[curr_component-1] = node_len;
+            //Also assign this node to a connected component
+            node_to_component[id - min_node_id] = component_num;
+
+
+            //Make a snarl index for it
+            size_t snarl_assignment = snarl_indexes.size();
+            primary_snarl_assignments[id-min_node_id] = snarl_assignment+1;
+            primary_snarl_ranks[id - min_node_id] = 1;
+
+            snarl_indexes.emplace_back(0, false, id, id, false, 0, 1, false);
+            snarl_indexes.back().distances[0]  = node_len + 1; 
+        }
+        return true;
+
+    };
+    graph->for_each_handle(add_single_nodes);
+
     #ifdef debugIndex
     //Every node should be assigned to a snarl
     auto check_assignments = [&](const handle_t& h)-> bool {

src/subcommand/index_main.cpp

@@ -86,7 +86,7 @@ void help_index(char** argv) {
          << "    -D, --dump             print the contents of the db to stdout" << endl
          << "    -C, --compact          compact the index into a single level (improves performance)" << endl
          << "snarl distance index options" << endl
-         << "    -s  --snarl-name FILE  load snarls from FILE" << endl
+         << "    -s  --snarl-name FILE  load snarls from FILE (snarls must include trivial snarls)" << endl
          << "    -j  --dist-name FILE   use this file to store a snarl-based distance index" << endl
          << "    -w  --max_dist N       cap beyond which the maximum distance is no longer accurate. If this is not included or is 0, don't build maximum distance index" << endl;
 }

src/unittest/min_distance.cpp

@@ -116,6 +116,30 @@ int64_t min_distance(VG* graph, pos_t pos1, pos_t pos2){
     return shortestDistance == -1 ? -1 : shortestDistance-1;
 };
 
+    TEST_CASE( "One node",
+                   "[min_dist]" ) {
+        VG graph;
+
+        Node* n1 = graph.create_node("GCATGGAGCGTTGAGTGCGGGTG");
+
+        CactusSnarlFinder bubble_finder(graph);
+        SnarlManager snarl_manager = bubble_finder.find_snarls(); 
+
+
+
+        SECTION( "Create distance index" ) {
+
+            MinimumDistanceIndex di (&graph, &snarl_manager);
+#ifdef print
+            di.print_self();
+#endif
+
+            REQUIRE(di.min_distance(make_pos_t(1, false, 0),
+                                   make_pos_t(1, false, 3) ) == 3);
+
+        }
+
+    }
     TEST_CASE( "Create min distance index for simple nested snarl",
                    "[min_dist]" ) {
         VG graph;