Merge pull request #472 from ksahlin/base-hash

marcelm · web-flow · commit b2105c67de79 · 2025-01-20T15:04:03.000+01:00
Always pick first syncmer as main
diff --git a/src/index.cpp b/src/index.cpp
@@ -21,7 +21,7 @@
 #include <sstream>
 
 static Logger& logger = Logger::get();
-static const uint32_t STI_FILE_FORMAT_VERSION = 4;
+static const uint32_t STI_FILE_FORMAT_VERSION = 5;
 
 
 namespace {
@@ -309,8 +309,7 @@ void StrobemerIndex::assign_randstrobes(size_t ref_index, size_t offset) {
                 randstrobe.hash,
                 randstrobe.strobe1_pos,
                 static_cast<uint32_t>(ref_index),
-                static_cast<uint8_t>(randstrobe.strobe2_pos - randstrobe.strobe1_pos),
-                randstrobe.first_strobe_is_main
+                static_cast<uint8_t>(randstrobe.strobe2_pos - randstrobe.strobe1_pos)
             };
         }
         chunk.clear();
diff --git a/src/index.hpp b/src/index.hpp
@@ -100,7 +100,7 @@ struct StrobemerIndex {
 
         auto pos = std::lower_bound(randstrobes.begin() + position_start,
                                     randstrobes.begin() + position_end,
-                                    RefRandstrobe{key, 0, 0, 0, 0},
+                                    RefRandstrobe{key, 0, 0, 0},
                                     cmp);
         if ((pos->hash() & hash_mask) == masked_key) return pos - randstrobes.begin();
         return end();
@@ -122,10 +122,6 @@ struct StrobemerIndex {
         }
     }
 
-    bool first_strobe_is_main(bucket_index_t position) const {
-        return randstrobes[position].first_strobe_is_main();
-    }
-    
     bool is_filtered(bucket_index_t position) const {
         return get_hash(position) == get_hash(position + filter_cutoff);
     }
