simplify / refine path cover table (but not finalized)

Author: glennhickey

src/deconstructor.cpp

@@ -1053,10 +1053,14 @@ bool Deconstructor::deconstruct_site(const handle_t& snarl_start, const handle_t
                 ref_info.parent_allele = 0;
                 ref_info.parent_len = v.alleles[0].length();
                 ref_info.parent_ref_len = v.alleles[0].length();
+                ref_info.parent_path_interval = trav_steps[ref_trav_idx];
+                ref_info.parent_ref_interval = trav_steps[ref_trav_idx];
+                ref_info.lv0_ref_name = v.sequenceName;
                 ref_info.lv0_ref_name = v.sequenceName;
                 ref_info.lv0_ref_start = v.position;
                 ref_info.lv0_ref_len = v.alleles[0].length();
                 ref_info.lv0_alt_len = v.alleles[ref_info.parent_allele].length();
+                ref_info.lv = 0;
             }            
             v.info["PA"].push_back(std::to_string(ref_info.parent_allele));
             v.info["PL"].push_back(std::to_string(ref_info.parent_len));
@@ -1101,14 +1105,20 @@ bool Deconstructor::deconstruct_site(const handle_t& snarl_start, const handle_t
                     }
                 }
             }
+            // update the path cover
+#pragma omp critical (off_ref_info)                
+            {
+                this->update_path_cover(trav_steps, trav_clusters, in_nesting_info->has_ref ? *in_nesting_info : ref_info);
+            }
+            
             // remember the reference path of this variant site
             // for fasta output, along with information about its parent
             int64_t ref_trav = ref_travs[i];
             const PathInterval& ref_path_interval = trav_steps[ref_trav];
             if (!this->ref_paths.count(graph->get_path_name(graph->get_path_handle_of_step(ref_path_interval.first)))) {
 #pragma omp critical (off_ref_info)                
             {
-                this->off_ref_sequences[ref_path_interval] = *in_nesting_info;
+                //this->off_ref_sequences[ref_path_interval] = *in_nesting_info;
             }
             }
         }
@@ -1480,6 +1490,62 @@ void Deconstructor::deconstruct(vector<string> ref_paths, const PathPositionHand
     write_variants(cout, snarl_manager);
 }
 
+void Deconstructor::update_path_cover(const vector<pair<step_handle_t, step_handle_t>>& trav_steps,
+                                      const vector<vector<int>>& traversal_clusters,
+                                      const NestingInfo& nesting_info) const {
+    // for every cluster, add off-reference sequences
+    // todo: are these in the best order? 
+    for (const vector<int>& trav_cluster : traversal_clusters) {
+        if (trav_cluster[0] >= trav_steps.size()) {
+            assert(this->star_allele);
+            continue;
+        }
+        const PathInterval& path_interval = trav_steps.at(trav_cluster[0]);
+        int64_t start = graph->get_position_of_step(path_interval.first);
+        int64_t end = graph->get_position_of_step(path_interval.second);
+        bool reversed = start > end;
+        
+        // scan the interval storing any uncovered sub-intervals
+        vector<PathInterval> sub_intervals;
+        bool open_interval = false;
+        step_handle_t prev_step;
+        string prev_name;
+        for (step_handle_t step = path_interval.first; step != path_interval.second;
+             step = (reversed ? graph->get_previous_step(step) : graph->get_next_step(step))) {
+            if (this->node_cover.count(graph->get_id(graph->get_handle_of_step(step)))) {
+                if (open_interval) {
+                    if (!this->ref_paths.count(prev_name)) {
+                        // expects inclusive interval
+                        step_handle_t end_step = reversed ? graph->get_next_step(step) : graph->get_previous_step(step);
+                        sub_intervals.push_back(make_pair(prev_step, end_step));
+                    }
+                    open_interval = false;
+                }
+            } else {
+                if (open_interval == false) {
+                    prev_step = step;
+                    prev_name = graph->get_path_name(graph->get_path_handle_of_step(prev_step));
+                    open_interval = true;
+                }
+                this->node_cover.insert(graph->get_id(graph->get_handle_of_step(step)));
+            }
+        }
+        
+        if (open_interval && !this->ref_paths.count(prev_name)) {
+            // expects inclusive interval
+            step_handle_t end_step = reversed ? graph->get_next_step(path_interval.second) :
+                graph->get_previous_step(path_interval.second);            
+            sub_intervals.push_back(make_pair(prev_step, end_step));
+        }
+
+        // update the path cover
+        for (const PathInterval& interval : sub_intervals) {
+            // todo: store something
+            this->off_ref_sequences[interval] = nesting_info;
+        }
+    }
+}           
+                                      
 static string resolve_path_name(const PathPositionHandleGraph* graph,
                                 const PathInterval& path_interval,
                                 int64_t& out_start,
@@ -1615,11 +1681,10 @@ void Deconstructor::save_off_ref_sequences(const string& out_fasta_filename) con
                                                          graph->get_handle_of_step(nesting_info.parent_ref_interval.second));
 
 
