Merge remote-tracking branch 'origin/master' into lr-giraffe

Author: adamnovak

Makefile

@@ -717,15 +717,17 @@ $(LIB_DIR)/libtabixpp.a: $(LIB_DIR)/libhts.a $(TABIXPP_DIR)/*.cpp $(TABIXPP_DIR)
 	+echo "Cflags: -I$(CWD)/$(INC_DIR)" >> $(LIB_DIR)/pkgconfig/tabixpp.pc
 	+echo "Libs: -L$(CWD)/$(LIB_DIR) -ltabixpp" >> $(LIB_DIR)/pkgconfig/tabixpp.pc
 
-# Build vcflib. Install the library and headers but not binaries or man pages.
-# We need to build as RelWithDebInfo to avoid vcflib using its own
+# Build vcflib. Build the library and any needed binarise. Install the library
+# and headers but not binaries or man pages. We need to build as RelWithDebInfo
+# to avoid vcflib using its own
 # -march=native, which would conflict with the -march that comes in through
-# CXXFLAGS from the vg Dockerfile.
-# We also need to use the magic path hint to let CMake find Mac OpenMP.
-# We need to use /usr first for CMake search or Ubuntu 22.04 will decide pybind11 is installed in / when actually it is only fully installed in /usr.
+# CXXFLAGS from the vg Dockerfile. We also need to use the magic path hint to
+# let CMake find Mac OpenMP. We need to use /usr first for CMake search or
+# Ubuntu 22.04 will decide pybind11 is installed in / when actually it is only
+# fully installed in /usr.
 $(LIB_DIR)/libvcflib.a: $(LIB_DIR)/libhts.a $(LIB_DIR)/libtabixpp.a $(VCFLIB_DIR)/src/*.cpp $(VCFLIB_DIR)/src/*.hpp $(VCFLIB_DIR)/contrib/*/*.cpp $(VCFLIB_DIR)/contrib/*/*.h
 	+rm -f $(VCFLIB_DIR)/contrib/WFA2-lib/VERSION
-	+cd $(VCFLIB_DIR) && rm -Rf build && mkdir build && cd build && PKG_CONFIG_PATH="$(CWD)/$(LIB_DIR)/pkgconfig:$(PKG_CONFIG_PATH)" cmake -DCMAKE_C_COMPILER="$(CC)" -DCMAKE_CXX_COMPILER="$(CXX)" -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON -DZIG=OFF -DCMAKE_C_FLAGS="$(CFLAGS)" -DCMAKE_CXX_FLAGS="$(CXXFLAGS) ${CPPFLAGS}" -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_INSTALL_PREFIX=$(CWD) -DCMAKE_PREFIX_PATH="/usr;$(OMP_PREFIXES)" .. && cmake --build .
+	+cd $(VCFLIB_DIR) && rm -Rf build && mkdir build && cd build && PKG_CONFIG_PATH="$(CWD)/$(LIB_DIR)/pkgconfig:$(PKG_CONFIG_PATH)" cmake -DCMAKE_C_COMPILER="$(CC)" -DCMAKE_CXX_COMPILER="$(CXX)" -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON -DZIG=OFF -DCMAKE_C_FLAGS="$(CFLAGS)" -DCMAKE_CXX_FLAGS="$(CXXFLAGS) ${CPPFLAGS}" -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_INSTALL_PREFIX=$(CWD) -DCMAKE_PREFIX_PATH="/usr;$(OMP_PREFIXES)" .. && cmake --build . --target vcflib vcf2tsv
 	+cp $(VCFLIB_DIR)/contrib/filevercmp/*.h* $(INC_DIR)
 	+cp $(VCFLIB_DIR)/contrib/fastahack/*.h* $(INC_DIR)
 	+cp $(VCFLIB_DIR)/contrib/smithwaterman/*.h* $(INC_DIR)

deps/libvgio

@@ -119,7 +119,7 @@ pair<vector<pair<string, int64_t>>, unordered_map<string, int64_t>> extract_path
     return make_pair(path_names_and_lengths, subpath_to_length);
 }
 
-vector<tuple<path_handle_t, size_t, size_t>> get_sequence_dictionary(const string& filename, const vector<string>& path_names, const PathPositionHandleGraph& graph) {
+vector<tuple<path_handle_t, size_t, size_t>> get_sequence_dictionary(const string& filename, const vector<string>& path_names, const std::unordered_set<std::string>& reference_samples, const PathPositionHandleGraph& graph) {
 
     // TODO: We assume we're using the one true default path metadata <-> name mapping
 
@@ -282,9 +282,21 @@ vector<tuple<path_handle_t, size_t, size_t>> get_sequence_dictionary(const strin
         // We got no paths in so we need to guess them.
         // We will deduplicate their base names, without subpath info.
         unordered_set<string> base_names;
+
+        // We will also track the distinct reference samples and complain if
+        // there are several, because people don't usually want that.
+        std::unordered_set<std::string> sample_names_in_use;
 
         // When we find a path or subpath we want, we will keep it.
         auto keep_path_or_subpath = [&](const path_handle_t& path) {
+            if (!reference_samples.empty()) {
+                if (!reference_samples.count(graph.get_sample_name(path))) {
+                    // This is not on any acceptable reference, so skip it.
+                    return;
+                }
+            }
+            sample_names_in_use.insert(graph.get_sample_name(path));
+
             string base_name = get_path_base_name(graph, path);
             int64_t base_length = -1;
             if (graph.get_subrange(path) == PathMetadata::NO_SUBRANGE) {
@@ -305,9 +317,35 @@ vector<tuple<path_handle_t, size_t, size_t>> get_sequence_dictionary(const strin
 
         // First look for reference sense paths and their subpaths
         graph.for_each_path_of_sense(PathSense::REFERENCE, keep_path_or_subpath);
+
+        if (sample_names_in_use.size() > 1 && reference_samples.empty()) {
+            // We auto-detected multiple references. Warn the user that they probably don't want this.
+            #pragma omp critical (cerr)
+            {
+                std::cerr << "warning:[vg::get_sequence_dictionary] Using multiple target reference assemblies (";
+                for (auto it = sample_names_in_use.begin(); it != sample_names_in_use.end(); ++it) {
+                    if (it != sample_names_in_use.begin()) {
+                        std::cerr << ", ";
+                    }
+                    std::cerr << *it;
+                }
+                std::cerr << "). Most tools that read SAM/BAM/CRAM will not support this. Consider providing a reference path list or dictionary, or reference sample name." << std::endl;
+            }
+        }
+
+        if (sample_names_in_use.size() < reference_samples.size()) {
+            for (auto& name : reference_samples) {
+                if (!sample_names_in_use.count(name)) {
+                    // TODO: We should learn to promote whole haplotype-sense samples on the fly. For now error out.
+                    cerr << "error:[vg::get_sequence_dictionary] Requested reference assembly " << name << " is not a reference assembly in the graph." << endl;
+                    exit(1);
+                }
+            }
+        }
 
         if (input_names_lengths.empty()) {
             // If none of those exist, try generic sense paths and their subpaths
+            #pragma omp critical (cerr)
             cerr << "warning:[vg::get_sequence_dictionary] No reference-sense paths available in the graph; falling back to generic paths." << endl;
             graph.for_each_path_of_sense(PathSense::GENERIC, [&](const path_handle_t& path) {
                 if (Paths::is_alt(graph.get_path_name(path))) {

doc/wiki

@@ -73,9 +73,11 @@ unique_ptr<AlignmentEmitter> get_alignment_emitter(const string& filename, const
  * those in the file, if any.
  *
  * If the filename is itself an empty string, and no path names are passed,
- * then all reference-sense paths from the graph will be collected in arbitrary
+ * then reference-sense paths from the graph will be collected in arbitrary
  * order. If there are none, all non-alt-allele generic sense paths from the
- * graph will be collected in arbitrary order.
+ * graph will be collected in arbitrary order. If reference_samples is
+ * nonempty, only reference-sense paths in those particular samples will be
+ * used.
  *
  * TODO: Be able to generate the autosomes human-sort, X, Y, MT order typical
  * of references.
@@ -85,7 +87,7 @@ unique_ptr<AlignmentEmitter> get_alignment_emitter(const string& filename, const
  * This information needs to come from the user in order to be correct, but 
  * if it's not specified, it'll be guessed from the graph
  */
