diff --git a/src/aln.cpp b/src/aln.cpp index 47571dd8..9496e4ba 100644 --- 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 index 7ff9b51e..300de9e0 100644 --- 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 index b5eb5986..30f5c368 100644 --- 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 index 102c481c..1e7fbb20 100644 --- a/src/nam.cpp +++ b/src/nam.cpp @@ -197,7 +197,8 @@ std::vector merge_hits_into_nams_forward_and_reverse( */ std::tuple> find_nams( const QueryRandstrobeVector &query_randstrobes, - const StrobemerIndex& index + const StrobemerIndex& index, + bool use_mcs ) { std::vector 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> 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> find_nams( std::pair> 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> 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 index 1926bd6d..1486825a 100644 --- a/src/nam.hpp +++ b/src/nam.hpp @@ -37,13 +37,15 @@ struct Nam { std::tuple> find_nams( const QueryRandstrobeVector &query_randstrobes, - const StrobemerIndex& index + const StrobemerIndex& index, + bool use_mcs ); std::pair> 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);