Optimised reverse complement

[jeff-k]

It should be possible to implement faster reverse complementing like jellyfish and the classic biostars post:

• https://github.com/gmarcais/Jellyfish/blob/master/include/jellyfish/mer_dna.hpp
• https://www.biostars.org/p/113640/

[jeff-k]

We'll want a fast reversing of usizes for codecs that don't complement. Codecs with a width that doesn't divide 64 could be problematic.

[werner291]

Hi, I'm exploring the Rust biology ecosystem - I'm very interested in type safety in these kinds of crates (while observing rather frequent use of rather unconstrained byte-strings in many cases) and saw this issue. The algorithm indeed works very nicely:

// Needs a pub `Seq::into_raw_mut(&mut self) -> &mut [usize]` companion
// to `into_raw`; `bv` is `pub(crate)`.
impl Seq<Dna> {
    fn revcomp_fast(&mut self) {
        assert_eq!(self.len() * 2 % usize::BITS as usize, 0, "aligned only");
        let words = self.into_raw_mut();
        words.reverse();
        for w in words.iter_mut() {
            let r = (!*w).reverse_bits();
            *w = ((r & 0x5555_5555_5555_5555usize) << 1)
                | ((r & 0xAAAA_AAAA_AAAA_AAAAusize) >> 1);
        }
    }
}

Against current Seq<Dna>::revcomp():

n=     1024  trait=  0.014 ms  word=  0.000 ms  ≈300x
n=    65536  trait=  1.727 ms  word=  0.003 ms  ≈570x
n=  1048576  trait= 27.538 ms  word=  0.048 ms  ≈570x
n= 16777216  trait=442.695 ms  word=  0.773 ms  ≈570x

~30× of that is the algorithmic gap; the rest is BitSlice::reverse going bit-by-bit on Lsb0.

Two open questions before a PR:

1. The blanket impl<A: Codec> {Reverse,Complement}Mut for Seq<A> blocks specialization on Seq<Dna>. Add Codec::rev_word/comp_word hooks with default decode-loops, override for Dna? Or drop the blankets?
2. Aligned-only here. Unaligned needs a fallback or a post-reverse multi-word shift (~10 lines) — include?

[jeff-k]

wow that's a huge difference! thanks a lot for the benchmarking.

If you're up to it, the benchmarks themselves would be a valuable contribution. I mentioned them here: #10 but haven't elaborated on it much.

This would be great to put in.

Before we start you should know that I've started work on an alternative sequence storage architecture- it preserves the current API but adds a layer of abstraction between the underlying bit sequence (currently provided by bitvec) which handles bitwise operations. So operations like this will be handled at that level. The branch is here: https://github.com/jeff-k/bio-seq/tree/seqstorage but isn't usable yet. However, it should make adding this improvement easier.

I'm going to revise the roadmap, make a dev branch, and make it clear what the plan is. It might be a lot easier to work off of that.