-vector<tuple<path_handle_t, size_t, size_t>> get_sequence_dictionary(const string& filename, const vector<string>& path_names, const PathPositionHandleGraph& graph);                                                   
+vector<tuple<path_handle_t, size_t, size_t>> get_sequence_dictionary(const string& filename, const vector<string>& path_names, const std::unordered_set<std::string>& reference_samples, const PathPositionHandleGraph& graph);                                                   
 
 /**
  * Given a list of path handles and size info (from get_sequence_dictionary), return two things:

src/hts_alignment_emitter.cpp

@@ -1082,13 +1082,15 @@ class MinimizerMapper : public AlignerClient {
      */
     static std::pair<size_t, size_t> align_sequence_between(const pos_t& left_anchor, const pos_t& right_anchor, size_t max_path_length, size_t max_gap_length, const HandleGraph* graph, const GSSWAligner* aligner, Alignment& alignment, const std::string* alignment_name = nullptr, size_t max_dp_cells = std::numeric_limits<size_t>::max(), const std::function<size_t(const Alignment&, const HandleGraph&)>& choose_band_padding = algorithms::pad_band_random_walk());
 
+public:
     /**
      * Version of align_sequence_between() that guarantees that you get the
      * same answer (modulo reverse-complementation) no matter whether the
      * sequence and anchors are reverse-complemented or not.
      */
     static std::pair<size_t, size_t> align_sequence_between_consistently(const pos_t& left_anchor, const pos_t& right_anchor, size_t max_path_length, size_t max_gap_length, const HandleGraph* graph, const GSSWAligner* aligner, Alignment& alignment, const std::string* alignment_name = nullptr, size_t max_dp_cells = std::numeric_limits<size_t>::max(), const std::function<size_t(const Alignment&, const HandleGraph&)>& choose_band_padding = algorithms::pad_band_random_walk());
