More complete rGFA support (experimental)

src/algorithms/coverage_depth.cpp

@@ -295,7 +295,8 @@ pair<double, double> sample_gam_depth(const HandleGraph& graph, const vector<Ali
     return combine_and_average_node_coverages(graph, node_coverages, min_coverage);
 }
 
-void path_depths(const PathHandleGraph& graph, const string& path_name, size_t min_coverage, bool count_cycles, ostream& out_stream) {
+void path_depths(const PathHandleGraph& graph, const string& path_name, size_t min_coverage, bool count_cycles,
+                 const unordered_set<path_handle_t>& ignore_paths, ostream& out_stream) {
     assert(graph.has_path(path_name));
 
     path_handle_t path_handle = graph.get_path_handle(path_name);
@@ -311,10 +312,14 @@ void path_depths(const PathHandleGraph& graph, const string& path_name, size_t m
             size_t step_count = 0;            
             handle_t handle = graph.get_handle_of_step(step_handle);
             graph.for_each_step_on_handle(handle, [&](step_handle_t step_handle_2) {
+                    path_handle_t step_path_handle = graph.get_path_handle_of_step(step_handle_2);
+                    if (ignore_paths.count(step_path_handle)) {
+                        return true;
+                    }
                     if (count_cycles) {
                         ++step_count;
                     } else {
-                        path_handle_t step_path_handle = graph.get_path_handle_of_step(step_handle_2);
+
                         auto it = path_to_name.find(step_path_handle);
                         if (it == path_to_name.end()) {
                             string step_path_name = graph.get_path_name(step_path_handle);
@@ -323,6 +328,7 @@ void path_depths(const PathHandleGraph& graph, const string& path_name, size_t m
                         }
                         path_set.insert(it->second);
                     }
+                    return true;
                 });
             size_t coverage = (count_cycles ? step_count : path_set.size()) - 1;
             size_t node_len = graph.get_length(handle);
@@ -339,7 +345,7 @@ void path_depths(const PathHandleGraph& graph, const string& path_name, size_t m
 
 pair<double, double> path_depth_of_bin(const PathHandleGraph& graph,
                                        step_handle_t start_step, step_handle_t end_plus_one_step,
-                                       size_t min_coverage, bool count_cycles) {
+                                       size_t min_coverage, bool count_cycles, const unordered_set<path_handle_t>& ignore_paths) {
 
     // compute the mean and variance of our base coverage across the bin
     size_t bin_length = 0;
@@ -357,18 +363,22 @@ pair<double, double> path_depth_of_bin(const PathHandleGraph& graph,
         unordered_set<string> path_set;
         size_t step_count = 0;
         graph.for_each_step_on_handle(cur_handle, [&](step_handle_t step_handle) {
+                path_handle_t step_path_handle = graph.get_path_handle_of_step(step_handle);
+                if (ignore_paths.count(step_path_handle)) {
+                    return true;
+                }
                 if (count_cycles) {
                     ++step_count;
                 } else {
-                        path_handle_t step_path_handle = graph.get_path_handle_of_step(step_handle);
                         auto it = path_to_name.find(step_path_handle);
                         if (it == path_to_name.end()) {
                             string step_path_name = graph.get_path_name(step_path_handle);
                             // disregard subpath tags when counting
                             it = path_to_name.insert(make_pair(step_path_handle, Paths::strip_subrange(step_path_name))).first;
                         }
                         path_set.insert(it->second);
-                }                    
+                }
+                return true;
             });
         size_t coverage = (count_cycles ? step_count : path_set.size()) - 1;        
 
@@ -386,7 +396,8 @@ vector<tuple<size_t, size_t, double, double>> binned_path_depth(const PathHandle
                                                                 const string& path_name,
                                                                 size_t bin_size,
                                                                 size_t min_coverage,
-                                                                bool count_cycles) {
+                                                                bool count_cycles,
+                                                                const unordered_set<path_handle_t>& ignore_paths) {
 
     path_handle_t path_handle = graph.get_path_handle(path_name);
 
@@ -414,7 +425,8 @@ vector<tuple<size_t, size_t, double, double>> binned_path_depth(const PathHandle
         step_handle_t bin_end_step = i < bins.size() - 1 ? bins[i+1].second : end_step;
         size_t bin_start = bins[i].first;
         size_t bin_end = i < bins.size() - 1 ? bins[i+1].first : offset;
-        pair<double, double> coverage = path_depth_of_bin(graph, bin_start_step, bin_end_step, min_coverage, count_cycles);
+        pair<double, double> coverage = path_depth_of_bin(graph, bin_start_step, bin_end_step, min_coverage, count_cycles,
+                                                          ignore_paths);
         binned_depths[i] = make_tuple(bin_start, bin_end, coverage.first, coverage.second);
     }
 

src/algorithms/coverage_depth.hpp

@@ -67,14 +67,16 @@ pair<double, double> sample_mapping_depth(const HandleGraph& graph, const vector
 /// print path-name offset base-coverage for every base on a path (just like samtools depth)
 /// ignoring things below min_coverage.  offsets are 1-based in output stream
 /// coverage here is the number of steps from (unique) other paths
-void path_depths(const PathHandleGraph& graph, const string& path_name, size_t min_coverage, bool count_cycles, ostream& out_stream);
+void path_depths(const PathHandleGraph& graph, const string& path_name, size_t min_coverage, bool count_cycles,
+                 const unordered_set<path_handle_t>& ignore_paths, ostream& out_stream);
 
 /// like packed_depth_of_bin (above), but use paths (as in path_depths) for measuring coverage
 pair<double, double> path_depth_of_bin(const PathHandleGraph& graph, step_handle_t start_step, step_handle_t end_plus_one_step,
-                                       size_t min_coverage, bool count_cycles);
+                                       size_t min_coverage, bool count_cycles, const unordered_set<path_handle_t>& ignore_paths);
 
-vector<tuple<size_t, size_t, double, double>> binned_path_depth(const PathHandleGraph& graph, const string& path_name, size_t bin_size,
-                                                                size_t min_coverage, bool count_cycles);
+vector<tuple<size_t, size_t, double, double>> binned_path_depth(const PathHandleGraph& graph, const string& path_name,
+                                                                size_t bin_size, size_t min_coverage, bool count_cycles,
+                                                                const unordered_set<path_handle_t>& ignore_paths);
 
 }
 }

src/algorithms/gfa_to_handle.cpp

@@ -1,6 +1,6 @@
 #include "gfa_to_handle.hpp"
 #include "../path.hpp"
-
+#include "../rgfa.hpp"
 #include <gbwtgraph/utils.h>
 
 namespace vg {
@@ -291,7 +291,6 @@ static void add_path_listeners(GFAParser& parser, MutablePathMutableHandleGraph*
         }
         if (found == rgfa_cache->end()) {
             // Need to make a new path, possibly with subrange start info.
-
             std::pair<int64_t, int64_t> subrange;
             if (offset == 0) {
                 // Don't send a subrange
@@ -300,14 +299,22 @@ static void add_path_listeners(GFAParser& parser, MutablePathMutableHandleGraph*
                 // Start later than 0
                 subrange = std::pair<int64_t, int64_t>(offset, PathMetadata::NO_END_POSITION);
             }
+
+            string rgfa_path_name;
+            if (path_rank > 0) {
+                // Special logic for off-reference paths, which get loaded into special rGFA cover paths
+                rgfa_path_name = RGFACover::make_rgfa_path_name(path_name, offset, length, false);
+            } else {
+                // TODO: See if we can split up the path name into a sample/haplotype/etc. to give it a ref sense.
+                rgfa_path_name = PathMetadata::create_path_name(PathSense::GENERIC,
+                                                                PathMetadata::NO_SAMPLE_NAME,
+                                                                path_name, 
+                                                                PathMetadata::NO_HAPLOTYPE,
+                                                                PathMetadata::NO_PHASE_BLOCK,
+                                                                subrange);
+            }
+            path_handle_t path = graph->create_path_handle(rgfa_path_name);
 
-            // TODO: See if we can split up the path name into a sample/haplotype/etc. to give it a ref sense.
-            path_handle_t path = graph->create_path(PathSense::GENERIC,
-                                                    PathMetadata::NO_SAMPLE_NAME,
-                                                    path_name, 
-                                                    PathMetadata::NO_HAPLOTYPE,
-                                                    PathMetadata::NO_PHASE_BLOCK,
-                                                    subrange);
             // Then cache it
             found = rgfa_cache->emplace_hint(found, path_name, std::make_pair(path, offset));
         }
@@ -667,10 +674,6 @@ tuple<string, size_t, string, pair<int64_t, int64_t>, GFAParser::chars_t, GFAPar
 
     // Grab the haplotype number
     int64_t haplotype_number = take_number(cursor, end, -1, "parsing haplotype number");
-    if (haplotype_number == -1) {
-        // This field is required
-        throw GFAFormatError("Missing haplotype number in W line", cursor);
-    }
     take_tab(cursor, end, "parsing end of haplotype number");
 
     // Grab the sequence/contig/locus name
@@ -1041,7 +1044,13 @@ void GFAParser::parse(istream& in) {
                             } else {
                                 // This path existed already. Make sure we aren't showing a conflicting rank
                                 if (rgfa_path_rank != get<0>(found->second)) {
-                                    throw GFAFormatError("rGFA path " + rgfa_path_name + " has conflicting ranks " + std::to_string(rgfa_path_rank) + " and " + std::to_string(get<0>(found->second)));
+                                  // vg paths can now write different fragments of the same sample/hap with different ranks
+                                  // (unlike minigraph where all fragments of teh same assembly share a rank)
+                                  // by consequence, this check needs to be disabled.
+                                  // TODO: we still want to check that each individual fragment has consisten ranks
+                                  //       so the check should be moved downstream
+                                  //
+                                  // throw GFAFormatError("rGFA path " + rgfa_path_name + " has conflicting ranks " + std::to_string(rgfa_path_rank) + " and " + std::to_string(get<0>(found->second)));
                                 }
                             }
                             auto& visit_queue = get<2>(found->second);

src/clip.cpp

@@ -481,7 +481,8 @@ void clip_low_depth_nodes_and_edges_generic(MutablePathMutableHandleGraph* graph
                                             function<void(function<void(handle_t, const Region*)>)> iterate_handles,
                                             function<void(function<void(edge_t, const Region*)>)> iterate_edges,                                            
                                             int64_t min_depth, const vector<string>& ref_prefixes,
-                                            int64_t min_fragment_len, bool verbose) {
+                                            int64_t min_fragment_len, const unordered_set<path_handle_t>& ignore_paths,
+                                            bool verbose) {
 
     // find all nodes in the snarl that are not on the reference interval (reference path name from containing interval)
     unordered_set<nid_t> to_delete;
@@ -502,14 +503,18 @@ void clip_low_depth_nodes_and_edges_generic(MutablePathMutableHandleGraph* graph
         bool on_ref = false;
         size_t depth = 0;
         graph->for_each_step_on_handle(handle, [&](step_handle_t step_handle) {
+                path_handle_t path_handle = graph->get_path_handle_of_step(step_handle);
+                if (ignore_paths.count(path_handle)) {
+                    return true;
+                }
                 ++depth;
                 if (depth > min_depth || on_ref) {
                     return false;
                 }
                 if (!ref_prefixes.empty() || region) {
                     // if we have a region, do exact comparison to it.
                     // otherwise, do a prefix check against ref_prefix
-                    string path_name = graph->get_path_name(graph->get_path_handle_of_step(step_handle));
+                    string path_name = graph->get_path_name(path_handle);
                     if ((region && region->seq == path_name) || (!region && check_prefixes(path_name))) {
                         on_ref = true;
                         return false;
@@ -540,6 +545,7 @@ void clip_low_depth_nodes_and_edges_generic(MutablePathMutableHandleGraph* graph
     edge_depths.resize(edge_count + 1);
 
     graph->for_each_path_handle([&](path_handle_t path_handle) {
+        if (!ignore_paths.count(path_handle)) {                
             bool is_ref_path = check_prefixes(graph->get_path_name(path_handle));
             handle_t prev_handle;
             bool first = true;
@@ -563,6 +569,7 @@ void clip_low_depth_nodes_and_edges_generic(MutablePathMutableHandleGraph* graph
                     }
                     prev_handle = handle;
                 });
+            }
         });
 
     unordered_set<edge_t> edges_to_delete;
@@ -622,7 +629,7 @@ void clip_low_depth_nodes_and_edges_generic(MutablePathMutableHandleGraph* graph
 }
 
 void clip_low_depth_nodes_and_edges(MutablePathMutableHandleGraph* graph, int64_t min_depth, const vector<string>& ref_prefixes,
-                          int64_t min_fragment_len, bool verbose) {
+                                    int64_t min_fragment_len, const unordered_set<path_handle_t>& ignore_paths, bool verbose) {
 
     function<void(function<void(handle_t, const Region*)>)> iterate_handles = [&] (function<void(handle_t, const Region*)> visit_handle) {
         graph->for_each_handle([&](handle_t handle) {
@@ -636,11 +643,13 @@ void clip_low_depth_nodes_and_edges(MutablePathMutableHandleGraph* graph, int64_
             });        
     };
 
-    clip_low_depth_nodes_and_edges_generic(graph, iterate_handles, iterate_edges, min_depth, ref_prefixes, min_fragment_len, verbose);
+    clip_low_depth_nodes_and_edges_generic(graph, iterate_handles, iterate_edges, min_depth, ref_prefixes, min_fragment_len,
+                                           ignore_paths, verbose);
 }
 
 void clip_contained_low_depth_nodes_and_edges(MutablePathMutableHandleGraph* graph, PathPositionHandleGraph* pp_graph, const vector<Region>& regions,
-                                    SnarlManager& snarl_manager, bool include_endpoints, int64_t min_depth, int64_t min_fragment_len, bool verbose) {
+                                    SnarlManager& snarl_manager, bool include_endpoints, int64_t min_depth, int64_t min_fragment_len,
+                                              const unordered_set<path_handle_t>& ignore_paths, bool verbose) {
 
     function<void(function<void(handle_t, const Region*)>)> iterate_handles = [&] (function<void(handle_t, const Region*)> visit_handle) {
 
@@ -678,7 +687,7 @@ void clip_contained_low_depth_nodes_and_edges(MutablePathMutableHandleGraph* gra
     }
     vector<string> ref_paths_from_regions(ref_path_set.begin(), ref_path_set.end());
 
-    clip_low_depth_nodes_and_edges_generic(graph, iterate_handles, iterate_edges, min_depth, ref_paths_from_regions, min_fragment_len, verbose);
+    clip_low_depth_nodes_and_edges_generic(graph, iterate_handles, iterate_edges, min_depth, ref_paths_from_regions, min_fragment_len, ignore_paths, verbose);
 
 }
 

src/clip.hpp

@@ -64,19 +64,22 @@ void clip_low_depth_nodes_and_edges_generic(MutablePathMutableHandleGraph* graph
                                             function<void(function<void(handle_t, const Region*)>)> iterate_handles,
                                             function<void(function<void(edge_t, const Region*)>)> iterate_edges,
                                             int64_t min_depth, const vector<string>& ref_prefixes,
-                                            int64_t min_fragment_len, bool verbose);
+                                            int64_t min_fragment_len, const unordered_set<path_handle_t>& ignore_paths,
+                                            bool verbose);
 
 /**
  * Run above function on graph
  */
 void clip_low_depth_nodes_and_edges(MutablePathMutableHandleGraph* graph, int64_t min_depth, const vector<string>& ref_prefixes,
-                          int64_t min_fragment_len, bool verbose);
+                                    int64_t min_fragment_len, const unordered_set<path_handle_t>& ignore_paths, bool verbose);
 
 /**
  * Or on contained snarls
  */
-void clip_contained_low_depth_nodes_and_edges(MutablePathMutableHandleGraph* graph, PathPositionHandleGraph* pp_graph, const vector<Region>& regions,
-                                    SnarlManager& snarl_manager, bool include_endpoints, int64_t min_depth, int64_t min_fragment_len, bool verbose);
+void clip_contained_low_depth_nodes_and_edges(MutablePathMutableHandleGraph* graph, PathPositionHandleGraph* pp_graph,
+                                              const vector<Region>& regions,  SnarlManager& snarl_manager, bool include_endpoints,
+                                              int64_t min_depth, int64_t min_fragment_len,
+                                              const unordered_set<path_handle_t>& ignore_paths, bool verbose);
 
 /**
  * clip out deletion edges