seq_hash/
lib.rs

1//! A crate for streaming hashing of k-mers via `KmerHasher`.
2//!
3//! This builds on [`packed_seq`] and is used by e.g. [`simd_minimizers`].
4//!
5//! The default [`NtHasher`] is canonical.
6//! If that's not needed, [`NtHasher<false>`] will be slightly faster.
7//! For non-DNA sequences with >2-bit alphabets, use [`MulHasher`] instead.
8//!
9//! Note that [`KmerHasher`] objects need `k` on their construction, so that they can precompute required constants.
10//! Prefer reusing the same [`KmerHasher`].
11//!
12//! This crate also includes [`AntiLexHasher`], see [this blogpost](https://curiouscoding.nl/posts/practical-minimizers/).
13//!
14//! ## Typical usage
15//!
16//! Construct a default [`NtHasher`] via `let hasher = <NtHasher>::new(k)`.
17//! Then call either `hasher.hash_kmers_simd(seq, context)`,
18//! or use the underlying 'mapper' via `hasher.in_out_mapper_simd(seq)`.
19//! ```
20//! use packed_seq::{AsciiSeqVec, PackedSeqVec, SeqVec};
21//! use seq_hash::{KmerHasher, NtHasher};
22//! let k = 3;
23//!
24//! // Default `NtHasher` is canonical.
25//! let hasher = <NtHasher>::new(k);
26//! let kmer = PackedSeqVec::from_ascii(b"ACG");
27//! let kmer_rc = PackedSeqVec::from_ascii(b"CGT");
28//! // Normally, prefer `hash_kmers_simd` over `hash_seq`.
29//! assert_eq!(
30//!     hasher.hash_seq(kmer.as_slice()),
31//!     hasher.hash_seq(kmer_rc.as_slice())
32//! );
33//!
34//! let fwd_hasher = NtHasher::<false>::new(k);
35//! assert_ne!(
36//!     fwd_hasher.hash_seq(kmer.as_slice()),
37//!     fwd_hasher.hash_seq(kmer_rc.as_slice())
38//! );
39//!
40//! let seq = b"ACGGCAGCGCATATGTAGT";
41//! let ascii_seq = AsciiSeqVec::from_ascii(seq);
42//! let packed_seq = PackedSeqVec::from_ascii(seq);
43//!
44//! // hasher.hash_kmers_scalar(seq.as_slice()); // Panics since `NtHasher` does not support ASCII.
45//! let hashes_1: Vec<_> = hasher.hash_kmers_scalar(ascii_seq.as_slice()).collect();
46//! let hashes_2: Vec<_> = hasher.hash_kmers_scalar(packed_seq.as_slice()).collect();
47//! // Hashes are equal for [`packed_seq::AsciiSeq`] and [`packed_seq::PackedSeq`].
48//! assert_eq!(hashes_1, hashes_2);
49//! assert_eq!(hashes_1.len(), seq.len() - (k-1));
50//!
51//! // Consider a 'context' of a single kmer.
52//! let hashes_3: Vec<_> = hasher.hash_kmers_simd(ascii_seq.as_slice(), 1).collect();
53//! let hashes_4: Vec<_> = hasher.hash_kmers_simd(packed_seq.as_slice(), 1).collect();
54//! assert_eq!(hashes_1, hashes_3);
55//! assert_eq!(hashes_1, hashes_4);
56//! ```
57
58mod anti_lex;
59mod intrinsics;
60mod nthash;
61#[cfg(test)]
62mod test;
63
64pub use anti_lex::AntiLexHasher;
65pub use nthash::{MulHasher, NtHasher};
66
67/// Re-export of the `packed-seq` crate.
68pub use packed_seq;
69
70use packed_seq::{BitSeq, ChunkIt, Delay, PaddedIt, Seq};
71use std::iter::{repeat, zip};
72
73type S = wide::u32x8;
74
75/// A hasher that can hash all k-mers in a string.
76///
77/// Note that a `KmerHasher` must be initialized with a specific `k`,
78/// so that it can precompute associated constants.
79pub trait KmerHasher {
80    /// True when the hash function is invariant under reverse-complement.
81    const CANONICAL: bool;
82
83    fn new(k: usize) -> Self;
84
85    /// Helper function returning [`Self::CANONICAL`].
86    #[inline(always)]
87    fn is_canonical(&self) -> bool {
88        Self::CANONICAL
89    }
90
91    /// The value of `k` for this hasher.
92    fn k(&self) -> usize;
93
94    /// The delay of the 'out' character passed to the `in_out_mapper` functions.
95    /// Defaults to `k-1`.
96    #[inline(always)]
97    fn delay(&self) -> Delay {
98        Delay(self.k() - 1)
99    }
100
101    /// A scalar mapper function that should be called with each `(in, out)` base.
102    ///
103    /// The delay should be [`Self::delay()`]. The first `delay` calls should have `out=0`.
104    /// `seq` is only used to ensure that the hasher can handle the underlying alphabet.
105    fn in_out_mapper_scalar<'s>(&self, seq: impl Seq<'s>) -> impl FnMut((u8, u8)) -> u32;
106    /// A SIMD mapper function that should be called with a `(in, out)` base per lane.
107    ///
108    /// The delay should be [`Self::delay()`]. The first `delay` calls should have `out=u32x8::splat(0)`.
109    /// `seq` is only used to ensure that the hasher can handle the underlying alphabet.
110    fn in_out_mapper_simd<'s>(&self, seq: impl Seq<'s>) -> impl FnMut((S, S)) -> S;
111
112    fn in_out_mapper_ambiguous_scalar<'s>(
113        &self,
114        seq: impl Seq<'s>,
115        ambiguous: &BitSeq<'s>,
116    ) -> impl FnMut((u8, u8)) -> u32 {
117        let mut mapper = self.in_out_mapper_scalar(seq);
118        let mut ambi = ambiguous.iter_kmer_ambiguity(self.k());
119        let k = self.k();
120        let mut i = 0;
121        move |(a, r)| {
122            let hash = mapper((a, r));
123            let ambi = if i > k - 1 {
124                ambi.next().unwrap()
125            } else {
126                false
127            };
128            i += 1;
129            if ambi { u32::MAX } else { hash }
130        }
131    }
132
133    fn in_out_mapper_ambiguous_simd<'s>(
134        &self,
135        seq: impl Seq<'s>,
136        ambiguous: &BitSeq<'s>,
137        context: usize,
138    ) -> impl FnMut((S, S)) -> S {
139        let mut mapper = self.in_out_mapper_simd(seq);
140        let mut ambi = ambiguous.par_iter_kmer_ambiguity(self.k(), context, 0);
141        move |(a, r)| {
142            let hash = mapper((a, r));
143            let ambi = ambi.it.next().unwrap();
144            ambi.blend(S::MAX, hash)
145        }
146    }
147
148    /// A scalar iterator over all k-mer hashes in `seq`.
149    #[inline(always)]
150    fn hash_kmers_scalar<'s>(&self, seq: impl Seq<'s>) -> impl ExactSizeIterator<Item = u32> {
151        let k = self.k();
152        let delay = self.delay();
153        let mut add = seq.iter_bp();
154        let mut remove = seq.iter_bp();
155        let mut mapper = self.in_out_mapper_scalar(seq);
156        zip(add.by_ref().take(delay.0), repeat(0)).for_each(|a| {
157            mapper(a);
158        });
159        zip(add.by_ref(), remove.by_ref())
160            .take(k - 1 - delay.0)
161            .for_each(|a| {
162                mapper(a);
163            });
164        zip(add, remove).map(mapper)
165    }
166
167    /// A SIMD-parallel iterator over all k-mer hashes in `seq`.
168    #[inline(always)]
169    fn hash_kmers_simd<'s>(&self, seq: impl Seq<'s>, context: usize) -> PaddedIt<impl ChunkIt<S>> {
170        let k = self.k();
171        let delay = self.delay();
172        seq.par_iter_bp_delayed(context + k - 1, delay)
173            .map(self.in_out_mapper_simd(seq))
174            .advance(k - 1)
175    }
176
177    /// An iterator over all k-mer hashes in `seq`.
178    /// Ambiguous kmers get hash `u32::MAX`.
179    #[inline(always)]
180    fn hash_valid_kmers_scalar<'s>(
181        &self,
182        seq: impl Seq<'s>,
183        ambiguous: &BitSeq<'s>,
184    ) -> impl ExactSizeIterator<Item = u32> {
185        let k = self.k();
186        let delay = self.delay();
187        assert!(delay.0 < k);
188
189        let mut mapper = self.in_out_mapper_scalar(seq);
190
191        let mut a = seq.iter_bp();
192        let mut r = seq.iter_bp();
193
194        a.by_ref().take(delay.0).for_each(|a| {
195            mapper((a, 0));
196        });
197
198        zip(a.by_ref(), r.by_ref())
199            .take((k - 1) - delay.0)
200            .for_each(|(a, r)| {
201                mapper((a, r));
202            });
203
204        zip(zip(a, r), ambiguous.iter_kmer_ambiguity(k)).map(move |(ar, ambi)| {
205            let hash = mapper(ar);
206            if ambi { u32::MAX } else { hash }
207        })
208    }
209
210    /// A SIMD-parallel iterator over all k-mer hashes in `seq`.
211    /// Ambiguous kmers get hash `u32::MAX`.
212    #[inline(always)]
213    fn hash_valid_kmers_simd<'s>(
214        &self,
215        seq: impl Seq<'s>,
216        ambiguous: &BitSeq<'s>,
217        context: usize,
218    ) -> PaddedIt<impl ChunkIt<S>> {
219        let k = self.k();
220        let delay = self.delay();
221        seq.par_iter_bp_delayed(context + k - 1, delay)
222            .map(self.in_out_mapper_simd(seq))
223            .zip(ambiguous.par_iter_kmer_ambiguity(k, context + k - 1, 0))
224            .map(|(hash, is_ambiguous)| is_ambiguous.blend(S::MAX, hash))
225            .advance(k - 1)
226    }
227
228    /// Hash a sequence one character at a time. Ignores `k`.
229    ///
230    /// `seq` is only used to ensure that the hasher can handle the underlying alphabet.
231    fn mapper<'s>(&self, seq: impl Seq<'s>) -> impl FnMut(u8) -> u32;
232
233    /// Hash the given sequence. Ignores `k`.
234    ///
235    /// This is slightly inefficient because it recomputes the constants based on the sequence length.
236    #[inline(always)]
237    fn hash_seq<'s>(&self, seq: impl Seq<'s>) -> u32 {
238        seq.iter_bp().map(self.mapper(seq)).last().unwrap_or(0)
239    }
240    /// Hash all non-empty prefixes of the given sequence. Ignores `k`.
241    #[inline(always)]
242    fn hash_prefixes<'s>(&self, seq: impl Seq<'s>) -> impl ExactSizeIterator<Item = u32> {
243        seq.iter_bp().map(self.mapper(seq))
244    }
245}
seq_hash/lib.rs

seq_hash/
lib.rs
