fix assertion fails that happen in full hprc graph

vgteam · Mar 24, 2023 · c8b9d55 · c8b9d55 · adamnovak · Mar 24, 2023
1 parent da89deb
commit c8b9d55
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 7 deletions.
diff --git a/src/gfa.cpp b/src/gfa.cpp
@@ -539,7 +539,9 @@ void rgfa_snarl_cover(const PathHandleGraph* graph,
         }
     }
 
-    if (ref_paths.empty() && !cover_node_to_fragment.count(snarl.start().node_id())) {
+    // note: checking both snarl endpoint here is actually necessary: if the reference path doesn't end in a tip,
+    // you can end up with a trivial snarl at its end which will crash on this test. 
+    if (ref_paths.empty() && (!cover_node_to_fragment.count(snarl.start().node_id()) || !cover_node_to_fragment.count(snarl.end().node_id()))) {      
         // we're not nested in a reference snarl, and we have no reference path
         // by the current logic, there's nothing to be done.
         cerr << "[rgfa] warning: No referene path through snarl " 
@@ -613,11 +615,20 @@ void rgfa_snarl_cover(const PathHandleGraph* graph,
         vector<pair<int64_t, int64_t>> uncovered_intervals = get_uncovered_intervals(trav);
 
         for (const auto& uncovered_interval : uncovered_intervals) {
+            unordered_set<nid_t> cycle_check;
+            bool cyclic = false;
             int64_t interval_length = 0;            
-            for (int64_t i = uncovered_interval.first; i < uncovered_interval.second; ++i) {
-                interval_length += graph->get_length(graph->get_handle_of_step(trav[i]));
+            for (int64_t i = uncovered_interval.first; i < uncovered_interval.second && !cyclic; ++i) {
+                handle_t handle = graph->get_handle_of_step(trav[i]);
+                interval_length += graph->get_length(handle);
+                nid_t node_id = graph->get_id(handle);
+                if (cycle_check.count(node_id)) {
+                    cyclic = true;
+                } else {
+                    cycle_check.insert(node_id);
+                }                
             }
-            if (interval_length >= minimum_length) {
+            if (!cyclic && interval_length >= minimum_length) {
                 auto trav_stats = rgfa_traversal_stats(graph, trav, uncovered_interval);
                 ranked_trav_fragments.push_back(make_pair(trav_stats, make_pair(trav_idx, uncovered_interval)));
             }

diff --git a/src/subcommand/paths_main.cpp b/src/subcommand/paths_main.cpp
@@ -140,7 +140,8 @@ int main_paths(int argc, char** argv) {
             {"extract-vg", no_argument, 0, 'V'},
             {"drop-paths", no_argument, 0, 'd'},
             {"retain-paths", no_argument, 0, 'r'},
-            {"rgfa-cover", required_argument, 0, 'R'},            
+            {"rgfa-cover", required_argument, 0, 'R'},
+            {"snarls", required_argument, 0, 's'},
             {"extract-gam", no_argument, 0, 'X'},
             {"extract-gaf", no_argument, 0, 'A'},            
             {"list", no_argument, 0, 'L'},
@@ -162,7 +163,7 @@ int main_paths(int argc, char** argv) {
         };
 
         int option_index = 0;
-        c = getopt_long (argc, argv, "hLXv:x:g:Q:VEMCFAS:Tq:drR:aGp:ct:",
+        c = getopt_long (argc, argv, "hLXv:x:g:Q:VEMCFAS:Tq:drR:s:aGp:ct:",
                 long_options, &option_index);
 
         // Detect the end of the options.
@@ -202,7 +203,11 @@ int main_paths(int argc, char** argv) {
             rgfa_min_len = parse<int64_t>(optarg);            
             output_formats++;
             break;
-
+
+        case 's':
+            snarl_filename = optarg;
+            break;
+
         case 'X':
             extract_as_gam = true;
             output_formats++;