From 7ca428d940d7ca8fe9d5648e6d7372961f0e7fca Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Fri, 9 Jun 2023 10:24:58 +0200 Subject: [PATCH 1/3] Asymmetric randstrobe hashes --- src/randstrobes.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/randstrobes.cpp b/src/randstrobes.cpp index 3cfb5a0a..cd78acca 100644 --- a/src/randstrobes.cpp +++ b/src/randstrobes.cpp @@ -45,8 +45,10 @@ static inline syncmer_hash_t syncmer_smer_hash(uint64_t packed) { return xxh64(packed); } -static inline randstrobe_hash_t randstrobe_hash(syncmer_hash_t hash1, syncmer_hash_t hash2) { - return hash1 + hash2; +static inline randstrobe_hash_t randstrobe_hash(syncmer_hash_t syncmer1_hash, syncmer_hash_t syncmer2_hash) { + // alternative: + // (syncmer1_hash << 1) ^ syncmer2_hash + return 2 * syncmer1_hash - syncmer2_hash; } std::ostream& operator<<(std::ostream& os, const Syncmer& syncmer) { From 1960befe6213a675db4aec49653bba127e1e9820 Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Wed, 28 Feb 2024 14:54:49 +0100 Subject: [PATCH 2/3] Randstrobe rescue --- src/aln.cpp | 11 +++++++++ src/randstrobes.cpp | 58 +++++++++++++++++++++++++++++++++++++++++++++ src/randstrobes.hpp | 1 + 3 files changed, 70 insertions(+) diff --git a/src/aln.cpp b/src/aln.cpp index b5105bb1..cf242425 100644 --- a/src/aln.cpp +++ b/src/aln.cpp @@ -972,6 +972,11 @@ void align_or_map_paired( // Find NAMs Timer nam_timer; auto [nonrepetitive_fraction, nams] = find_nams(query_randstrobes, index); + + if (nams.empty()) { + query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters); + std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index); + } statistics.tot_find_nams += nam_timer.duration(); if (map_param.rescue_level > 1) { @@ -1092,6 +1097,12 @@ void align_or_map_single( // Find NAMs Timer nam_timer; auto [nonrepetitive_fraction, nams] = find_nams(query_randstrobes, index); + + if (nams.empty()) { + query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters); + std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index); + } + statistics.tot_find_nams += nam_timer.duration(); if (map_param.rescue_level > 1) { diff --git a/src/randstrobes.cpp b/src/randstrobes.cpp index cd78acca..a2a0ad6c 100644 --- a/src/randstrobes.cpp +++ b/src/randstrobes.cpp @@ -253,3 +253,61 @@ QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexP } return randstrobes; } + +/* + * Generate all possible strobes using a pre-computed vector of syncmers + */ +void randstrobes_rescue(const std::vector& syncmers, std::vector& randstrobes, bool is_rc, RandstrobeParameters parameters, int k) { + const unsigned w_min = parameters.w_min; + const unsigned w_max = parameters.w_max; + const unsigned int max_dist = parameters.max_dist; + if (w_min > w_max) { + throw std::invalid_argument("w_min is greater than w_max"); + } + unsigned int strobe1_index = 0; + + while (strobe1_index + w_min < syncmers.size()) { + unsigned int w_end = std::min(static_cast(strobe1_index + w_max), syncmers.size() - 1); + auto strobe1 = syncmers[strobe1_index]; + auto max_position = strobe1.position + max_dist; + unsigned int w_start = strobe1_index + w_min; + Syncmer strobe2 = strobe1; + + for (auto i = w_start; i <= w_end && syncmers[i].position <= max_position; i++) { + strobe2 = syncmers[i]; + randstrobes.push_back( + QueryRandstrobe{ + randstrobe_hash(strobe1.hash, strobe2.hash), static_cast(strobe1.position), static_cast(strobe2.position + k), is_rc + } + ); + } + strobe1_index++; + } +} + +/* + * Generate randstrobes for a query sequence and its reverse complement. + */ +QueryRandstrobeVector randstrobes_query_rescue(const std::string_view seq, const IndexParameters& parameters) { + QueryRandstrobeVector randstrobes; + if (seq.length() < parameters.randstrobe.w_max) { + return randstrobes; + } + + // Generate syncmers for the forward sequence + auto syncmers = canonical_syncmers(seq, parameters.syncmer); + if (syncmers.empty()) { + return randstrobes; + } + + randstrobes_rescue(syncmers, randstrobes, false, parameters.randstrobe, parameters.syncmer.k); + + std::reverse(syncmers.begin(), syncmers.end()); + for (size_t i = 0; i < syncmers.size(); i++) { + syncmers[i].position = seq.length() - syncmers[i].position - parameters.syncmer.k; + } + + randstrobes_rescue(syncmers, randstrobes, true, parameters.randstrobe, parameters.syncmer.k); + + return randstrobes; +} diff --git a/src/randstrobes.hpp b/src/randstrobes.hpp index 7a79e8ad..15ecfc67 100644 --- a/src/randstrobes.hpp +++ b/src/randstrobes.hpp @@ -69,6 +69,7 @@ std::ostream& operator<<(std::ostream& os, const QueryRandstrobe& randstrobe); using QueryRandstrobeVector = std::vector; QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexParameters& parameters); +QueryRandstrobeVector randstrobes_query_rescue(const std::string_view seq, const IndexParameters& parameters); struct Randstrobe { randstrobe_hash_t hash; From ae4fe29de1575c54f44069e1f54b2b780da8dbe5 Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Wed, 28 Feb 2024 22:02:28 +0100 Subject: [PATCH 3/3] Swap order of rescues --- src/aln.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/aln.cpp b/src/aln.cpp index cf242425..e8594dfa 100644 --- a/src/aln.cpp +++ b/src/aln.cpp @@ -973,10 +973,6 @@ void align_or_map_paired( Timer nam_timer; auto [nonrepetitive_fraction, nams] = find_nams(query_randstrobes, index); - if (nams.empty()) { - query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters); - std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index); - } statistics.tot_find_nams += nam_timer.duration(); if (map_param.rescue_level > 1) { @@ -984,6 +980,10 @@ void align_or_map_paired( if (nams.empty() || nonrepetitive_fraction < 0.7) { nams = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff); details[is_revcomp].nam_rescue = true; + if (nams.empty()) { + query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters); + std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index); + } } statistics.tot_time_rescue += rescue_timer.duration(); } @@ -1098,11 +1098,6 @@ void align_or_map_single( Timer nam_timer; auto [nonrepetitive_fraction, nams] = find_nams(query_randstrobes, index); - if (nams.empty()) { - query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters); - std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index); - } - statistics.tot_find_nams += nam_timer.duration(); if (map_param.rescue_level > 1) { @@ -1110,6 +1105,10 @@ void align_or_map_single( if (nams.empty() || nonrepetitive_fraction < 0.7) { details.nam_rescue = true; nams = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff); + if (nams.empty()) { + query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters); + std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index); + } } statistics.tot_time_rescue += rescue_timer.duration(); }