Less stringent alignment and seed filtering

metagraph/CMakeLists.txt

@@ -288,7 +288,7 @@ target_include_directories(caches INTERFACE external-libraries/caches/include)
 
 add_library(eigen INTERFACE)
 target_include_directories(eigen INTERFACE external-libraries/eigen)
-target_compile_options(eigen INTERFACE -Wno-unused-but-set-variable)
+target_compile_options(eigen INTERFACE -Wno-unused-but-set-variable -Wno-unused-parameter)
 
 
 set(Boost_USE_STATIC_LIBS ON)

metagraph/src/cli/config/config.cpp

@@ -518,12 +518,6 @@ Config::Config(int argc, char *argv[]) {
         print_usage_and_exit = true;
     }
 
-    // only the best alignment is used in query
-    // |alignment_num_alternative_paths| must be set to 1
-    if (identity == QUERY && align_sequences
-                          && alignment_num_alternative_paths != 1)
-        print_usage_and_exit = true;
-
     if (identity == ALIGN && infbase.empty())
         print_usage_and_exit = true;
 
@@ -1063,7 +1057,7 @@ if (advanced) {
             fprintf(stderr, "\t   --align-only-forwards \t\t\tdo not align backwards from a seed on basic-mode graphs [off]\n");
             fprintf(stderr, "\t   --align-no-seed-complexity-filter \t\t\t\tdisable the filter for low-complexity seeds. [off]\n");
 }
-            fprintf(stderr, "\t   --align-alternative-alignments \t\tthe number of alternative paths to report per seed [1]\n");
+            fprintf(stderr, "\t   --align-alternative-alignments \t\tthe maximum number of paths to report per seed [inf]\n");
             fprintf(stderr, "\t   --align-chain \t\t\t\tconstruct seed chains before alignment. Useful for long error-prone reads. [off]\n");
             fprintf(stderr, "\t   --align-post-chain \t\t\tperform multiple local alignments and chain them together into a single alignment. Useful for long error-prone reads. [off]\n");
             fprintf(stderr, "\t         \t\t\t\t\t\tA '$' inserted into the reference sequence indicates a jump in the graph.\n");

metagraph/src/cli/config/config.hpp

