Merge remote-tracking branch 'origin/master' into deconstruct

Author: glennhickey

Makefile

@@ -762,9 +762,13 @@ $(LIB_DIR)/libtabixpp.a: $(LIB_DIR)/libhts.a $(TABIXPP_DIR)/*.cpp $(TABIXPP_DIR)
 # 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.
+# We also have to make sure to point at a Python that probably has its headers,
+# because there's no way to turn off vcflib's pybind11 build and pybind11 will
+# try and use the latest installed Python over the default one that probably
+# has headers.
 $(LIB_DIR)/libvcflib%a $(LIB_DIR)/libwfa2%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 . --target vcflib vcf2tsv wfa2_static
+	+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)" -DPYTHON_EXECUTABLE="$(shell which python3)" .. && cmake --build . --target vcflib vcf2tsv wfa2_static
 	+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/libbdsg

@@ -22,7 +22,9 @@ namespace vg {
 
 const char* const BAM_DNA_LOOKUP = "=ACMGRSVTWYHKDBN";
 
-// htslib-based alignment read functions
+// htslib-based alignment read functions.
+// When encountering read records that don't agree with the graph (i.e. go off
+// path ends, etc.), these stop the program and print a useful error message.
 
 int hts_for_each(string& filename, function<void(Alignment&)> lambda);
 int hts_for_each_parallel(string& filename, function<void(Alignment&)> lambda);
@@ -99,11 +101,15 @@ string mapping_string(const string& source, const Mapping& mapping);
 
 void cigar_mapping(const bam1_t *b, Mapping& mapping);
 
+/// Convert a BAM record to an Alignment.
+/// May throw AlignmentEmbeddingError if the BAM record is inconsistent with
+/// the provided graph.
 Alignment bam_to_alignment(const bam1_t *b,
                            const map<string, string>& rg_sample,
                            const map<int, path_handle_t>& tid_path_handle,
                            const bam_hdr_t *bh,
                            const PathPositionHandleGraph* graph);
+/// Convert a BAM record to an Alignment without a graph.
 Alignment bam_to_alignment(const bam1_t *b, const map<string, string>& rg_sample, const map<int, path_handle_t>& tid_path_handle);
 
 // the CIGAR string of the graph alignment
@@ -286,7 +292,11 @@ void reverse_complement_alignment_in_place(Alignment* aln, const function<int64_
 vector<Alignment> reverse_complement_alignments(const vector<Alignment>& alns, const function<int64_t(nid_t)>& node_length);
 int non_match_start(const Alignment& alignment);
 int non_match_end(const Alignment& alignment);
+/// Get the leading softclip from an Alignment, assuming it is coalesced into a
+/// single Edit
 int softclip_start(const Alignment& alignment);
+/// Get the trailing softclip from an Alignment, assuming it is coalesced into a
+/// single Edit
 int softclip_end(const Alignment& alignment);
 int softclip_trim(Alignment& alignment);
 int query_overlap(const Alignment& aln1, const Alignment& aln2);
@@ -400,15 +410,43 @@ struct AlignmentValidity {
 /// node lengths or ids. Result can be used like a bool or inspected for
 /// further details. Does not log anything itself about bad alignments.
 AlignmentValidity alignment_is_valid(const Alignment& aln, const HandleGraph* hgraph, bool check_sequence = false);
-    
+
+/**
+ * Represents a problem when trying to find a path region in a graph as an
+ * Alignment, or when trying to inject a linear CIGAR-based alignment into the
+ * graph along an embedded path.
+ *
+ * This could be a problem like the alignment trying to go out of range on the
+ * embedded linear path/reference, or trying to go across the junction of a
+ * path that isn't really circular.
+ *
+ * We expect the user to be able to cause this with bad inputs, so this
+ * exception should be handled and reported in a helpful way, rather than as a
+ * crash.
+ */
+class AlignmentEmbeddingError : public std::runtime_error {
+public:
+    using std::runtime_error::runtime_error;
+};
+
 /// Make an Alignment corresponding to a subregion of a stored path.
 /// Positions are 0-based, and pos2 is excluded.
 /// Respects path circularity, so pos2 < pos1 is not a problem.
 /// If pos1 == pos2, returns an empty alignment.
+///
+/// Throws AlignmentEmbeddingError if the region goes out of range, or tries to
+/// go across the junction of a non-circular path. Despite taking 0-based
+/// coordinates, error messages will describe 1-based coordinates.
 Alignment target_alignment(const PathPositionHandleGraph* graph, const path_handle_t& path, size_t pos1, size_t pos2,
                            const string& feature, bool is_reverse);
 /// Same as above, but uses the given Mapping, translated directly form a CIGAR string, as a source of edits.
 /// The edits are inserted into the generated Alignment, cut as necessary to fit into the Alignment's Mappings.
