No secret paths

Makefile

@@ -43,6 +43,7 @@ include $(wildcard $(IO_SHARED_OBJ_DIR)/*.d)
 include $(wildcard $(SUBCOMMAND_OBJ_DIR)/*.d)
 include $(wildcard $(UNITTEST_OBJ_DIR)/*.d)
 include $(wildcard $(UNITTEST_BIN_DIR)/*.d)
+include $(wildcard test/*.d)
 
 # What pkg-config-controlled system dependencies should we use compile and link flags from?
 # Use PKG_CONFIG_PATH to point the build system at the right versions of these, if they aren't picked up automatically.
@@ -606,8 +607,8 @@ else
 	ln -s `which shuf` $(BIN_DIR)/shuf
 endif
 
-test/build_graph: test/build_graph.cpp $(LIB_DIR)/libvg.a $(SRC_DIR)/vg.hpp
-	$(CXX) $(INCLUDE_FLAGS) $(CPPFLAGS) $(CXXFLAGS) -o test/build_graph test/build_graph.cpp $(LD_LIB_DIR_FLAGS) $(LDFLAGS) $(LIB_DIR)/libvg.a $(LD_LIB_FLAGS) $(START_STATIC) $(LD_STATIC_LIB_FLAGS) $(END_STATIC) $(FILTER)
+test/build_graph: test/build_graph.o $(LIB_DIR)/libvg.a
+	$(CXX) $(INCLUDE_FLAGS) $(CPPFLAGS) $(CXXFLAGS) -o test/build_graph $(PRE_LINK_DEPS) test/build_graph.o $(LD_LIB_DIR_FLAGS) $(LDFLAGS) $(LIB_DIR)/libvg.a $(LD_LIB_FLAGS) $(START_STATIC) $(LD_STATIC_LIB_FLAGS) $(END_STATIC) $(LD_STATIC_LIB_DEPS) $(LD_EXE_LIB_FLAGS)
 
 $(LIB_DIR)/mimalloc.o: $(MIMALLOC_DIR)/src/*.c $(MIMALLOC_DIR)/src/*/*.c $(MIMALLOC_DIR)/src/*/*/*.c $(MIMALLOC_DIR)/include/*.h $(MIMALLOC_DIR)/include/*/*.h $(MIMALLOC_DIR)/CMakeLists.txt
 	+rm -f $(LIB_DIR)/mimalloc.o && rm -Rf $(INC_DIR)/mimalloc $(INC_DIR)/mimalloc*.h && cd $(MIMALLOC_DIR) && rm -Rf build && mkdir build && cd build && cmake -DCMAKE_C_COMPILER="$(CC)" -DCMAKE_CXX_COMPILER="$(CXX)" -DCMAKE_C_FLAGS="$(CFLAGS)" -DCMAKE_CXX_FLAGS="$(CXXFLAGS)" ..  && $(MAKE) $(FILTER) &&  cp -r ../include/* $(CWD)/$(INC_DIR)/ && cp mimalloc.o $(CWD)/$(LIB_DIR)/
@@ -708,9 +709,11 @@ $(LIB_DIR)/cleaned_old_elfutils:
 $(LIB_DIR)/libvgio.a: $(LIB_DIR)/libhts.a $(LIB_DIR)/libhandlegraph.a $(LIB_DIR)/pkgconfig/htslib.pc $(LIB_DIR)/cleaned_old_protobuf_v003 $(LIBVGIO_DIR)/CMakeLists.txt $(LIBVGIO_DIR)/src/*.cpp $(LIBVGIO_DIR)/include/vg/io/*.hpp $(LIBVGIO_DIR)/deps/vg.proto
 	+rm -f $(CWD)/$(INC_DIR)/vg.pb.h $(CWD)/$(INC_DIR)/vg/vg.pb.h
 	+rm -Rf $(CWD)/$(INC_DIR)/vg/io/
-	+export CXXFLAGS="$(CPPFLAGS) $(CXXFLAGS)" && export LDFLAGS="$(LD_LIB_DIR_FLAGS) $(LDFLAGS)" && cd $(LIBVGIO_DIR) && rm -Rf CMakeCache.txt CMakeFiles *.cmake install_manifest.txt *.pb.cc *.pb.h *.a && rm -rf build-vg && mkdir build-vg && cd build-vg && PKG_CONFIG_PATH=$(CWD)/$(LIB_DIR)/pkgconfig:$(PKG_CONFIG_PATH) cmake -DCMAKE_C_COMPILER="$(CC)" -DCMAKE_CXX_COMPILER="$(CXX)" -DCMAKE_CXX_STANDARD=$(CXX_STANDARD) -DCMAKE_VERBOSE_MAKEFILE=ON -DCMAKE_PREFIX_PATH="/usr;$(OMP_PREFIXES)" -DCMAKE_INSTALL_PREFIX=$(CWD) -DCMAKE_INSTALL_LIBDIR=lib .. $(FILTER) && $(MAKE) clean && VERBOSE=1 $(MAKE) $(FILTER) && $(MAKE) install
+	+export CXXFLAGS="$(CPPFLAGS) $(CXXFLAGS)" && export LDFLAGS="$(LD_LIB_DIR_FLAGS) $(LDFLAGS)" && cd $(LIBVGIO_DIR) && rm -Rf CMakeCache.txt CMakeFiles *.cmake install_manifest.txt *.pb.cc *.pb.h *.a && rm -rf build-vg && mkdir build-vg && cd build-vg && PKG_CONFIG_PATH=$(CWD)/$(LIB_DIR)/pkgconfig:$(PKG_CONFIG_PATH) cmake -DCMAKE_C_COMPILER="$(CC)" -DCMAKE_CXX_COMPILER="$(CXX)" -DCMAKE_CXX_STANDARD=$(CXX_STANDARD) -DCMAKE_VERBOSE_MAKEFILE=ON -DCMAKE_PREFIX_PATH="/usr;$(OMP_PREFIXES)" -DCMAKE_INSTALL_PREFIX=$(CWD) -DCMAKE_INSTALL_LIBDIR=lib -DUSE_INSTALLED_LIBHANDLEGRAPH_ONLY=ON .. $(FILTER) && $(MAKE) clean && VERBOSE=1 $(MAKE) $(FILTER) && $(MAKE) install
 
 $(LIB_DIR)/libhandlegraph.a: $(LIBHANDLEGRAPH_DIR)/src/include/handlegraph/*.hpp $(LIBHANDLEGRAPH_DIR)/src/*.cpp
+	+rm -f $(LIB_DIR)/libhandlegraph.a $(LIB_DIR)/libhandlegraph.$(SHARED_SUFFIX)
+	+rm -Rf $(INC_DIR)/handlegraph $(LIB_DIR)/cmake/libhandlegraph
 	+cd $(LIBHANDLEGRAPH_DIR) && rm -Rf build CMakeCache.txt CMakeFiles && mkdir build && cd build && CXXFLAGS="$(CXXFLAGS) $(CPPFLAGS)" cmake -DCMAKE_C_COMPILER="$(CC)" -DCMAKE_CXX_COMPILER="$(CXX)" -DCMAKE_VERBOSE_MAKEFILE=ON -DCMAKE_INSTALL_PREFIX=$(CWD) -DCMAKE_INSTALL_LIBDIR=lib .. && $(MAKE) $(FILTER) && $(MAKE) install
 
 
@@ -1001,6 +1004,9 @@ $(UNITTEST_OBJ): $(UNITTEST_OBJ_DIR)/%.o : $(UNITTEST_SRC_DIR)/%.cpp $(UNITTEST_
 $(UNITTEST_SUPPORT_OBJ): $(UNITTEST_SUPPORT_OBJ_DIR)/%.o : $(UNITTEST_SUPPORT_SRC_DIR)/%.cpp $(UNITTEST_SUPPORT_OBJ_DIR)/%.d $(DEPS)
 	$(CXX) $(INCLUDE_FLAGS) $(CPPFLAGS) $(CXXFLAGS) $(DEPGEN_FLAGS) -c -o $@ $< $(FILTER)
 	@touch $@
+test/build_graph.o: test/%.o : test/%.cpp test/%.d $(DEPS)
+	$(CXX) $(INCLUDE_FLAGS) $(CPPFLAGS) $(CXXFLAGS) $(DEPGEN_FLAGS) -c -o $@ $< $(FILTER)
+	@touch $@
 
 # Config objects get individual rules
 $(CONFIG_OBJ_DIR)/allocator_config_mimalloc.o: $(CONFIG_SRC_DIR)/allocator_config_mimalloc.cpp $(CONFIG_OBJ_DIR)/allocator_config_mimalloc.d $(DEPS) $(LIB_DIR)/mimalloc.o
@@ -1023,10 +1029,10 @@ $(CONFIG_OBJ_DIR)/%.d: ;
 $(IO_OBJ_DIR)/%.d: ;
 $(SUBCOMMAND_OBJ_DIR)/%.d: ;
 $(UNITTEST_OBJ_DIR)/%.d: ;
+test/%.d: ;
 
 # Don't delete them.
-.PRECIOUS: $(OBJ_DIR)/%.d $(ALGORITHMS_OBJ_DIR)/%.d $(CONFIG_OBJ_DIR)/%.d $(IO_OBJ_DIR)/%.d $(SUBCOMMAND_OBJ_DIR)/%.d $(UNITTEST_OBJ_DIR)/%.d
-
+.PRECIOUS: $(OBJ_DIR)/%.d $(ALGORITHMS_OBJ_DIR)/%.d $(CONFIG_OBJ_DIR)/%.d $(IO_OBJ_DIR)/%.d $(SUBCOMMAND_OBJ_DIR)/%.d $(UNITTEST_OBJ_DIR)/%.d test/%.d
 # Use no implicit rules
 .SUFFIXES:
 

src/algorithms/alignment_path_offsets.hpp

@@ -23,6 +23,8 @@ using namespace std;
 ///
 /// If path_filter is set, and it returns false for a path, that path is not
 /// used to annotate the read.
+///
+/// Doesn't consider haplotype paths.
 unordered_map<path_handle_t, vector<pair<size_t, bool> > >
 alignment_path_offsets(const PathPositionHandleGraph& graph,
                        const Alignment& aln,
@@ -37,7 +39,8 @@ alignment_path_offsets(const PathPositionHandleGraph& graph,
 /// is disconnected or fans out toward the sources or sinks.
 ///
 /// If path_filter is set, and it returns false for a path, that path is not
-/// used to annotate the read.
+/// used to annotate the read. If it is not set, all reference and
+/// generic-sense paths are used.
 unordered_map<path_handle_t, vector<pair<size_t, bool> > >
 multipath_alignment_path_offsets(const PathPositionHandleGraph& graph,
                                  const multipath_alignment_t& mp_aln,
@@ -53,6 +56,8 @@ multipath_alignment_path_offsets(const PathPositionHandleGraph& graph,
 ///
 /// If path_filter is set, and it returns false for a path, that path is not
 /// used to annotate the read.
+///
+/// Doesn't consider haplotype paths.
 void annotate_with_initial_path_positions(const PathPositionHandleGraph& graph, Alignment& aln, int64_t search_limit = 0, const std::function<bool(const path_handle_t&)>* path_filter = nullptr);
 
 /// Use the graph to annotate an Alignment with the first
@@ -67,6 +72,8 @@ void annotate_with_initial_path_positions(const PathPositionHandleGraph& graph,
 ///
 /// If path_filter is set, and it returns false for a path, that path is not
 /// used to annotate the read.
+///
+/// Doesn't consider haplotype paths.
 void annotate_with_node_path_positions(const PathPositionHandleGraph& graph, Alignment& aln, int64_t search_limit = 0, const std::function<bool(const path_handle_t&)>* path_filter = nullptr);
 
 /// Use the graph to annotate an Alignment with positions on each reference
@@ -81,6 +88,8 @@ void annotate_with_node_path_positions(const PathPositionHandleGraph& graph, Ali
 ///
 /// If path_filter is set, and it returns false for a path, that path is not
 /// used to annotate the read.
+///
+/// Doesn't consider haplotype paths.
 void annotate_with_path_positions(const PathPositionHandleGraph& graph, Alignment& aln, bool just_min, int64_t search_limit = 0, const std::function<bool(const path_handle_t&)>* path_filter = nullptr);
 
 /// Use the graph annotate Alignments with the first position
@@ -93,6 +102,8 @@ void annotate_with_path_positions(const PathPositionHandleGraph& graph, Alignmen
 ///
 /// If path_filter is set, and it returns false for a path, that path is not
 /// used to annotate the read.
+///
+/// Doesn't consider haplotype paths.
 void annotate_with_initial_path_positions(const PathPositionHandleGraph& graph, vector<Alignment>& aln, int64_t search_limit = 0, const std::function<bool(const path_handle_t&)>* path_filter = nullptr);
 
 

src/algorithms/nearest_offsets_in_paths.cpp

@@ -47,8 +47,8 @@ path_offset_collection_t nearest_offsets_in_paths(const PathPositionHandleGraph*
         cerr << "traversing " << graph->get_id(here) << (graph->get_is_reverse(here) ? "-" : "+")
         << " in " << (search_left ? "leftward" : "rightward") << " direction at distance " << dist << endl;
 #endif
-        
-        for (const step_handle_t& step : graph->steps_of_handle(here)) {
+
+        graph->for_each_step_of_sense(here, {PathSense::REFERENCE, PathSense::GENERIC}, [&](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;
@@ -61,9 +61,9 @@ path_offset_collection_t nearest_offsets_in_paths(const PathPositionHandleGraph*
 #ifdef debug
                 cerr << "handle is on ignored path " << graph->get_path_name(path_handle) << endl;
 #endif
-                continue;
+                return;
             }
-            
+
             // flip the handle back to the orientation it started in
             handle_t oriented = search_left ? graph->flip(here) : here;
 
@@ -89,7 +89,7 @@ path_offset_collection_t nearest_offsets_in_paths(const PathPositionHandleGraph*
 
             // add in the search distance and add the result to the output
             return_val[path_handle].emplace_back(path_offset, rev_on_path);
-        }
+        });
 
         if (!return_val.empty()) {
             // we found the closest, we're done

src/algorithms/nearest_offsets_in_paths.hpp

@@ -32,6 +32,8 @@ using path_offset_collection_t = unordered_map<path_handle_t, vector<pair<size_t
 /// Stops search when path(s) are ancountered.
 ///
 /// If path_filter is set, ignores paths for which it returns false.
+///
+/// Doesn't consider haplotype paths.
 path_offset_collection_t nearest_offsets_in_paths(const PathPositionHandleGraph* graph,
                                                   const pos_t& pos, int64_t max_search,
                                                   const std::function<bool(const path_handle_t&)>* path_filter = nullptr);

src/clip.cpp

@@ -123,7 +123,7 @@ void visit_contained_snarls(const PathPositionHandleGraph* graph, const vector<R
         path_name_set.insert(region.seq);
     }
     unordered_set<string> graph_path_name_set;
-    graph->for_each_path_handle([&](path_handle_t path_handle) {
+    graph->for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](path_handle_t path_handle) {
             string graph_path_name = graph->get_path_name(path_handle);
             if (path_name_set.count(Paths::strip_subrange(graph_path_name))) {
                 graph_path_name_set.insert(graph_path_name);
@@ -614,7 +614,7 @@ void clip_low_depth_nodes_and_edges_generic(MutablePathMutableHandleGraph* graph
     bdsg::PackedVector<> edge_depths;
     edge_depths.resize(edge_count + 1);
 
-    graph->for_each_path_handle([&](path_handle_t path_handle) {
+    graph->for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](path_handle_t path_handle) {
             bool is_ref_path = check_prefixes(graph->get_path_name(path_handle));
             handle_t prev_handle;
             bool first = true;
@@ -1053,7 +1053,7 @@ void clip_deletion_edges(MutablePathMutableHandleGraph* graph, int64_t max_delet
 
     // load up the reference paths and their ids
     unordered_set<path_handle_t> ref_paths;
-    graph->for_each_path_handle([&](path_handle_t path_handle) {
+    graph->for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](path_handle_t path_handle) {
             string path_name = graph->get_path_name(path_handle);
             for (const string& ref_prefix : ref_prefixes) {
                 if (path_name.compare(0, ref_prefix.length(), ref_prefix) == 0) {
@@ -1312,7 +1312,7 @@ void clip_contained_stubs(MutablePathMutableHandleGraph* graph, PathPositionHand
 void stubbify_ref_paths(MutablePathMutableHandleGraph* graph, const vector<string>& ref_prefixes, int64_t min_fragment_len, bool verbose) {
     unordered_set<edge_t> edges_to_delete;
     int64_t stubbified_path_count = 0; // just for logging
-    graph->for_each_path_handle([&](path_handle_t path_handle) {
+    graph->for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](path_handle_t path_handle) {
         string path_name = graph->get_path_name(path_handle);
         for (const string& ref_prefix : ref_prefixes) {
             bool was_stubbified = false;

src/deconstructor.cpp

@@ -1158,14 +1158,31 @@ string Deconstructor::get_vcf_header() {
     if (!long_ref_contig) {
         long_ref_contig = ref_samples.size() > 1 || ref_haplotypes.size() > 1 || nested_decomposition;
     }
-    this->long_ref_contig = long_ref_contig;
     sample_names.clear();
     unordered_map<string, set<int>> sample_to_haps;
 
-    // find sample names from non-reference paths
-    graph->for_each_path_handle([&](const path_handle_t& path_handle) {
+    // We always want generic and reference paths, since we do our own designated deconstruction reference dropping.
+    std::unordered_set<PathSense> wanted_senses{PathSense::GENERIC, PathSense::REFERENCE};
+
+    // If we're using a GBWT, we don't want to include any haplotypes in the
+    // input graph, because we want the haplotypes from the GBWT instead (even
+    // if the haplotypes in ther input graph are actually backed by the
+    // provided GBWT). If we aren't, we want the haplotypes from the input
+    // graph included.
+    if (!gbwt) {
+        wanted_senses.insert(PathSense::HAPLOTYPE);
+    }
+
+    // TODO: Should we unify around using haplotype-sense paths and not the
+    // GBWT API?
+
+    // find sample names from paths
+    graph->for_each_path_of_sense(wanted_senses, [&](const path_handle_t& path_handle) {
         string path_name = graph->get_path_name(path_handle);
         if (!this->ref_paths.count(path_name)) {
+            // This isn't a designated decosntruction reference path.
+            // Note that we allow alt paths through here.
+
             string sample_name = graph->get_sample_name(path_handle);
             // for backward compatibility
             if (sample_name == PathMetadata::NO_SAMPLE_NAME) {

src/gbwtgraph_helper.cpp

@@ -362,14 +362,17 @@ size_t estimate_hash_table_size(const gbwtgraph::GBZ& gbz, bool progress) {
     size_t genome_size = 0;
 
     if (gbz.graph.get_path_count() > 0) {
-        gbz.graph.for_each_path_handle([&](const path_handle_t& path_handle) {
-            gbz.graph.for_each_step_in_path(path_handle, [&](const step_handle_t& step_handle) {
-                handle_t handle = gbz.graph.get_handle_of_step(step_handle);
-                genome_size += gbz.graph.get_length(handle);
-            });
+        gbz.graph.for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](const path_handle_t& path_handle) {
+            std::string path_name = gbz.graph.get_path_name(path_handle);
+            if (!Paths::is_alt(path_name)) {
+                gbz.graph.for_each_step_in_path(path_handle, [&](const step_handle_t& step_handle) {
+                    handle_t handle = gbz.graph.get_handle_of_step(step_handle);
+                    genome_size += gbz.graph.get_length(handle);
+                });
+            }
         });
         if (progress) {
-            std::cerr << "Estimated size based on reference / generic paths: " << genome_size << std::endl;
+            std::cerr << "Estimated size based on reference / non-alt generic paths: " << genome_size << std::endl;
         }
     }
 

src/graph_synchronizer.cpp

@@ -12,10 +12,10 @@ GraphSynchronizer::GraphSynchronizer(VG& graph) : graph(graph) {
     // build a path index after a path has been modified (since we don't keep
     // the ranks up to date internally), we need to build all the indexes up
     // front, even if we're just working on a single path.
-    graph.for_each_path_handle([&](const path_handle_t& path) {
+    graph.for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](const path_handle_t& path) {
         string name = graph.get_path_name(path);
         if (!Paths::is_alt(name)) {
-            // We only care about reference paths.
+            // We only care about reference and non-alt generic paths.
             get_path_index(name);
         }
     });

src/haplotype_indexer.cpp

@@ -30,9 +30,9 @@ HaplotypeIndexer::HaplotypeIndexer() {
 
 std::vector<std::string> HaplotypeIndexer::parse_vcf(const std::string& filename, const PathHandleGraph& graph, const std::string& job_name) const {
 
-    // Parse all non-alt paths.
+    // Parse all non-alt, non-haplotype paths.
     std::vector<path_handle_t> path_handles;
-    graph.for_each_path_handle([&](path_handle_t path_handle) {
+    graph.for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](path_handle_t path_handle) {
         std::string path_name = graph.get_path_name(path_handle);
         if (!Paths::is_alt(path_name)) {
             path_handles.push_back(path_handle);

src/index_registry.cpp

@@ -2606,7 +2606,7 @@ IndexRegistry VGIndexes::get_vg_index_registry() {
         // need a job to do them.
         unordered_set<path_handle_t> broadcast_graph_paths_to_do;
         if (include_named_paths && broadcast_graph) {
-            broadcast_graph->for_each_path_handle([&](const path_handle_t& path_handle) {
+            broadcast_graph->for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](const path_handle_t& path_handle) {
                 // Look at all the paths in advance
                 if (broadcast_graph->is_empty(path_handle) || Paths::is_alt(broadcast_graph->get_path_name(path_handle))) {
                     // Skip empty paths and alt allele paths

src/mcmc_genotyper.cpp

@@ -397,11 +397,14 @@ namespace vg {
         unique_ptr<PhasedGenome> genome(new PhasedGenome(snarls));
         vector<NodeTraversal> haplotype; //will add twice  
 
-        graph.for_each_path_handle([&](const path_handle_t& path){
+        // TODO: This just concatenates all the reference and generic paths.
+        // We're going to need to change this to handle multiple contigs, or
+        // multiple references in the graph, for this to really work.
+        graph.for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](const path_handle_t& path){
         // capture all variables (paths) in scope by reference 
 
             if(!Paths::is_alt(graph.get_path_name(path))) {
-            // If it isn't an alt path, we want to trace it
+            // If it isn't an alt or haplotype path, we want to trace it
 
                 for (handle_t handle : graph.scan_path(path)) {
                 // For each occurrence from start to end

src/rare_variant_simplifier.cpp

@@ -12,11 +12,11 @@ void RareVariantSimplifier::simplify() {
     // This holds the IDs of all the nodes we want to keep around
     unordered_set<id_t> to_keep;
 
-    graph.for_each_path_handle([&](const path_handle_t& path) {
+    graph.for_each_path_of_sense({PathSense::REFERENCE, PathSense::GENERIC}, [&](const path_handle_t& path) {
         // For each path
 
         if (!Paths::is_alt(graph.get_path_name(path))) {
-            // If it isn't an alt path, we want to trace it
+            // If it isn't an alt or haplotype path, we want to trace it
 
             // For each occurrence from start to end
             // Put the ID of the node we are visiting in the to-keep set