-    
+
+protected:
     /**
      * Produce a WFAAlignment of the given sequence between the given points
      * that will be the same (modulo reverse-complementation) no matter whether

src/hts_alignment_emitter.hpp

@@ -1011,7 +1011,6 @@ vector<Alignment> MinimizerMapper::map_from_chains(Alignment& aln) {
     return mappings;
 }
 
-#define debug
 void MinimizerMapper::do_fragmenting_on_trees(Alignment& aln, const ZipCodeForest& zip_code_forest, 
         const std::vector<Seed>& seeds, const VectorView<MinimizerMapper::Minimizer>& minimizers,
         const vector<algorithms::Anchor>& seed_anchors,
@@ -1618,7 +1617,6 @@ void MinimizerMapper::do_fragmenting_on_trees(Alignment& aln, const ZipCodeFores
     }
 
 }
-#undef debug
 
 void MinimizerMapper::do_chaining_on_fragments(Alignment& aln, const ZipCodeForest& zip_code_forest, 
         const std::vector<Seed>& seeds, const VectorView<MinimizerMapper::Minimizer>& minimizers,
@@ -3587,14 +3585,19 @@ void MinimizerMapper::with_dagified_local_graph(const pos_t& left_anchor, const
     nid_t local_left_anchor_id = 0;
     nid_t local_right_anchor_id = 0;
     for (auto& kv : local_to_base) {
-        if (kv.second == id(left_anchor) && kv.second == id(right_anchor)) {
+        auto& local_id = kv.first;
+        auto& base_id = kv.second;
+        if (base_id == id(left_anchor) && base_id == id(right_anchor)) {
             // The left and right anchors are on the same node, and this is a copy of it.
             // It could be that the anchors face each other, and we extracted one intervening piece of node.
             // In which case we go through this section once.
             if (local_left_anchor_id == 0 && local_right_anchor_id == 0) {
                 // First time through, say we probably cut out the middle piece of a node
-                local_left_anchor_id = kv.first;
-                local_right_anchor_id = kv.first;
+                local_left_anchor_id = local_id;
+                local_right_anchor_id = local_id;
+#ifdef debug
+                std::cerr << "Assume left and right anchors are both node " << local_id << " representing " << base_id << std::endl;
+#endif
             } else {
                 // Or it could be that we have two pieces of the original
                 // shared node represented as separate nodes, because the
@@ -3619,10 +3622,13 @@ void MinimizerMapper::with_dagified_local_graph(const pos_t& left_anchor, const
                 }
                 // Whichever copy has the lower ID is the left one and
                 // whichever copy has the higher ID is the right one.
-                local_left_anchor_id = std::min(local_left_anchor_id, kv.first);
-                local_right_anchor_id = std::max(local_right_anchor_id, kv.second);
+                local_left_anchor_id = std::min(local_left_anchor_id, local_id);
+                local_right_anchor_id = std::max(local_right_anchor_id, local_id);
+#ifdef debug
+                std::cerr << "Second shared anchor copy as " << local_id << " representing " << base_id << "; left is now " << local_left_anchor_id << " and right is " << local_right_anchor_id << std::endl;
+#endif
             }
-        } else if (kv.second == id(left_anchor)) {
+        } else if (base_id == id(left_anchor)) {
             if (local_left_anchor_id != 0) {
                 // We thought we already figured out the start node; there are
                 // too many copies of our start node to find it. 
@@ -3635,8 +3641,11 @@ void MinimizerMapper::with_dagified_local_graph(const pos_t& left_anchor, const
                 local_graph.serialize("crashdump.vg");
                 throw std::runtime_error(ss.str());
             }
-            local_left_anchor_id = kv.first;
-        } else if (kv.second == id(right_anchor)) {
+            local_left_anchor_id = local_id;
+#ifdef debug
+            std::cerr << "Left anchor is " << local_left_anchor_id << std::endl;
+#endif
+        } else if (base_id == id(right_anchor)) {
             if (local_right_anchor_id != 0) {
                 // We thought we already figured out the end node; there are
                 // too many copies of our end node to find it. 
@@ -3649,18 +3658,15 @@ void MinimizerMapper::with_dagified_local_graph(const pos_t& left_anchor, const
                 local_graph.serialize("crashdump.vg");
                 throw std::runtime_error(ss.str());
             }
-            local_right_anchor_id = kv.first;
+            local_right_anchor_id = local_id;
+#ifdef debug
+            std::cerr << "Right anchor is " << local_right_anchor_id << std::endl;
+#endif
         }
         // TODO: Stop early when we found them all.
     }
 
     if (!is_empty(left_anchor) && local_left_anchor_id == 0) {
-        #pragma omp critical (cerr)
-        {
-            for (auto& kv : local_to_base) {
-                std::cerr << "Local ID " << kv.first << " = base graph ID " << kv.second << std::endl;
-            }
-        }
         // Somehow the left anchor didn't come through. Complain.
         std::stringstream ss;
         ss << "Extracted graph of " << local_graph.get_node_count() << " nodes";
@@ -3716,6 +3722,10 @@ void MinimizerMapper::with_dagified_local_graph(const pos_t& left_anchor, const
 
         // And get the node that that orientation of it is in the strand-split graph
         handle_t overlay_handle = split_graph.get_overlay_handle(local_handle);
+
+#ifdef debug
+        std::cerr << "Left anchor " << local_graph.get_id(local_handle) << (local_graph.get_is_reverse(local_handle) ? "-" : "+") << " in local graph is " << split_graph.get_id(overlay_handle) << (split_graph.get_is_reverse(overlay_handle) ? "-" : "+") << " in strand-split graph." << std::endl;
+#endif
 
         // And use that
         bounding_handles.push_back(overlay_handle);
@@ -3742,6 +3752,10 @@ void MinimizerMapper::with_dagified_local_graph(const pos_t& left_anchor, const
         // And get the node that that orientation of it is in the strand-split graph
         // But flip it because we want to dagify going inwards from the right
         handle_t overlay_handle = split_graph.flip(split_graph.get_overlay_handle(local_handle));
+
+#ifdef debug
+        std::cerr << "Right anchor " << local_graph.get_id(local_handle) << (local_graph.get_is_reverse(local_handle) ? "-" : "+") << " in local graph is " << split_graph.get_id(overlay_handle) << (split_graph.get_is_reverse(overlay_handle) ? "-" : "+") << " in strand-split graph." << std::endl;
+#endif
 
         // And use that
         bounding_handles.push_back(overlay_handle);
@@ -3856,20 +3870,20 @@ std::pair<size_t, size_t> MinimizerMapper::align_sequence_between(const pos_t& l
                     // This is a head in the subgraph, and either matches the
                     // left anchor or we don't have any, so keep it.
 #ifdef debug
-                    std::cerr << "Dagified graph node " << dagified_graph.get_id(h) << " " << dagified_graph.get_is_reverse(h) << " is an acceptable source (" << base_coords.first << " " << base_coords.second << ")" << std::endl;
+                    std::cerr << "Dagified graph node " << dagified_graph.get_id(h) << " " << dagified_graph.get_is_reverse(h) << " is an acceptable source (" << base_coords.first << " " << base_coords.second << ") length " << dagified_graph.get_length(h) << std::endl;
 #endif
                 } else if (dagified_graph.get_is_reverse(h) && (is_empty(right_anchor) || h == dagified_graph.flip(right_anchor_handle))) {
                     // Tip is inward reverse, so it's a sink.
                     // This is a tail in the subgraph, and either matches the right
                     // anchor (which faces outwards) or we don't have any, so keep it.
 #ifdef debug
-                    std::cerr << "Dagified graph node " << dagified_graph.get_id(h) << " " << dagified_graph.get_is_reverse(h) << " is an acceptable sink (" << base_coords.first << " " << base_coords.second << ")" << std::endl;
+                    std::cerr << "Dagified graph node " << dagified_graph.get_id(h) << " " << dagified_graph.get_is_reverse(h) << " is an acceptable sink (" << base_coords.first << " " << base_coords.second << ") length " << dagified_graph.get_length(h) << std::endl;
 #endif
                 } else {
                     // This is a wrong orientation or other copy of an anchoring node, or some other tip.
                     // We don't want to keep this handle
 #ifdef debug
-                    std::cerr << "Dagified graph node " << dagified_graph.get_id(h) << " " << dagified_graph.get_is_reverse(h) << " is an unacceptable tip (" << base_coords.first << " " << base_coords.second << ")" << std::endl;
+                    std::cerr << "Dagified graph node " << dagified_graph.get_id(h) << " " << dagified_graph.get_is_reverse(h) << " is an unacceptable tip (" << base_coords.first << " " << base_coords.second << ") length " << dagified_graph.get_length(h) << std::endl;
 #endif
                     nid_t dagified_id = dagified_graph.get_id(h);
                     if (!to_remove_ids.count(dagified_id)) {

src/minimizer_mapper.hpp

@@ -288,7 +288,8 @@ std::pair<string, size_t> PrimerFinder::get_graph_coordinates_from_sequence(cons
     //TODO: These are empty but they could be command line arguments
     string path_file;
     vector<string> path_names;
-    vector<tuple<path_handle_t, size_t, size_t>> sequence_dictionary = get_sequence_dictionary(path_file, path_names, *graph);
+    std::unordered_set<std::string> reference_assembly_names;
+    vector<tuple<path_handle_t, size_t, size_t>> sequence_dictionary = get_sequence_dictionary(path_file, path_names, reference_assembly_names, *graph);
     unordered_set<path_handle_t> reference_paths;
     reference_paths.reserve(sequence_dictionary.size());
     for (auto& entry : sequence_dictionary) {