-        out_metadata_file << path_name << "\t"
+        out_metadata_file << Paths::strip_subrange(path_name) << "\t"
+                          << pos1 << "\t"
+                          << pos2 << "\t"
                           << snarl_name << "\t"
-                          << par_path_name << "\t"
-                          << par_snarl_name << "\t"
-                          << nesting_info.lv << "\t"
                           << nesting_info.lv0_ref_name << "\t"
                           << nesting_info.lv0_ref_start << "\t"
                           << (nesting_info.lv0_ref_start + nesting_info.lv0_ref_len) << "\t"

src/deconstructor.hpp

@@ -142,6 +142,12 @@ class Deconstructor : public VCFOutputCaller {
                                               const vector<string>& trav_to_name,
                                               const vector<int>& gbwt_phases) const;
 
+    // update off-ref path cover
+    void update_path_cover(const vector<pair<step_handle_t, step_handle_t>>& trav_steps,
+                           const vector<vector<int>>& traversal_clusters,
+                           const NestingInfo& nesting_info) const;
+
+
     // the underlying context-getter
     vector<nid_t> get_context(
         step_handle_t start_step,
@@ -217,6 +223,9 @@ class Deconstructor : public VCFOutputCaller {
 
     // keep track of off-reference reference sequences to print to fasta at the end
     mutable unordered_map<PathInterval, NestingInfo> off_ref_sequences;
+
+    // keep track of path cover for computing off_ref_sequences
+    mutable unordered_set<nid_t> node_cover;
 };
 
 

src/traversal_finder.cpp

@@ -40,6 +40,12 @@ string graph_interval_to_string(const HandleGraph* graph, const handle_t& start_
     return handle_to_string(start_handle) + handle_to_string(end_handle);
 }
 
+string path_interval_to_string(const PathHandleGraph* graph, const PathInterval& path_interval) {
+    string s = graph->get_path_name(graph->get_path_handle_of_step(path_interval.first));
+    return s + ":" + graph_interval_to_string(graph, graph->get_handle_of_step(path_interval.first),
+                                              graph->get_handle_of_step(path_interval.second));   
+}
+
 PathBasedTraversalFinder::PathBasedTraversalFinder(const PathHandleGraph& g, SnarlManager& sm) : graph(g), snarlmanager(sm){
 }
 

src/traversal_finder.hpp

@@ -42,6 +42,7 @@ using PathInterval = pair<step_handle_t, step_handle_t>;
 string traversal_to_string(const PathHandleGraph* graph, const Traversal& traversal, int64_t max_steps = 10);
 // replaces pb2json(snarl)
 string graph_interval_to_string(const HandleGraph* graph, const handle_t& start_handle, const handle_t& end_handle);
+string path_interval_to_string(const PathHandleGraph* graph, const PathInterval& path_interval);
 
 /**
  * Represents a strategy for finding traversals of (nested) sites. Polymorphic