Re-use syncmers of the forward sequence for the revcomp sequence

Author: marcelm

src/index.rs

@@ -545,16 +545,17 @@ impl<'a> StrobemerIndex<'a> {
         if seq.len() < self.parameters.randstrobe.w_max {
             return;
         }
-        let mut syncmer_iter = SyncmerIterator::new(
+        let syncmer_iter = SyncmerIterator::new(
             seq,
             self.parameters.syncmer.k,
             self.parameters.syncmer.s,
             self.parameters.syncmer.t,
         );
+
         let randstrobe_iter =
-            RandstrobeIterator::new(&mut syncmer_iter, &self.parameters.randstrobe);
+            RandstrobeIterator::new(syncmer_iter, self.parameters.randstrobe.clone());
 
-        let mut n= 0;
+        let mut n = 0;
         for (i, randstrobe) in randstrobe_iter.enumerate() {
             n += 1;
             let offset = randstrobe.strobe2_pos - randstrobe.strobe1_pos;

src/mapper.rs

@@ -84,34 +84,53 @@ pub struct QueryRandstrobe {
 /// Generate randstrobes for a query sequence and its reverse complement.
 /// TODO move to strobes.rs?
 pub fn randstrobes_query(seq: &[u8], parameters: &IndexParameters) -> [Vec<QueryRandstrobe>; 2] {
-    let mut randstrobes= [Vec::<QueryRandstrobe>::new(), Vec::<QueryRandstrobe>::new()];
+    let mut randstrobes = [vec![], vec![]];
     if seq.len() < parameters.randstrobe.w_max {
         return randstrobes;
     }
 
-    // TODO
-    // For the reverse complement, we could re-use the syncmers of the forward
+    // Generate syncmers for the forward sequence
+    let mut syncmer_iter = SyncmerIterator::new(seq, parameters.syncmer.k, parameters.syncmer.s, parameters.syncmer.t);
+    let mut syncmers: Vec<_> = syncmer_iter.collect();
+
+    // Generate randstrobes for the forward sequence
+    let randstrobe_iter = RandstrobeIterator::new(syncmers.iter().cloned(), parameters.randstrobe.clone());
+
+    for randstrobe in randstrobe_iter {
+        randstrobes[0].push(
+            QueryRandstrobe {
+                hash: randstrobe.hash,
+                hash_revcomp: randstrobe.hash_revcomp,
+                start: randstrobe.strobe1_pos,
+                end: randstrobe.strobe2_pos + parameters.syncmer.k,
+            }
+        );
+    }
+
+    // For the reverse complement, we can re-use the syncmers of the forward
     // sequence because canonical syncmers are invariant under reverse
     // complementing. Only the coordinates need to be adjusted.
-
-    let seq_rc = reverse_complement(seq);
-    for (s, is_revcomp) in [(seq, false), (&seq_rc, true)] {
-        // Generate randstrobes for the forward sequence
-        let mut syncmer_iter = SyncmerIterator::new(s, parameters.syncmer.k, parameters.syncmer.s, parameters.syncmer.t);
-        let randstrobe_iter = RandstrobeIterator::new(&mut syncmer_iter, &parameters.randstrobe);
-
-        for randstrobe in randstrobe_iter {
-            randstrobes[is_revcomp as usize].push(
-                QueryRandstrobe {
-                    hash: randstrobe.hash,
-                    hash_revcomp: randstrobe.hash_revcomp,
-                    start: randstrobe.strobe1_pos,
-                    end: randstrobe.strobe2_pos + parameters.syncmer.k,
-                }
-            );
-        }
+    syncmers.reverse();
+    for i in 0..syncmers.len() {
+        syncmers[i].position = seq.len() - syncmers[i].position - parameters.syncmer.k;
     }
 
+    // Randstrobes cannot be re-used for the reverse complement:
+    // If in the forward direction, syncmer[i] and syncmer[j] were paired up, it
+    // is not necessarily the case that syncmer[j] is going to be paired with
+    // syncmer[i] in the reverse direction because i is fixed in the forward
+    // direction and j is fixed in the reverse direction.
+    let rc_randstrobe_iter = RandstrobeIterator::new(syncmers.into_iter(), parameters.randstrobe.clone());
+    for randstrobe in rc_randstrobe_iter {
+        randstrobes[1].push(
+            QueryRandstrobe {
+                hash: randstrobe.hash,
+                hash_revcomp: randstrobe.hash_revcomp,
+                start: randstrobe.strobe1_pos,
+                end: randstrobe.strobe2_pos + parameters.syncmer.k,
+            }
+        );
+    }
     randstrobes
 }
 

src/strobes.rs

@@ -57,14 +57,14 @@ impl Randstrobe {
     }
 }
 
-pub struct RandstrobeIterator<'a> {
-    parameters: &'a RandstrobeParameters,
+pub struct RandstrobeIterator<I: Iterator<Item = Syncmer>> {
+    parameters: RandstrobeParameters,
     syncmers: VecDeque<Syncmer>,
-    syncmer_iterator: &'a mut SyncmerIterator<'a>,
+    syncmer_iterator: I,
 }
 
-impl<'a> RandstrobeIterator<'a> {
-    pub fn new(syncmer_iterator: &'a mut SyncmerIterator<'a>, parameters: &'a RandstrobeParameters) -> RandstrobeIterator<'a> {
+impl<I: Iterator<Item = Syncmer>> RandstrobeIterator<I> {
+    pub fn new(syncmer_iterator: I, parameters: RandstrobeParameters) -> RandstrobeIterator<I> {
         RandstrobeIterator {
             parameters,
             syncmers: VecDeque::<Syncmer>::new(),
@@ -73,7 +73,7 @@ impl<'a> RandstrobeIterator<'a> {
     }
 }
 
-impl Iterator for RandstrobeIterator<'_> {
+impl<SI: Iterator<Item = Syncmer>> Iterator for RandstrobeIterator<SI> {
     type Item = Randstrobe;
     fn next(&mut self) -> Option<Self::Item> {
         while self.syncmers.len() <= self.parameters.w_max {
@@ -128,8 +128,8 @@ mod test {
     fn test_randstrobe_iterator() {
         let refseq = read_phix().sequence;
         let parameters = IndexParameters::default_from_read_length(300);
-        let mut syncmer_iter = SyncmerIterator::new(&refseq, parameters.syncmer.k, parameters.syncmer.s, parameters.syncmer.t);
-        let randstrobe_iter = RandstrobeIterator::new(&mut syncmer_iter, &parameters.randstrobe);
+        let syncmer_iter = SyncmerIterator::new(&refseq, parameters.syncmer.k, parameters.syncmer.s, parameters.syncmer.t);
+        let randstrobe_iter = RandstrobeIterator::new(syncmer_iter, parameters.randstrobe.clone());
 
         for randstrobe in randstrobe_iter {
             assert!(randstrobe.hash > 0);
@@ -148,8 +148,8 @@ mod test {
         let syncmer_iter = SyncmerIterator::new(&refseq, parameters.syncmer.k, parameters.syncmer.s, parameters.syncmer.t);
         let syncmer_count = syncmer_iter.count();
 
-        let mut syncmer_iter = SyncmerIterator::new(&refseq, parameters.syncmer.k, parameters.syncmer.s, parameters.syncmer.t);
-        let randstrobe_iter = RandstrobeIterator::new(&mut syncmer_iter, &parameters.randstrobe);
+        let syncmer_iter = SyncmerIterator::new(&refseq, parameters.syncmer.k, parameters.syncmer.s, parameters.syncmer.t);
+        let randstrobe_iter = RandstrobeIterator::new(syncmer_iter, parameters.randstrobe.clone());
         let randstrobe_count = randstrobe_iter.count();
 
         assert_eq!(randstrobe_count, syncmer_count);