Make use of mcs optional and disable it by default

ksahlin · Oct 1, 2024 · ae8d35e · ae8d35e
1 parent 6a0b344
commit ae8d35e
Show file tree

Hide file tree

Showing 5 changed files with 16 additions and 10 deletions.
diff --git a/src/aln.cpp b/src/aln.cpp
@@ -1038,7 +1038,7 @@ void align_or_map_paired(
 
         // Find NAMs
         Timer nam_timer;
-        auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index);
+        auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index, map_param.use_mcs);
         statistics.tot_find_nams += nam_timer.duration();
         statistics.n_hits += n_hits;
         details[is_revcomp].nams = nams.size();
@@ -1047,7 +1047,7 @@ void align_or_map_paired(
             Timer rescue_timer;
             if (nams.empty() || nonrepetitive_fraction < 0.7) {
                 int n_rescue_hits;
-                std::tie(n_rescue_hits, nams) = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff);
+                std::tie(n_rescue_hits, nams) = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff, map_param.use_mcs);
                 details[is_revcomp].nam_rescue = true;
                 details[is_revcomp].rescue_nams = nams.size();
                 statistics.n_rescue_hits += n_rescue_hits;
@@ -1181,7 +1181,7 @@ void align_or_map_single(
 
     // Find NAMs
     Timer nam_timer;
-    auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index);
+    auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index, map_param.use_mcs);
     statistics.tot_find_nams += nam_timer.duration();
     statistics.n_hits += n_hits;
     details.nams = nams.size();
@@ -1190,7 +1190,7 @@ void align_or_map_single(
         Timer rescue_timer;
         if (nams.empty() || nonrepetitive_fraction < 0.7) {
             int n_rescue_hits;
-            std::tie(n_rescue_hits, nams) = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff);
+            std::tie(n_rescue_hits, nams) = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff, map_param.use_mcs);
             statistics.n_rescue_hits += n_rescue_hits;
             details.rescue_nams = nams.size();
             details.nam_rescue = true;

diff --git a/src/aln.hpp b/src/aln.hpp
@@ -26,6 +26,7 @@ struct MappingParameters {
     int rescue_level { 2 };
     int max_tries { 20 };
     int rescue_cutoff;
+    bool use_mcs{false};  // multi-context seeds
     OutputFormat output_format {OutputFormat::SAM};
     CigarOps cigar_ops{CigarOps::M};
     bool output_unmapped { true };

diff --git a/src/main.cpp b/src/main.cpp
@@ -196,6 +196,7 @@ int run_strobealign(int argc, char **argv) {
     map_param.dropoff_threshold = opt.dropoff_threshold;
     map_param.rescue_level = opt.rescue_level;
     map_param.max_tries = opt.max_tries;
+    map_param.use_mcs = false; //opt.r < 200;
     map_param.output_format = (
             opt.is_abundance_out ? OutputFormat::Abundance :
             opt.is_sam_out ? OutputFormat::SAM :

diff --git a/src/nam.cpp b/src/nam.cpp
@@ -197,7 +197,8 @@ std::vector<Nam> merge_hits_into_nams_forward_and_reverse(
  */
 std::tuple<float, int, std::vector<Nam>> find_nams(
     const QueryRandstrobeVector &query_randstrobes,
-    const StrobemerIndex& index
+    const StrobemerIndex& index,
+    bool use_mcs
 ) {
     std::vector<PartialSeed> partial_queried; // TODO: is a small set more efficient than linear search in a small vector?
     partial_queried.reserve(10);
@@ -216,7 +217,7 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
             nr_good_hits++;
             add_to_hits_per_ref_full(hits_per_ref[q.is_reverse], q.start, q.end, index, position);
         }
-        else {
+        else if (use_mcs) {
             PartialSeed ph{q.hash >> index.get_aux_len(), q.partial_start, q.is_reverse};
             if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
                 // already queried
@@ -249,7 +250,8 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
 std::pair<int, std::vector<Nam>> find_nams_rescue(
     const QueryRandstrobeVector &query_randstrobes,
     const StrobemerIndex& index,
-    unsigned int rescue_cutoff
+    unsigned int rescue_cutoff,
+    bool use_mcs
 ) {
     struct RescueHit {
         size_t position;
@@ -284,7 +286,7 @@ std::pair<int, std::vector<Nam>> find_nams_rescue(
                 hits_fw.push_back(rh);
             }
         }
-        else {
+        else if (use_mcs) {
             PartialSeed ph = {qr.hash >> index.get_aux_len(), qr.partial_start, qr.is_reverse};
             if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
                 // already queried

diff --git a/src/nam.hpp b/src/nam.hpp
@@ -37,13 +37,15 @@ struct Nam {
 
 std::tuple<float, int, std::vector<Nam>> find_nams(
     const QueryRandstrobeVector &query_randstrobes,
-    const StrobemerIndex& index
+    const StrobemerIndex& index,
+    bool use_mcs
 );
 
 std::pair<int, std::vector<Nam>> find_nams_rescue(
     const QueryRandstrobeVector &query_randstrobes,
     const StrobemerIndex& index,
-    unsigned int rescue_cutoff
+    unsigned int rescue_cutoff,
+    bool use_mcs
 );
 
 std::ostream& operator<<(std::ostream& os, const Nam& nam);