@@ -145,9 +141,6 @@ struct StrobemerIndex {
     std::pair<int, int> strobe_extent_partial(bucket_index_t position) const {
         // Construct the match from the strobe that was selected as the main part of the hash
         int ref_start = get_strobe1_position(position);
-        if (!first_strobe_is_main(position)) {
-            ref_start += strobe2_offset(position);
-        }
         return {ref_start, ref_start + k()};
     }
 
@@ -205,7 +198,7 @@ struct StrobemerIndex {
 
         auto pos = std::upper_bound(randstrobes.begin() + position,
                                     randstrobes.begin() + position_end,
-                                    RefRandstrobe{key, 0, 0, 0, 0},
+                                    RefRandstrobe{key, 0, 0, 0},
                                     cmp);
         return (pos - randstrobes.begin() - 1) - position + 1;
     }
diff --git a/src/nam.cpp b/src/nam.cpp
@@ -13,19 +13,6 @@ bool operator==(const Match& lhs, const Match& rhs) {
     return (lhs.query_start == rhs.query_start) && (lhs.query_end == rhs.query_end) && (lhs.ref_start == rhs.ref_start) && (lhs.ref_end == rhs.ref_end);
 }
 
-/*
- * A partial hit is a hit where not the full randstrobe hash could be found in
- * the index but only the "main" hash (only the first aux_len bits).
- */
-struct PartialHit {
-    randstrobe_hash_t hash;
-    unsigned int start;  // position in strobemer index
-    bool is_reverse;
-    bool operator==(const PartialHit& rhs) const {
-        return (hash == rhs.hash) && (start == rhs.start) && (is_reverse == rhs.is_reverse);
-    }
-};
-
 inline void add_to_matches_map_full(
     robin_hood::unordered_map<unsigned int, std::vector<Match>>& matches_map,
     int query_start,
@@ -47,25 +34,6 @@ inline void add_to_matches_map_full(
     }
 }
 
-/*
-This function often produces the same Match multiple times.
-This happens when multiple randstrobes involve a syncmer that has a low hash
-value (and which therefore gets chosen as main hash multiple times).
-
-For example, assume we have syncmers starting at positions 14, 21, 39, 51
-and randstrobes like this:
-
-14, 21, 39, 51
-14------39
-    21--39
-        39--51
-
-If 39 now has the lowest hash value, it will be chosen as the main hash for all
-three randstrobes. Thus, when looking up that main hash, all three randstrobes
-are found. After the start coordinate of the randstrobe is converted to the
-start coordinate of the single syncmer (using strobe_extent_partial()), they all
-result in a hit that starts at 39 (and has length k).
-*/
 inline void add_to_matches_map_partial(
     robin_hood::unordered_map<unsigned int, std::vector<Match>>& matches_map,
     int query_start,
@@ -223,10 +191,6 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
     const StrobemerIndex& index,
     bool use_mcs
 ) {
-    std::vector<PartialHit> partial_queried; // TODO: is a small set more efficient than linear search in a small vector?
-    if (use_mcs) {
-        partial_queried.reserve(10);
-    }
     std::array<robin_hood::unordered_map<unsigned int, std::vector<Match>>, 2> matches_map;
     matches_map[0].reserve(100);
     matches_map[1].reserve(100);
@@ -243,24 +207,18 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
             add_to_matches_map_full(matches_map[q.is_reverse], q.start, q.end, index, position);
         }
         else if (use_mcs) {
-            PartialHit ph{q.hash & index.get_main_hash_mask(), q.partial_start, q.is_reverse};
-            if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
-                // already queried
-                continue;
-            }
             size_t partial_pos = index.find_partial(q.hash);
             if (partial_pos != index.end()) {
                 total_hits++;
                 if (index.is_partial_filtered(partial_pos)) {
-                    partial_queried.push_back(ph);
                     continue;
                 }
                 nr_good_hits++;
-                add_to_matches_map_partial(matches_map[q.is_reverse], q.partial_start, q.partial_end, index, partial_pos);
+                add_to_matches_map_partial(matches_map[q.is_reverse], q.start, q.start + index.k(), index, partial_pos);
             }
-            partial_queried.push_back(ph);
         }
     }
+
     float nonrepetitive_fraction = total_hits > 0 ? ((float) nr_good_hits) / ((float) total_hits) : 1.0;
     auto nams = merge_matches_into_nams_forward_and_reverse(matches_map, index.k(), use_mcs);
     return {nonrepetitive_fraction, nr_good_hits, nams};
@@ -290,8 +248,6 @@ std::pair<int, std::vector<Nam>> find_nams_rescue(
                 < std::tie(rhs.count, rhs.query_start, rhs.query_end);
         }
     };
-    std::vector<PartialHit> partial_queried;  // TODO: is a small set more efficient than linear search in a small vector?
-    partial_queried.reserve(10);
     std::array<robin_hood::unordered_map<unsigned int, std::vector<Match>>, 2> matches_map;
     std::vector<RescueHit> hits_fw;
     std::vector<RescueHit> hits_rc;
@@ -312,22 +268,16 @@ std::pair<int, std::vector<Nam>> find_nams_rescue(
             }
         }
         else if (use_mcs) {
-            PartialHit ph = {qr.hash & index.get_main_hash_mask(), qr.partial_start, qr.is_reverse};
-            if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
-                // already queried
-                continue;
-            }
             size_t partial_pos = index.find_partial(qr.hash);
             if (partial_pos != index.end()) {
                 unsigned int partial_count = index.get_count_partial(partial_pos);
-                RescueHit rh{partial_pos, partial_count, qr.partial_start, qr.partial_end, true};
+                RescueHit rh{partial_pos, partial_count, qr.start, qr.start + index.k(), true};
                 if (qr.is_reverse){
                     hits_rc.push_back(rh);
                 } else {
                     hits_fw.push_back(rh);
                 }
             }
-            partial_queried.push_back(ph);
         }
     }
 