+///
+/// Throws AlignmentEmbeddingError if the region goes out of range, or tries to
+/// go across the junction of a non-circular path, or if cigar_mapping
+/// describes edits that are impossible, like matches past the end of the
+/// described region. Despite taking 0-based coordinates, error messages will
+/// describe 1-based coordinates.
 Alignment target_alignment(const PathPositionHandleGraph* graph, const path_handle_t& path, size_t pos1, size_t pos2,
                            const string& feature, bool is_reverse, Mapping& cigar_mapping);
 

deps/libvgio

@@ -5203,6 +5203,10 @@ IndexingPlan IndexRegistry::make_plan(const IndexGroup& end_products) const {
 
         // map dependency priority to requesters
         map<size_t, set<size_t>, greater<size_t>> queue;
+
+        // To help the user, we track index identifiers that, if we had them,
+        // we would have been able to make a deeper plan.
+        set<IndexGroup> missing_index_sets;
 
         auto num_recipes = [&](const IndexGroup& indexes) {
             int64_t num = 0;
@@ -5371,6 +5375,9 @@ IndexingPlan IndexRegistry::make_plan(const IndexGroup& end_products) const {
 #ifdef debug_index_registry
                 cerr << "index " << to_string(index_group) << " cannot be made from existing inputs, need to backtrack" << endl;
 #endif
+
+                // Remember that if we had had this, we could have proceeded.
+                missing_index_sets.insert(index_group);
 
                 // prune to requester and advance to its next recipe, as many times as necessary until
                 // requester has remaining un-tried recipes
@@ -5380,7 +5387,7 @@ IndexingPlan IndexRegistry::make_plan(const IndexGroup& end_products) const {
 
                     if (get<1>(plan_path.back()).empty()) {
                         // this is the product of the plan path, and we're out of recipes for it
-                        throw InsufficientInputException(product, *this);
+                        throw InsufficientInputException(product, *this, missing_index_sets);
                     }
 
                     // remove items off the plan path until we get to the first index that requested
@@ -5745,16 +5752,34 @@ vector<pair<IndexName, vector<IndexName>>> AliasGraph::non_intermediate_aliases(
 }
 
 InsufficientInputException::InsufficientInputException(const IndexName& target,
-                                                       const IndexRegistry& registry) noexcept :
-    runtime_error("Insufficient input to create " + target), target(target), inputs(registry.completed_indexes())
+                                                       const IndexRegistry& registry,
+                                                       const set<IndexGroup>& missing_index_sets) noexcept :
+    runtime_error("Insufficient input to create " + target), target(target), inputs(registry.completed_indexes()), missing_index_sets(missing_index_sets)
 {
     // nothing else to do
     stringstream ss;
-    ss << "Inputs" << endl;
-    for (const auto& input : inputs) {
-        ss << "\t" << input << endl;
+    ss << "Inputs:" << endl;
+    if (inputs.empty()) {
+        ss << "\t<no inputs provided!>" << endl;
+    } else {
+        for (const auto& input : inputs) {
+            ss << "\t" << input << endl;
+        }
+    }
+    ss << "are insufficient to create target index " << target << "." << endl;
+    if (!missing_index_sets.empty()) {
+        ss << "Hint: more progress would be possible if you provided one of these:" << endl;
+        for (auto& index_set : missing_index_sets) {
+            ss << "\t";
+            for (auto name_iter = index_set.begin(); name_iter != index_set.end(); ++name_iter) {
+                if (name_iter != index_set.begin()) {
+                    ss << ", ";
+                }
+                ss << *name_iter;
+            }
+            ss << endl;
+        }
     }
-    ss << "are insufficient to create target index " << target << endl;
     msg = ss.str();
 }
 

deps/sdsl-lite

@@ -452,12 +452,14 @@ class InsufficientInputException : public runtime_error {
 public:
     InsufficientInputException() = delete;
     InsufficientInputException(const IndexName& target,
-                               const IndexRegistry& registry) noexcept;
+                               const IndexRegistry& registry,
+                               const set<IndexGroup>& missing_index_sets = {}) noexcept;
     const char* what() const noexcept;
 private:
     string msg;
     IndexName target;
     vector<IndexName> inputs;
+    set<IndexGroup> missing_index_sets;
 };
 
 

doc/wiki

@@ -432,6 +432,10 @@ class MinimizerMapper : public AlignerClient {
     static constexpr int default_wfa_max_distance = WFAExtender::ErrorModel::default_distance().max;
     int wfa_max_distance = default_wfa_max_distance;
 
+    /// How much should candidate alignment scores be penalized for softclipped bases?
+    static constexpr double default_softclip_penalty = 0.0;
+    double softclip_penalty = default_softclip_penalty;
+
     /// Should alignments be ranked by chain score instead of base-level score?
     static constexpr bool default_sort_by_chain_score = false;
     bool sort_by_chain_score = default_sort_by_chain_score;