@@ -124,7 +124,7 @@ class Config {
     int32_t alignment_min_path_score = 0;
     int32_t alignment_xdrop = 27;
 
-    size_t alignment_num_alternative_paths = 1;
+    size_t alignment_num_alternative_paths = std::numeric_limits<size_t>::max();
     size_t alignment_min_seed_length = 19;
     size_t alignment_max_seed_length = std::numeric_limits<size_t>::max();
     size_t alignment_max_num_seeds_per_locus = 1000;

metagraph/src/cli/query.cpp

@@ -1172,9 +1172,6 @@ int query_graph(Config *config) {
 
     std::unique_ptr<align::DBGAlignerConfig> aligner_config;
     if (config->align_sequences) {
-        assert(config->alignment_num_alternative_paths == 1u
-                && "only the best alignment is used in query");
-
         aligner_config.reset(new align::DBGAlignerConfig(
             initialize_aligner_config(*config, *graph)
         ));
@@ -1233,10 +1230,10 @@ Alignment align_sequence(std::string *seq,
         + revised_config.left_end_bonus + revised_config.right_end_bonus;
     auto alignments = aligner.align(*seq);
 
-    assert(alignments.size() <= 1 && "Only the best alignment is needed");
-
     if (alignments.size()) {
+        // TODO: incorporate multiple alignments
         auto &match = alignments[0];
+
         // modify sequence for querying with the best alignment
         if (match.get_offset()) {
             *seq = graph.get_node_sequence(match.get_nodes()[0]).substr(0, match.get_offset())

metagraph/src/cli/server.cpp

@@ -151,13 +151,13 @@ std::string process_align_request(const std::string &received_message,
 
     config.alignment_num_alternative_paths = json.get(
         "max_alternative_alignments",
-        (uint64_t)config.alignment_num_alternative_paths).asInt();
+        (uint64_t)config.alignment_num_alternative_paths).asUInt64();
 
     if (!config.alignment_num_alternative_paths) {
         // TODO: better throw an exception and send an error response to the client
         logger->warn("[Server] Got invalid value of alignment_num_alternative_paths = {}."
-                     " The default value of 1 will be used instead...", config.alignment_num_alternative_paths);
-        config.alignment_num_alternative_paths = 1;
+                     " The default value of inf will be used instead...", config.alignment_num_alternative_paths);
+        config.alignment_num_alternative_paths = std::numeric_limits<size_t>::max();
     }
 
     config.alignment_min_exact_match

metagraph/src/graph/alignment/aligner_aggregator.hpp

@@ -46,7 +46,6 @@ class AlignmentAggregator {
     bool add_alignment(Alignment&& alignment);
 
     score_t get_global_cutoff() const;
-    score_t get_score_cutoff(const Columns &labels) const;
 
     std::vector<Alignment> get_alignments();
 
@@ -59,8 +58,6 @@ class AlignmentAggregator {
     VectorMap<Column, PathQueue> path_queue_;
     PathQueue unlabeled_;
     ValCmp cmp_;
-
-    score_t get_label_cutoff(Column label) const;
 };
 
 // return true if the alignment was added
@@ -148,34 +145,6 @@ ::get_global_cutoff() const -> score_t {
     return cur_max > 0 ? cur_max * config_.rel_score_cutoff : cur_max;
 }
 
-// TODO: define it the same way as in get_global_cutoff()?
-template <class AlignmentCompare>
-inline auto AlignmentAggregator<AlignmentCompare>
-::get_score_cutoff(const Vector<Column> &labels) const -> score_t {
-    assert(labels.size());
-
-    score_t global_min = get_global_cutoff();
-
-    score_t min_score = std::numeric_limits<score_t>::max();
-    for (Column label : labels) {
-        min_score = std::min(min_score, get_label_cutoff(label));
-        if (min_score < global_min)
-            return global_min;
-    }
-    return min_score;
-}
-
-template <class AlignmentCompare>
-inline auto AlignmentAggregator<AlignmentCompare>
-::get_label_cutoff(Column label) const -> score_t {
-    auto find = path_queue_.find(label);
-    return find == path_queue_.end()
-            || find->second.size() < config_.num_alternative_paths
-            || config_.post_chain_alignments
-        ? config_.ninf
-        : find->second.minimum()->get_score();
-}
-
 template <class AlignmentCompare>
 inline std::vector<Alignment> AlignmentAggregator<AlignmentCompare>::get_alignments() {
     // move all alignments to one vector

metagraph/src/graph/alignment/aligner_labeled.cpp

@@ -108,6 +108,11 @@ void LabeledExtender::flush() {
         auto cur_labels = annotation_buffer_.get_labels(table_elem.node);
         assert(cur_labels);
 
+        if (cur_labels->empty()) {
+            assert(table_elem.offset - this->seed_->get_offset() < graph_->get_k());
+            continue;
+        }
+
 #ifndef NDEBUG
         if (table[parent_i].offset >= 0
                 && static_cast<size_t>(table[parent_i].offset) >= graph_->get_k() - 1) {
@@ -227,6 +232,13 @@ ::call_outgoing(node_index node,
             auto next_labels = annotation_buffer_.get_labels(next);
             assert(next_labels);
 
+            if (next_labels->empty()) {
+                assert(next_offset < graph_->get_k());
+                node_labels_.push_back(node_labels_[table_i]);
+                callback(next, c, score);
+                continue;
+            }
+
             Columns intersect_labels;
             std::set_intersection(columns.begin(), columns.end(),
                                   next_labels->begin(), next_labels->end(),
@@ -263,6 +275,13 @@ ::call_outgoing(node_index node,
             std::swap(base_coords, next_coords);
         }
 
+        if (next_labels->empty()) {
+            assert(next_offset < graph_->get_k());
+            node_labels_.push_back(node_labels_[table_i]);
+            callback(next, c, score);
+            continue;
+        }
+
         // check if at least one label has consistent coordinates
         Columns intersect_labels;
 
@@ -454,7 +473,8 @@ ::LabeledAligner(const DeBruijnGraph &graph,
                  const DBGAlignerConfig &config,
                  const Annotator &annotator)
       : DBGAligner<Seeder, Extender, AlignmentCompare>(graph, config),
-        annotation_buffer_(graph, annotator) {
+        annotation_buffer_(graph, annotator),
+        max_seed_length_(this->config_.max_seed_length) {
     // do not use a global xdrop cutoff since we need separate cutoffs for each label
     if (annotation_buffer_.has_coordinates()) {
         logger->trace("Coordinates detected. Enabling seed chaining");
@@ -478,9 +498,12 @@ LabeledAligner<Seeder, Extender, AlignmentCompare>::~LabeledAligner() {
 template <class Seeder, class Extender, class AlignmentCompare>
 auto LabeledAligner<Seeder, Extender, AlignmentCompare>
 ::build_seeders(const std::vector<IDBGAligner::Query> &seq_batch,
-                const std::vector<AlignmentResults> &wrapped_seqs) const -> BatchSeeders {
+                const std::vector<AlignmentResults> &wrapped_seqs,
+                std::vector<std::pair<std::vector<Seed>, std::vector<Seed>>> &discarded_seeds) const -> BatchSeeders {
     BatchSeeders seeders
-        = DBGAligner<Seeder, Extender, AlignmentCompare>::build_seeders(seq_batch, wrapped_seqs);
+        = DBGAligner<Seeder, Extender, AlignmentCompare>::build_seeders(seq_batch, wrapped_seqs, discarded_seeds);
+
+    assert(discarded_seeds.size() == seq_batch.size());
 
     // now we're going to filter the seeds
     logger->trace("Filtering seeds by label. Cur mem usage {} MB", get_curr_RSS() / 1e6);
@@ -532,7 +555,7 @@ ::build_seeders(const std::vector<IDBGAligner::Query> &seq_batch,
         auto &[seeder, seeder_rc] = seeders[i];
         auto &[seeds, num_matching] = counted_seeds[i];
         if (seeds.size()) {
-            num_matching = filter_seeds(seeds);
+            num_matching = filter_seeds(seeds, discarded_seeds[i].first);
             num_seeds_left += seeds.size();
         }
 
@@ -542,7 +565,7 @@ ::build_seeders(const std::vector<IDBGAligner::Query> &seq_batch,
         if (has_rc[i]) {
             auto &[seeds, num_matching] = counted_seeds_rc[i];
             if (seeds.size()) {
-                num_matching = filter_seeds(seeds);
+                num_matching = filter_seeds(seeds, discarded_seeds[i].second);
                 num_seeds_rc_left += seeds.size();
             }
 
@@ -610,14 +633,16 @@ void matched_intersection(AIt a_begin, AIt a_end, BIt a_c_begin,
 
 template <class Seeder, class Extender, class AlignmentCompare>
 size_t LabeledAligner<Seeder, Extender, AlignmentCompare>
-::filter_seeds(std::vector<Seed> &seeds) const {
+::filter_seeds(std::vector<Seed> &seeds,
+               std::vector<Seed> &discarded_seeds) const {
     if (seeds.empty())
         return 0;
 
     size_t query_size = seeds[0].get_clipping() + seeds[0].get_end_clipping()
                             + seeds[0].get_query_view().size();
 
     Columns labels;
+    Columns discarded_labels;
 
     {
         VectorMap<Column, sdsl::bit_vector> label_mapper;
@@ -660,47 +685,88 @@ ::filter_seeds(std::vector<Seed> &seeds) const {
         DEBUG_LOG("Keeping {} / {} labels",
                   std::distance(label_counts.begin(), it), label_counts.size());
 
-        label_counts.erase(it, label_counts.end());
+        labels.reserve(it - label_counts.begin());
+        discarded_labels.reserve(label_counts.end() - it);
 
-        labels.reserve(label_counts.size());
-        for (const auto &[label, count] : label_counts) {
-            labels.push_back(label);
-        }
-    }
-
-    if (labels.empty()) {
-        seeds.clear();
-        return 0;
+        std::transform(label_counts.begin(), it, std::back_inserter(labels),
+                       [&](const auto &a) { return a.first; });
+        std::transform(it, label_counts.end(), std::back_inserter(discarded_labels),
+                       [&](const auto &a) { return a.first; });
     }
 
     std::sort(labels.begin(), labels.end());
+    std::sort(discarded_labels.begin(), discarded_labels.end());
 
     for (size_t j = 0; j < seeds.size(); ++j) {
         Seed &seed = seeds[j];
         const std::vector<node_index> &nodes = seed.get_nodes();
         assert(nodes.size() == 1);
-        if (!seed.label_encoder) {
-            seed.label_columns.clear();
-            auto [fetch_labels, fetch_coords] = annotation_buffer_.get_labels_and_coords(nodes[0]);
-            assert(fetch_labels);
-            if (annotation_buffer_.has_coordinates()) {
-                assert(fetch_coords);
-                matched_intersection(fetch_labels->begin(), fetch_labels->end(),
-                                     fetch_coords->begin(),
-                                     labels.begin(), labels.end(),
-                                     std::back_inserter(seed.label_columns),
-                                     std::back_inserter(seed.label_coordinates));
-                if (seed.get_offset() && seed.label_coordinates.size()) {
-                    for (auto &tuple : seed.label_coordinates) {
+
+        // if a seed already as labels, use them
+        if (seed.label_encoder)
+            continue;
+
+        seed.label_columns.clear();
+        auto [fetch_labels, fetch_coords] = annotation_buffer_.get_labels_and_coords(nodes[0]);
+        assert(fetch_labels);
+        if (annotation_buffer_.has_coordinates()) {
+            Alignment::Columns discarded_columns;
+            Alignment::CoordinateSet discarded_coords;
+            bool added_discarded = false;
+            assert(fetch_coords);
+            assert(seed.label_coordinates.empty());
+            matched_intersection(fetch_labels->begin(), fetch_labels->end(),
+                                 fetch_coords->begin(),
+                                 labels.begin(), labels.end(),
+                                 std::back_inserter(seed.label_columns),
+                                 std::back_inserter(seed.label_coordinates));
+            matched_intersection(fetch_labels->begin(), fetch_labels->end(),
+                                 fetch_coords->begin(),
+                                 discarded_labels.begin(), discarded_labels.end(),
+                                 std::back_inserter(discarded_columns),
+                                 std::back_inserter(discarded_coords));
+
+            if (seed.label_columns.size())
+                seed.label_encoder = &annotation_buffer_.get_annotator().get_label_encoder();
+
+            if (discarded_columns.size()) {
+                added_discarded = true;
+                auto &discarded_seed = discarded_seeds.emplace_back(seed);
+                discarded_seed.label_encoder = &annotation_buffer_.get_annotator().get_label_encoder();
+                std::swap(discarded_seed.label_columns, discarded_columns);
+                std::swap(discarded_seed.label_coordinates, discarded_coords);
+            }
+
+            if (seed.get_offset()) {
+                for (auto &tuple : seed.label_coordinates) {
+                    for (auto &coord : tuple) {
+                        coord += seed.get_offset();
+                    }
+                }
+
+                if (added_discarded) {
+                    for (auto &tuple : discarded_seeds.back().label_coordinates) {
                         for (auto &coord : tuple) {
-                            coord += seed.get_offset();
+                            coord += discarded_seeds.back().get_offset();
                         }
                     }
                 }
-            } else {
-                std::set_intersection(fetch_labels->begin(), fetch_labels->end(),
-                                      labels.begin(), labels.end(),
-                                      std::back_inserter(seed.label_columns));
+            }
+        } else {
+            Alignment::Columns discarded_columns;
+
+            std::set_intersection(fetch_labels->begin(), fetch_labels->end(),
+                                  labels.begin(), labels.end(),
+                                  std::back_inserter(seed.label_columns));
+
+            std::set_intersection(fetch_labels->begin(), fetch_labels->end(),
+                                  discarded_labels.begin(), discarded_labels.end(),
+                                  std::back_inserter(discarded_columns));
+
+            if (discarded_columns.size()) {
+                auto &discarded_seed = discarded_seeds.emplace_back(seed);
+                discarded_seed.label_encoder = &annotation_buffer_.get_annotator().get_label_encoder();
+                std::swap(discarded_seed.label_columns, discarded_columns);
             }
 
             if (seed.label_columns.size())
@@ -718,6 +784,19 @@ ::filter_seeds(std::vector<Seed> &seeds) const {
         return a.get_query_view().size() >= this->config_.min_seed_length;
     }));
 
+    auto end = merge_into_unitig_mums(this->graph_, seeds.data(), seeds.data() + seeds.size(),
+                                      this->config_.min_seed_length, max_seed_length_);
+    seeds.erase(seeds.begin() + (end - seeds.data()), seeds.end());
+    assert(seeds.size());
+
+    if (discarded_seeds.size()) {
+        end = merge_into_unitig_mums(this->graph_, discarded_seeds.data(),
+                                     discarded_seeds.data() + discarded_seeds.size(),
+                                     this->config_.min_seed_length);
+        discarded_seeds.erase(discarded_seeds.begin() + (end - discarded_seeds.data()),
+                              discarded_seeds.end());
+    }
+
     return get_num_char_matches_in_seeds(seeds.begin(), seeds.end());
 }
 

metagraph/src/graph/alignment/aligner_labeled.hpp

@@ -149,10 +149,12 @@ class LabeledAligner : public DBGAligner<Seeder, Extender, AlignmentCompare>, pu
     typedef typename DBGAligner<Seeder, Extender, AlignmentCompare>::BatchSeeders BatchSeeders;
     BatchSeeders
     virtual build_seeders(const std::vector<IDBGAligner::Query> &seq_batch,
-                          const std::vector<AlignmentResults> &wrapped_seqs) const override final;
+                          const std::vector<AlignmentResults> &wrapped_seqs,
+                          std::vector<std::pair<std::vector<Seed>, std::vector<Seed>>> &discarded_seeds) const override final;
 
     // helper for the build_seeders method
-    size_t filter_seeds(std::vector<Seed> &seeds) const;
+    size_t filter_seeds(std::vector<Seed> &seeds,
+                        std::vector<Seed> &discarded_seeds) const;
 };
 
 } // namespace align