@@ -356,6 +306,6 @@ std::pair<int, std::vector<Nam>> find_nams_rescue(
 }
 
 std::ostream& operator<<(std::ostream& os, const Nam& n) {
-    os << "Nam(ref_id=" << n.ref_id << ", query: " << n.query_start << ".." << n.query_end << ", ref: " << n.ref_start << ".." << n.ref_end << ", score=" << n.score << ")";
+    os << "Nam(ref_id=" << n.ref_id << ", query: " << n.query_start << ".." << n.query_end << ", ref: " << n.ref_start << ".." << n.ref_end << ", rc=" << static_cast<int>(n.is_rc) << ", score=" << n.score << ")";
     return os;
 }
diff --git a/src/randstrobes.cpp b/src/randstrobes.cpp
@@ -49,8 +49,7 @@ static inline syncmer_hash_t syncmer_smer_hash(uint64_t packed) {
  * This function combines two individual syncmer hashes into a single hash
  * for the randstrobe.
  *
- * The syncmer with the smaller hash is designated as the "main", the other is
- * the "auxiliary".
+ * One syncmer is designated as the "main", the other is the "auxiliary".
  * The combined hash is obtained by setting the top bits to the bits of
  * the main hash and the bottom bits to the bits of the auxiliary
  * hash. Since entries in the index are sorted by randstrobe hash, this allows
@@ -59,10 +58,6 @@ static inline syncmer_hash_t syncmer_smer_hash(uint64_t packed) {
 static inline randstrobe_hash_t randstrobe_hash(
     syncmer_hash_t hash1, syncmer_hash_t hash2, randstrobe_hash_t main_hash_mask
 ) {
-    // Make the function symmetric
-    if (hash1 > hash2) {
-        std::swap(hash1, hash2);
-    }
     return ((hash1 & main_hash_mask) | (hash2 & ~main_hash_mask)) & RANDSTROBE_HASH_MASK;
 }
 
@@ -146,7 +141,7 @@ std::vector<Syncmer> canonical_syncmers(
 
 std::ostream& operator<<(std::ostream& os, const Randstrobe& randstrobe) {
     os << "Randstrobe(hash=" << randstrobe.hash << ", strobe1_pos=" << randstrobe.strobe1_pos << ", strobe2_pos="
-       << randstrobe.strobe2_pos << ", first_strobe_is_main=" << randstrobe.first_strobe_is_main << ")";
+       << randstrobe.strobe2_pos << ")";
     return os;
 }
 
@@ -161,12 +156,10 @@ std::ostream& operator<<(std::ostream& os, const QueryRandstrobe& randstrobe) {
 }
 
 Randstrobe make_randstrobe(Syncmer strobe1, Syncmer strobe2, randstrobe_hash_t main_hash_mask) {
-    bool first_strobe_is_main = strobe1.hash < strobe2.hash;
     return Randstrobe{
         randstrobe_hash(strobe1.hash, strobe2.hash, main_hash_mask),
         static_cast<uint32_t>(strobe1.position),
-        static_cast<uint32_t>(strobe2.position),
-        first_strobe_is_main
+        static_cast<uint32_t>(strobe2.position)
     };
 }
 
@@ -246,11 +239,9 @@ QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexP
     RandstrobeIterator randstrobe_fwd_iter{syncmers, parameters.randstrobe};
     while (randstrobe_fwd_iter.has_next()) {
         auto randstrobe = randstrobe_fwd_iter.next();
-        const unsigned int partial_start = randstrobe.first_strobe_is_main ? randstrobe.strobe1_pos : randstrobe.strobe2_pos;
         randstrobes.push_back(
             QueryRandstrobe {
-                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k,
-                partial_start, partial_start + parameters.syncmer.k,  false
+                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k, false
             }
         );
     }
@@ -271,12 +262,9 @@ QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexP
     RandstrobeIterator randstrobe_rc_iter{syncmers, parameters.randstrobe};
     while (randstrobe_rc_iter.has_next()) {
         auto randstrobe = randstrobe_rc_iter.next();
-        bool first_strobe_is_main = randstrobe.first_strobe_is_main;
-        const unsigned int partial_start = first_strobe_is_main ? randstrobe.strobe1_pos : randstrobe.strobe2_pos;
         randstrobes.push_back(
             QueryRandstrobe {
-                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k,
-                partial_start, partial_start + parameters.syncmer.k, true
+                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k, true
             }
         );
     }
diff --git a/src/randstrobes.hpp b/src/randstrobes.hpp
@@ -17,20 +17,20 @@
 using syncmer_hash_t = uint64_t;
 using randstrobe_hash_t = uint64_t;
 
-static constexpr uint64_t RANDSTROBE_HASH_MASK = 0xFFFFFFFFFFFFFE00;
+static constexpr uint64_t RANDSTROBE_HASH_MASK = 0xFFFFFFFFFFFFFF00;
 
 struct RefRandstrobe {
 private:
-    // packed representation of hash, offset and first_strobe_is_main
+    // packed representation of hash and offset
     randstrobe_hash_t m_hash_offset_flag;
     uint32_t m_position;
     uint32_t m_ref_index;
 
 public:
     RefRandstrobe() : m_hash_offset_flag(0), m_position(0), m_ref_index(0) { }
 
-    RefRandstrobe(randstrobe_hash_t hash, uint32_t position, uint32_t ref_index, uint8_t offset, bool first_strobe_is_main)
-        : m_hash_offset_flag((hash & RANDSTROBE_HASH_MASK) | (offset << 1) | first_strobe_is_main)
+    RefRandstrobe(randstrobe_hash_t hash, uint32_t position, uint32_t ref_index, uint8_t offset)
+        : m_hash_offset_flag((hash & RANDSTROBE_HASH_MASK) | offset)
         , m_position(position)
         , m_ref_index(ref_index)
     { }
@@ -45,16 +45,12 @@ struct RefRandstrobe {
         return lhs < rhs;
     }
 
-    bool first_strobe_is_main() const {
-        return m_hash_offset_flag & 1;
-    }
-
     unsigned reference_index() const {
         return m_ref_index;
     }
 
     unsigned strobe2_offset() const {
-        return (m_hash_offset_flag >> 1) & 0xff;
+        return m_hash_offset_flag & 0xff;
     }
 
     randstrobe_hash_t hash() const {
@@ -72,11 +68,6 @@ struct QueryRandstrobe {
     randstrobe_hash_t hash;
     unsigned int start;
     unsigned int end;
-    /* Start and end of the main syncmer (relevant if the randstrobe couldn’t
-     * be found in the index and we fall back to a partial hit)
-     */
-    unsigned int partial_start;
-    unsigned int partial_end;
     bool is_reverse;
 };
 
@@ -90,7 +81,6 @@ struct Randstrobe {
     randstrobe_hash_t hash;
     unsigned int strobe1_pos;
     unsigned int strobe2_pos;
-    bool first_strobe_is_main;
 
     bool operator==(const Randstrobe& other) const {
         return hash == other.hash && strobe1_pos == other.strobe1_pos && strobe2_pos == other.strobe2_pos;
@@ -188,7 +178,7 @@ class RandstrobeGenerator {
     { }
 
     Randstrobe next();
-    Randstrobe end() const { return Randstrobe{0, 0, 0, false}; }
+    Randstrobe end() const { return Randstrobe{0, 0, 0}; }
 
 private:
     SyncmerIterator syncmer_iterator;
diff --git a/tests/baseline-commit.txt b/tests/baseline-commit.txt
@@ -1 +1 @@
-92c6105d71549e19e8fb263975086c8ffc77e9c8
+7f5ac330f68afbe32e0cde4a32bbb786cdd4d01e
diff --git a/tests/compare-baseline.sh b/tests/compare-baseline.sh
@@ -14,14 +14,18 @@ python3 -c 'import pysam'
 
 ends="pe"
 threads=4
-while getopts "st:" opt; do
+mcs=0
+while getopts "st:m" opt; do
   case "${opt}" in
     t)
       threads=$OPTARG
       ;;
     s)
       ends=se  # single-end reads
       ;;
+    m)
+      mcs=1
+      ;;
     \?)
       exit 1
       ;;
@@ -39,10 +43,15 @@ tests/download.sh
 
 baseline_commit=$(< tests/baseline-commit.txt)
 
-baseline_bam=baseline/bam/${baseline_commit}.${ends}.bam
 baseline_binary=baseline/strobealign-${baseline_commit}
 cmake_options=-DCMAKE_BUILD_TYPE=RelWithDebInfo
+extra_ext=""
 strobealign_options="-t ${threads}"
+if [[ ${mcs} == 1 ]]; then
+  extra_ext=".mcs"
+  strobealign_options="${strobealign_options} --mcs"
+fi
+baseline_bam=baseline/bam/${baseline_commit}.${ends}${extra_ext}.bam
 
 # Generate the baseline BAM if necessary
 mkdir -p baseline/bam
diff --git a/tests/phix.mcs.se.paf b/tests/phix.mcs.se.paf
@@ -1,11 +1,11 @@
-SRR1377138.32	301	2	293	+	NC_001422.1	5386	1434	1725	47	291	255
-SRR1377138.33	301	2	297	+	NC_001422.1	5386	3818	4113	47	295	255
-SRR1377138.34	301	33	299	-	NC_001422.1	5386	844	1110	40	266	255
-SRR1377138.35	301	5	298	-	NC_001422.1	5386	4041	4334	47	293	255
-SRR1377138.36	301	3	301	+	NC_001422.1	5386	4997	5295	51	298	255
-SRR1377138.37	301	2	299	-	NC_001422.1	5386	794	1091	42	297	255
-SRR1377138.38	301	32	284	-	NC_001422.1	5386	4971	5223	42	252	255
-SRR1377138.39/1	301	1	295	-	NC_001422.1	5386	1791	2085	48	294	255
-SRR1377138.40	301	4	293	-	NC_001422.1	5386	3020	3309	44	289	255
-rescuable.42	301	4	293	-	NC_001422.1	5386	3020	3309	44	289	255
-not.rescuable	301	4	293	-	NC_001422.1	5386	3020	3309	44	289	255
+SRR1377138.32	301	2	293	+	NC_001422.1	5386	1434	1725	49	291	255
+SRR1377138.33	301	2	297	+	NC_001422.1	5386	3818	4113	51	295	255
+SRR1377138.34	301	33	299	-	NC_001422.1	5386	844	1110	43	266	255
+SRR1377138.35	301	5	298	-	NC_001422.1	5386	4041	4334	49	293	255
+SRR1377138.36	301	3	301	+	NC_001422.1	5386	4997	5295	53	298	255
+SRR1377138.37	301	2	299	-	NC_001422.1	5386	794	1091	43	297	255
+SRR1377138.38	301	32	284	-	NC_001422.1	5386	4971	5223	44	252	255
+SRR1377138.39/1	301	1	295	-	NC_001422.1	5386	1791	2085	50	294	255
+SRR1377138.40	301	4	293	-	NC_001422.1	5386	3020	3309	48	289	255
+rescuable.42	301	4	293	-	NC_001422.1	5386	3020	3309	48	289	255
+not.rescuable	301	4	293	-	NC_001422.1	5386	3020	3309	48	289	255
diff --git a/tests/samdiff.py b/tests/samdiff.py
diff --git a/tests/test_randstrobes.cpp b/tests/test_randstrobes.cpp

Original file line number	Diff line number	Diff line change
`@@ -1 +1 @@`
`1`		`-92c6105d71549e19e8fb263975086c8ffc77e9c8`
	`1`	`+7f5ac330f68afbe32e0cde4a32bbb786cdd4d01e`