bio_seq/
kmer.rs

1// Copyright 2021-2024 Jeff Knaggs
2// Licensed under the MIT license (http://opensource.org/licenses/MIT)
3// This file may not be copied, modified, or distributed
4// except according to those terms.
5
6//! Encoded sequences of static length
7//!
8//! Generally, the underlying storage type of `Kmer` should lend itself to optimisation. The default `Kmer` instance is packed into a `usize`, which can be efficiently `Copy`ed on the stack.
9//!
10//! `k * A::BITS` must fit in the storage type, e.g. `usize` (64 bits).
11//!
12//! ```
13//! use bio_seq::prelude::*;
14//!
15//! for (amino_kmer, amino_string) in Seq::<Amino>::try_from("SSLMNHKKL").unwrap()
16//!         .kmers::<3>()
17//!         .zip(["SSL", "SLM", "LMN", "MNH", "NHK", "HKK", "KKL"])
18//!     {
19//!         assert_eq!(amino_kmer, amino_string);
20//!     }
21//! ```
22//!
23//! Kmers can be copied from other sequence types:
24//!
25//! ```
26//! # use bio_seq::prelude::*;
27//! let kmer: Kmer<Dna, 8> = dna!("AGTTGGCA").try_into().unwrap();
28//! ```
29
30use crate::Bs;
31use crate::codec::{self, Codec};
32use crate::prelude::ParseBioError;
33use crate::seq::{Seq, SeqArray, SeqSlice};
34use crate::{
35    Complement, ComplementMut, Reverse, ReverseComplement, ReverseComplementMut, ReverseMut,
36};
37use bitvec::field::BitField;
38use bitvec::view::BitView;
39use core::fmt;
40use core::hash::{Hash, Hasher};
41use core::marker::PhantomData;
42use core::ops::Deref;
43use core::ptr;
44use core::str::FromStr;
45
46pub(crate) mod integral;
47
48#[cfg(feature = "serde")]
49use serde_derive::{Deserialize, Serialize};
50
51const fn make_2bit_table() -> [u8; 256] {
52    let mut table = [0u8; 256];
53    let mut i: usize = 0;
54    while i < 256 {
55        #[expect(
56            clippy::cast_possible_truncation,
57            reason = "the loop bounds i below 256"
58        )]
59        let byte = i as u8;
60        let b0: u8 = (byte & 0b11_00_00_00) >> 6;
61        let b1: u8 = (byte & 0b00_11_00_00) >> 2;
62        let b2: u8 = (byte & 0b00_00_11_00) << 2;
63        let b3: u8 = (byte & 0b00_00_00_11) << 6;
64
65        table[i] = b3 | b2 | b1 | b0;
66        i += 1;
67    }
68    table
69}
70
71const REV_2BIT: [u8; 256] = make_2bit_table();
72
73pub(crate) mod sealed {
74    use crate::Bs;
75    use crate::codec::Codec;
76
77    pub trait KmerStorage: Copy + Clone + PartialEq + std::fmt::Debug {
78        const BITS: usize;
79        type BaN: AsRef<Bs> + AsMut<Bs>;
80
81        fn to_bitarray(self) -> Self::BaN;
82        fn from_bitslice(bs: &Bs) -> Self;
83
84        fn shiftr(&mut self, n: u32);
85
86        fn complement(&mut self, mask: usize);
87        fn rev_blocks<A: Codec, const K: usize>(&mut self);
88    }
89}
90
91pub trait KmerStorage: sealed::KmerStorage {}
92
93impl KmerStorage for usize {}
94
95impl KmerStorage for u64 {}
96
97impl KmerStorage for u128 {}
98
99/// By default k-mers are backed by `usize` and `Codec::BITS` * `K` must be <= 64 on 64-bit platforms
100///
101/// ```compile_fail
102/// # use bio_seq::prelude::*;
103/// // 40 * 2 bits does not fit in a usize
104/// let kmer: Kmer<Dna, 40> = "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT".parse().unwrap();
105/// ```
106///
107/// ```
108/// # use bio_seq::prelude::*;
109/// // ..but 80 bits does fit into a u128
110/// let kmer: Kmer<Dna, 40, u128> = "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT".parse().unwrap();
111/// assert_eq!(kmer.len(), 40);
112/// ```
113#[derive(Debug, PartialEq, Eq, PartialOrd, Ord, Copy, Clone)]
114#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
115#[repr(transparent)]
116pub struct Kmer<C: Codec, const K: usize, S: KmerStorage = usize> {
117    pub(crate) _p: PhantomData<C>,
118    pub(crate) bs: S,
119}
120
121impl<A: Codec, const K: usize, S: KmerStorage> Kmer<A, K, S> {
122    // This error message can be formatted with constants in nightly (const_format)
123    const ASSERT_K: () = assert!(
124        K * A::BITS as usize <= S::BITS,
125        "`KmerStorage` not large enough for `Kmer`",
126    );
127
128    const ASSERT_K_NONZERO: () = assert!(K > 0, "`K` must be greater than 0");
129
130    const fn assert_k() {
131        let () = Self::ASSERT_K;
132        let () = Self::ASSERT_K_NONZERO;
133    }
134
135    const BITS: usize = K * A::BITS as usize;
136
137    #[must_use]
138    pub const fn len(&self) -> usize {
139        K
140    }
141
142    #[must_use]
143    pub const fn is_empty(&self) -> bool {
144        // This is recommended by clippy since we have `len`
145        // Kmers are never empty if K > 0
146        false
147    }
148
149    #[must_use]
150    pub fn rotated_left(&self, n: u32) -> Self {
151        let n: usize = (n as usize % K) * A::BITS as usize;
152        let mut ba = self.bs.to_bitarray();
153        let bs: &mut Bs = ba.as_mut();
154        bs[..Self::BITS].rotate_left(n);
155
156        Kmer {
157            _p: PhantomData,
158            bs: S::from_bitslice(&bs[..Self::BITS]),
159        }
160    }
161
162    #[must_use]
163    pub fn rotated_right(&self, n: u32) -> Self {
164        let n: usize = (n as usize % K) * A::BITS as usize;
165        let mut ba = self.bs.to_bitarray();
166        let bs: &mut Bs = ba.as_mut();
167        bs[..Self::BITS].rotate_right(n);
168
169        Kmer {
170            _p: PhantomData,
171            bs: S::from_bitslice(&bs[..Self::BITS]),
172        }
173    }
174
175    /// Shift bases to the right and push a base onto the end.
176    ///
177    /// ```
178    /// use bio_seq::prelude::*;
179    /// use bio_seq::codec::dna::Dna;
180    ///
181    /// let k = kmer!("ACGAT");
182    /// assert_eq!(k.pushr(Dna::T).to_string(), "CGATT");
183    /// ```
184    #[must_use]
185    pub fn pushr(self, base: A) -> Self {
186        let mut ba = self.rotated_left(1).bs.to_bitarray();
187        let bs: &mut Bs = ba.as_mut();
188
189        let start = Self::BITS - A::BITS as usize;
190        let end = start + A::BITS as usize;
191
192        bs[start..end].store(base.to_bits());
193
194        Kmer {
195            _p: PhantomData,
196            bs: S::from_bitslice(bs),
197        }
198    }
199
200    /// Push a base from the left
201    #[must_use]
202    pub fn pushl(self, base: A) -> Self {
203        let mut ba = self.rotated_right(1).bs.to_bitarray();
204        let bs: &mut Bs = ba.as_mut();
205
206        bs[..A::BITS as usize].store(base.to_bits());
207
208        Kmer {
209            _p: PhantomData,
210            bs: S::from_bitslice(bs),
211        }
212    }
213
214    /// Create Kmer from sequence without checking length
215    #[must_use]
216    pub fn unsafe_from_seqslice(seq: &SeqSlice<A>) -> Self {
217        Self::assert_k();
218        debug_assert!(K == seq.len(), "K != seq.len()");
219        Kmer {
220            _p: PhantomData,
221            bs: S::from_bitslice(&seq.bs),
222        }
223    }
224
225    fn complement(&mut self) {
226        self.bs.complement(K * A::BITS as usize);
227    }
228
229    fn rev_blocks(&mut self) {
230        self.bs.rev_blocks::<A, K>();
231        let shift = u32::try_from(S::BITS - Self::BITS).expect("k-mer storage is at most 128 bits");
232        self.bs.shiftr(shift);
233    }
234}
235
236impl<A: Codec, const K: usize> From<usize> for Kmer<A, K, usize> {
237    fn from(i: usize) -> Kmer<A, K, usize> {
238        Self::assert_k();
239        Kmer {
240            _p: PhantomData,
241            bs: i,
242        }
243    }
244}
245
246impl<A: Codec, const K: usize> From<u64> for Kmer<A, K, u64> {
247    fn from(i: u64) -> Kmer<A, K, u64> {
248        Self::assert_k();
249        Kmer {
250            _p: PhantomData,
251            bs: i,
252        }
253    }
254}
255
256impl<A: Codec, const K: usize> From<usize> for Kmer<A, K, u64> {
257    fn from(i: usize) -> Kmer<A, K, u64> {
258        Self::assert_k();
259        Kmer {
260            _p: PhantomData,
261            bs: i as u64,
262        }
263    }
264}
265
266/*
267impl<A: Codec, const K: usize, S: KmerStorage> From<&SeqSlice<A>> for Kmer<A, K, S> {
268    fn from(seq: &SeqSlice<A>) -> Self {
269        Kmer {
270            _p: PhantomData,
271            bs: S::from_bitslice(&seq.bs),
272        }
273    }
274}
275*/
276
277impl<S: KmerStorage + Into<usize>, A: Codec, const K: usize> From<&Kmer<A, K, S>> for usize {
278    fn from(kmer: &Kmer<A, K, S>) -> usize {
279        kmer.bs.into()
280    }
281}
282
283impl<A: Codec, const K: usize> Deref for Kmer<A, K, usize> {
284    type Target = SeqSlice<A>;
285
286    fn deref(&self) -> &Self::Target {
287        let bs: &Bs = &self.bs.view_bits()[0..(K * A::BITS as usize)];
288        let bs: *const Bs = ptr::from_ref::<Bs>(bs);
289        unsafe { &*(bs as *const SeqSlice<A>) }
290    }
291}
292
293impl<A: Codec, const K: usize> AsRef<SeqSlice<A>> for Kmer<A, K, usize> {
294    fn as_ref(&self) -> &SeqSlice<A> {
295        self
296    }
297}
298
299impl<A: Codec, const K: usize, S: KmerStorage> fmt::Display for Kmer<A, K, S> {
300    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
301        let mut s = String::new();
302        let ba = self.bs.to_bitarray();
303        let bs: &Bs = &ba.as_ref()[0..(K * A::BITS as usize)];
304
305        bs.chunks(A::BITS as usize).for_each(|chunk| {
306            s.push(A::unsafe_from_bits(chunk.load_le::<u8>()).to_char());
307        });
308        write!(f, "{s}")
309    }
310}
311
312/// An iterator over all kmers of a sequence with a specified length
313#[must_use = "iterators are lazy and do nothing unless consumed"]
314pub struct KmerIter<'a, A: Codec, const K: usize> {
315    pub(crate) slice: &'a SeqSlice<A>,
316    pub(crate) index: usize,
317    pub(crate) len: usize,
318    pub(crate) _p: PhantomData<A>,
319}
320
321impl<A: Codec, const K: usize, S: KmerStorage> Kmer<A, K, S> {
322    fn unsafe_from(seq: &SeqSlice<A>) -> Self {
323        Self::assert_k();
324        debug_assert!(K == seq.len(), "K != seq.len()");
325        Kmer {
326            _p: PhantomData,
327            bs: S::from_bitslice(&seq.bs),
328        }
329    }
330}
331
332impl<A: Codec, const K: usize> Iterator for KmerIter<'_, A, K> {
333    type Item = Kmer<A, K>;
334    fn next(&mut self) -> Option<Kmer<A, K>> {
335        let i = self.index;
336        if self.index + K > self.len {
337            return None;
338        }
339        self.index += 1;
340        Some(Kmer::<A, K>::unsafe_from(&self.slice[i..i + K]))
341    }
342
343    fn size_hint(&self) -> (usize, Option<usize>) {
344        let n = (self.len.saturating_sub(self.index) + 1).saturating_sub(K);
345        (n, Some(n))
346    }
347}
348
349impl<A: Codec, const K: usize> ExactSizeIterator for KmerIter<'_, A, K> {}
350
351/// ```
352/// use bio_seq::prelude::*;
353/// use std::hash::{Hash, Hasher, DefaultHasher};
354///
355/// let mut hasher1 = DefaultHasher::new();
356/// kmer!("AAA").hash(&mut hasher1);
357/// let hash1 = hasher1.finish();
358///
359/// let mut hasher2 = DefaultHasher::new();
360/// kmer!("AAAA").hash(&mut hasher2);
361/// let hash2 = hasher2.finish();
362///
363/// assert_ne!(hash1, hash2);
364/// ```
365impl<A: Codec, const K: usize, S: KmerStorage> Hash for Kmer<A, K, S> {
366    fn hash<H: Hasher>(&self, state: &mut H) {
367        // Length prefix as 8 little-endian bytes (see `SeqSlice` for rationale).
368        state.write(&(K as u64).to_le_bytes());
369        let ba = self.bs.to_bitarray();
370        let bs: &Bs = ba.as_ref();
371        // Only the meaningful bits are hashed; the unused high bits of the
372        // storage word/array are padding and must not influence the digest.
373        let n_bits = K * A::BITS as usize;
374        crate::hash::hash_bits(&bs[..n_bits], state);
375    }
376}
377
378impl<A: Codec, const K: usize, S: KmerStorage> TryFrom<&SeqSlice<A>> for Kmer<A, K, S> {
379    type Error = ParseBioError;
380
381    fn try_from(seq: &SeqSlice<A>) -> Result<Self, Self::Error> {
382        if seq.len() == K {
383            Ok(Kmer::<A, K, S>::unsafe_from(&seq[0..K]))
384        } else {
385            Err(ParseBioError::MismatchedLength(K, seq.len()))
386        }
387    }
388}
389
390impl<A: Codec, const K: usize> TryFrom<Seq<A>> for Kmer<A, K> {
391    type Error = ParseBioError;
392
393    fn try_from(seq: Seq<A>) -> Result<Self, Self::Error> {
394        Self::try_from(seq.as_ref())
395    }
396}
397
398impl<A: Codec, const K: usize, S: KmerStorage> PartialEq<SeqArray<A, K, 1>> for Kmer<A, K, S> {
399    fn eq(&self, seq: &SeqArray<A, K, 1>) -> bool {
400        if seq.len() != K {
401            return false;
402        }
403        &Kmer::<A, K, S>::unsafe_from(seq.as_ref()) == self
404    }
405}
406
407impl<A: Codec, const K: usize, S: KmerStorage> PartialEq<&SeqArray<A, K, 1>> for Kmer<A, K, S> {
408    fn eq(&self, seq: &&SeqArray<A, K, 1>) -> bool {
409        if seq.len() != K {
410            return false;
411        }
412        &Kmer::<A, K, S>::unsafe_from(seq.as_ref()) == self
413    }
414}
415
416impl<A: Codec, const K: usize> PartialEq<Seq<A>> for Kmer<A, K> {
417    fn eq(&self, seq: &Seq<A>) -> bool {
418        if seq.len() != K {
419            return false;
420        }
421        &Kmer::<A, K>::unsafe_from(seq.as_ref()) == self
422    }
423}
424
425impl<A: Codec, const K: usize, S: KmerStorage> PartialEq<SeqSlice<A>> for Kmer<A, K, S> {
426    fn eq(&self, seq: &SeqSlice<A>) -> bool {
427        if seq.len() != K {
428            return false;
429        }
430        &Kmer::<A, K, S>::unsafe_from(seq) == self
431    }
432}
433
434impl<A: Codec, const K: usize, S: KmerStorage> PartialEq<&SeqSlice<A>> for Kmer<A, K, S> {
435    fn eq(&self, seq: &&SeqSlice<A>) -> bool {
436        if seq.len() != K {
437            return false;
438        }
439        &Kmer::<A, K, S>::unsafe_from(seq) == self
440    }
441}
442
443impl<A: Codec, const K: usize> PartialEq<&str> for Kmer<A, K> {
444    fn eq(&self, seq: &&str) -> bool {
445        &self.to_string() == seq
446    }
447}
448
449impl<A: Codec, const K: usize, S: KmerStorage> FromStr for Kmer<A, K, S> {
450    type Err = ParseBioError;
451
452    fn from_str(s: &str) -> Result<Self, Self::Err> {
453        if s.len() != K {
454            return Err(ParseBioError::MismatchedLength(K, s.len()));
455        }
456        let seq: Seq<A> = Seq::from_str(s)?;
457        Kmer::<A, K, S>::try_from(seq.as_ref())
458    }
459}
460
461impl<A: Codec, const K: usize> From<Kmer<A, K, usize>> for Seq<A> {
462    fn from(kmer: Kmer<A, K, usize>) -> Self {
463        let mut seq: Seq<A> = Seq::with_capacity(K);
464        seq.extend(kmer.iter());
465        seq
466    }
467}
468
469impl<const K: usize, S: KmerStorage> ComplementMut for Kmer<codec::dna::Dna, K, S> {
470    fn comp(&mut self) {
471        self.complement();
472    }
473}
474
475impl<const K: usize, S: KmerStorage> Complement for Kmer<codec::dna::Dna, K, S> {}
476
477impl<A: Codec, const K: usize, S: KmerStorage> ReverseMut for Kmer<A, K, S> {
478    fn rev(&mut self) {
479        self.rev_blocks();
480    }
481}
482
483impl<A: Codec, const K: usize, S: KmerStorage> Reverse for Kmer<A, K, S> {}
484
485impl<const K: usize> ReverseComplementMut for Kmer<codec::dna::Dna, K, usize> {}
486
487impl<const K: usize> ReverseComplement for Kmer<codec::dna::Dna, K, usize> {}
488
489/// Convenient compile time kmer constructor
490///
491/// This is a wrapper for the `dna!` macro that returns a `Kmer`:
492/// ```
493/// # use bio_seq::prelude::*;
494/// let kmer: Kmer<Dna, 8> = kmer!("ACGTACGT");
495/// ```
496#[macro_export]
497macro_rules! kmer {
498    ($seq:expr) => {
499        $crate::kmer::Kmer::<$crate::codec::dna::Dna, { $seq.len() }>::unsafe_from_seqslice(dna!($seq))
500    };
501    ($seq:expr, $storage:ty) => {
502        $crate::kmer::Kmer::<$crate::codec::dna::Dna, { $seq.len() }, $storage>::unsafe_from_seqslice(
503            dna!($seq),
504        )
505    };
506}
507
508#[cfg(test)]
509mod tests {
510    use crate::prelude::*;
511    use crate::seq::SeqArray;
512
513    #[cfg(target_pointer_width = "64")]
514    #[test]
515    fn kmer_to_usize() {
516        let s: &'static SeqSlice<Dna> = dna!("AACTT");
517        println!("{s}");
518
519        for (kmer, index) in s.kmers::<2>().zip([0b00_00, 0b01_00, 0b11_01, 0b11_11]) {
520            println!("{kmer}");
521            assert_eq!(index, usize::from(&kmer));
522        }
523    }
524    #[test]
525    fn pushl_test() {
526        let k = kmer!("ACGT");
527        let k1 = k.pushl(Dna::G);
528        let k2 = k1.pushl(Dna::A);
529        let k3 = k2.pushl(Dna::T);
530        let k4 = k3.pushl(Dna::C);
531        let k5 = k4.pushl(Dna::C);
532
533        assert_eq!(k1, kmer!("GACG"));
534        assert_eq!(k2, kmer!("AGAC"));
535        assert_eq!(k3, kmer!("TAGA"));
536        assert_eq!(k4, kmer!("CTAG"));
537        assert_eq!(k5, kmer!("CCTA"));
538    }
539    #[test]
540    fn pushr_test() {
541        let k = kmer!("ACGT");
542        let k1 = k.pushr(Dna::G);
543        let k2 = k1.pushr(Dna::A);
544        let k3 = k2.pushr(Dna::T);
545        let k4 = k3.pushr(Dna::C);
546        let k5 = k4.pushr(Dna::C);
547
548        println!("{k1}");
549        assert_eq!(k1, kmer!("CGTG"));
550        assert_eq!(k2, kmer!("GTGA"));
551        assert_eq!(k3, kmer!("TGAT"));
552        assert_eq!(k4, kmer!("GATC"));
553        assert_eq!(k5, kmer!("ATCC"));
554    }
555
556    #[test]
557    fn amino_kmer_to_usize() {
558        for (kmer, index) in Seq::<Amino>::try_from("SRY")
559            .unwrap()
560            .kmers::<2>()
561            .zip([0b0010_0001_1000, 0b0100_1100_1000])
562        {
563            assert_eq!(index, usize::from(&kmer));
564        }
565    }
566    #[test]
567    fn big_kmer_shiftr() {
568        let mut kmer: Kmer<Dna, 32, u64> = kmer!("AATTTGTGGGTTCGTCTGCGGCTCCGCCCTTA", u64);
569        for base in dna!("TACTATGAGGACGATCAGCACCATAAGAACAAA") {
570            kmer = kmer.pushr(base);
571        }
572        assert_eq!(kmer!("ACTATGAGGACGATCAGCACCATAAGAACAAA", u64), kmer);
573    }
574
575    #[test]
576    fn big_kmer_shiftl() {
577        let mut kmer: Kmer<Dna, 32, u64> = kmer!("AATTTGTGGGTTCGTCTGCGGCTCCGCCCTTA", u64);
578        for base in dna!("GTACTATGAGGACGATCAGCACCATAAGAACAAA") {
579            kmer = kmer.pushl(base);
580        }
581        assert_eq!(kmer!("AAACAAGAATACCACGACTAGCAGGAGTATCA", u64), kmer);
582    }
583
584    #[test]
585    fn amino_kmer_iter() {
586        for (kmer, target) in Seq::<Amino>::try_from("SSLMNHKKL")
587            .unwrap()
588            .kmers::<3>()
589            .zip(["SSL", "SLM", "LMN", "MNH", "NHK", "HKK", "KKL"])
590        {
591            assert_eq!(kmer, target);
592        }
593    }
594
595    #[test]
596    fn test_rotations() {
597        let kmer: Kmer<Dna, 9> = Kmer::try_from(dna!("ACTGCGATG")).unwrap();
598
599        for (rotation, shift) in [
600            "ACTGCGATG",
601            "CTGCGATGA",
602            "TGCGATGAC",
603            "GCGATGACT",
604            "CGATGACTG",
605            "GATGACTGC",
606            "ATGACTGCG",
607            "TGACTGCGA",
608            "GACTGCGAT",
609            "ACTGCGATG",
610            "CTGCGATGA",
611            "TGCGATGAC",
612        ]
613        .into_iter()
614        .zip(0u32..)
615        {
616            assert_eq!(kmer.rotated_left(shift), rotation);
617        }
618
619        for (rotation, shift) in [
620            "ACTGCGATG",
621            "GACTGCGAT",
622            "TGACTGCGA",
623            "ATGACTGCG",
624            "GATGACTGC",
625            "CGATGACTG",
626            "GCGATGACT",
627            "TGCGATGAC",
628            "CTGCGATGA",
629            "ACTGCGATG",
630            "GACTGCGAT",
631        ]
632        .into_iter()
633        .zip(0u32..)
634        {
635            assert_eq!(kmer.rotated_right(shift), rotation);
636        }
637
638        let kmer: Kmer<Dna, 8> = Kmer::try_from(dna!("ACTGCGAT")).unwrap().rotated_left(1);
639
640        assert_ne!(kmer.to_string(), "ACTGCGAT");
641        assert_eq!(kmer.to_string(), "CTGCGATA");
642
643        let kmer: Kmer<Dna, 8> = kmer!("ACTGCGAT").rotated_right(1);
644
645        assert_ne!(kmer.to_string(), "ACTGCGAT");
646        assert_eq!(kmer.to_string(), "TACTGCGA");
647
648        let kmer: Kmer<Dna, 9> = Kmer::from_str("ACTGCGATG").unwrap().rotated_left(0);
649
650        assert_eq!(kmer.to_string(), "ACTGCGATG");
651        assert_ne!(kmer.to_string(), "ACTGCGATGA");
652
653        let kmer: Kmer<Dna, 9> = Kmer::try_from(dna!("ACTGCGATG")).unwrap().rotated_right(0);
654
655        assert_eq!(kmer.to_string(), "ACTGCGATG");
656        assert_ne!(kmer.to_string(), "ACTGCGATGA");
657
658        let kmer: Kmer<Dna, 9> = Kmer::from_str("ACTGCGATG").unwrap().rotated_left(9 * 3307);
659
660        assert_eq!(kmer.to_string(), "ACTGCGATG");
661        assert_ne!(kmer.to_string(), "ACTGCGATGA");
662
663        let kmer: Kmer<Dna, 9> = kmer!("ACTGCGATG").rotated_right(9 * 3307);
664
665        assert_eq!(kmer.to_string(), "ACTGCGATG");
666        assert_ne!(kmer.to_string(), "ACTGCGATGA");
667    }
668
669    #[test]
670    fn eq_functions() {
671        assert_eq!(kmer!("ACGT"), dna!("ACGT"));
672
673        // this should be a compiler error:
674        // assert_ne!(kmer!("ACGT"), dna!("ACGTA"));
675
676        let kmer: Kmer<Iupac, 4> = Kmer::from_str("ACGT").unwrap();
677        assert_eq!(kmer, iupac!("ACGT"));
678        assert_ne!(kmer, iupac!("NCGT"));
679    }
680
681    #[test]
682    fn kmer_iter() {
683        //let seq = dna!("ACTGA");
684        let cs: Vec<Kmer<Dna, 3>> = dna!("ACTGA").kmers().collect();
685        assert_eq!(cs[0], "ACT");
686        assert_eq!(cs[1], "CTG");
687        assert_eq!(cs[2], "TGA");
688        assert_eq!(cs.len(), 3);
689    }
690
691    #[test]
692    fn k_check() {
693        let _kmer = Kmer::<Dna, 32>::from(0);
694        let _kmer = Kmer::<Amino, 10>::from(0);
695        let _kmer = Kmer::<Iupac, 14>::from(0);
696    }
697
698    #[test]
699    fn kmer_revcomp() {
700        assert_eq!(kmer!("ACGT"), kmer!("ACGT").to_revcomp());
701        assert_ne!(kmer!("GTCGTA"), kmer!("TACGAC"));
702
703        let rc = kmer!("GTCGTA").to_revcomp();
704
705        assert_eq!(rc, kmer!("TACGAC"));
706
707        assert_eq!(
708            kmer!("GCTATCGATCTGATCG"),
709            kmer!("CGATCAGATCGATAGC").to_revcomp()
710        );
711    }
712
713    #[test]
714    fn kmer_deref() {
715        let kmer: Kmer<Dna, 3> = kmer!("ACG");
716        let seq: &SeqSlice<Dna> = &kmer;
717
718        assert_eq!(*seq, *kmer);
719        assert_eq!(seq.to_string(), "ACG");
720
721        let kmer: Kmer<Dna, 8> = kmer!("AAAAAAAA");
722        let seq: &SeqSlice<Dna> = &kmer;
723
724        assert_eq!(seq.to_string(), "AAAAAAAA");
725
726        let kmer: Kmer<Dna, 8> = kmer!("TTTTTTTT");
727        let seq: &SeqSlice<Dna> = &kmer;
728
729        assert_eq!(seq.to_string(), "TTTTTTTT");
730
731        let kmer: Kmer<Dna, 16> = kmer!("AGCTAGCTAGCTAGCT");
732        let seq: &SeqSlice<Dna> = &kmer;
733
734        assert_eq!(seq.to_string(), "AGCTAGCTAGCTAGCT");
735    }
736
737    #[test]
738    fn kmer_as_ref() {
739        let kmer: Kmer<Dna, 4> = kmer!("ACGT");
740        let seq: &SeqSlice<Dna> = kmer.as_ref();
741
742        assert_eq!(seq.to_string(), "ACGT");
743
744        let kmer: Kmer<Dna, 16> = kmer!("AGCTAGCTAGCTAGCT");
745        let seq: &SeqSlice<Dna> = kmer.as_ref();
746
747        assert_eq!(seq.to_string(), "AGCTAGCTAGCTAGCT");
748    }
749
750    #[test]
751    fn kmer_rev() {
752        let mut kmer: Kmer<Dna, 4> = kmer!("ACGT");
753
754        kmer.rev();
755
756        assert_eq!(kmer.to_string(), "TGCA");
757
758        kmer.rev();
759
760        assert_eq!(kmer.to_string(), "ACGT");
761
762        kmer.rev();
763
764        assert_eq!(kmer.to_string(), "TGCA");
765    }
766
767    #[test]
768    fn kmer_storage_types() {
769        let s1 = "AACGTAGCCGCGAACTTACGTAGCCGCGAAAA";
770        let s2 = "AACGTAGCCGCGAACTTACGTAGCCGCGAAA";
771        let s3 = "ACGTAGCCGCGAACTTACGTAGCCGCGAAAA";
772
773        let s4 = "AACGTAGCCGCGAACTTACGTAGCCGCGAAAAAACGTAGCCGCGAACTTACGTAGCCGCGAAAA";
774        let s5 = "AACGTAGCCGCGAACTTACGTAGCCGCGAAAAAACGTAGCCGCGAACTTACGTAGCCGCGAAAAA";
775
776        assert_eq!(s1.len(), 32);
777        assert_eq!(s2.len(), 31);
778        assert_eq!(s3.len(), 31);
779        assert_eq!(s4.len(), 64);
780        assert_eq!(s5.len(), 65);
781
782        let kmer1_64 = Kmer::<Dna, 32, u64>::from_str(s1).unwrap();
783        let kmer2_64 = Kmer::<Dna, 31, u64>::from_str(s2).unwrap();
784        let kmer3_64 = Kmer::<Dna, 31, u64>::from_str(s3).unwrap();
785
786        let kmer1 = Kmer::<Dna, 32, u64>::from_str(s1).unwrap();
787        let kmer2 = Kmer::<Dna, 31, u64>::from_str(s2).unwrap();
788        let kmer3 = Kmer::<Dna, 31, u64>::from_str(s3).unwrap();
789
790        let kmer4_128 = Kmer::<Dna, 64, u128>::from_str(s4).unwrap();
791
792        let seq5: Seq<Dna> = s5.try_into().unwrap();
793
794        assert_eq!(kmer4_128, &seq5[..64]);
795        assert_ne!(kmer4_128, &seq5[1..]);
796
797        assert_eq!(kmer1, &seq5[..32]);
798        assert_eq!(kmer1, &seq5[32..64]);
799
800        assert_eq!(kmer1_64, &seq5[..32]);
801        assert_eq!(kmer1_64, &seq5[32..64]);
802
803        assert_ne!(kmer1, &seq5[..31]);
804        assert_ne!(kmer1, &seq5[32..]);
805
806        assert_ne!(kmer1_64, &seq5[1..33]);
807        assert_ne!(kmer1_64, &seq5[33..]);
808
809        assert_eq!(kmer4_128, kmer4_128);
810        assert_eq!(kmer1_64, kmer1_64);
811        assert_eq!(kmer2, kmer2);
812
813        assert_ne!(kmer2, kmer3);
814        assert_ne!(kmer2_64, kmer3_64);
815        // PartialEq is not implemented for different storgage types
816        /*
817                assert_ne!(kmer2, kmer3_64);
818                assert_eq!(kmer2, kmer2_64);
819                assert_eq!(kmer1, kmer1_64);
820        */
821    }
822
823    #[test]
824    fn try_from_seq() {
825        let seq: Seq<Dna> = Seq::try_from("ACACACACACACGT").unwrap();
826        assert_eq!(
827            Kmer::<Dna, 8>::try_from(&seq[..8]).unwrap().to_string(),
828            "ACACACAC"
829        );
830        assert_eq!(
831            Kmer::<Dna, 8>::try_from(&seq[1..9]).unwrap().to_string(),
832            "CACACACA"
833        );
834
835        let err: Result<Kmer<Dna, 8>, ParseBioError> = Kmer::try_from(&seq[2..9]);
836        assert_eq!(err, Err(ParseBioError::MismatchedLength(8, 7)));
837
838        let err: Result<Kmer<Dna, 8>, ParseBioError> = Kmer::try_from(seq);
839        assert_eq!(err, Err(ParseBioError::MismatchedLength(8, 14)));
840
841        let seq: Seq<Dna> = Seq::try_from("ACACACACACACGT").unwrap();
842
843        assert_eq!(
844            Kmer::<Dna, 14>::try_from(seq).unwrap().to_string(),
845            "ACACACACACACGT"
846        );
847    }
848
849    #[test]
850    fn two_bit_reversal_table() {
851        let generate = std::hint::black_box(super::make_2bit_table as fn() -> [u8; 256]);
852        let table = generate();
853
854        for input in 0u8..=u8::MAX {
855            let mut rem = input;
856            let mut expected = 0u8;
857
858            for _ in 0..4 {
859                expected = (expected << 2) | (rem & 0b11);
860                rem >>= 2;
861            }
862
863            let i = usize::from(input);
864            assert_eq!(table[i], expected, "input={input:#010b}");
865            assert_eq!(super::REV_2BIT[i], expected, "input={input:#010b}");
866        }
867    }
868
869    #[test]
870    fn kmer_is_nonempty() {
871        let k: Kmer<Dna, 1> = "A".parse().unwrap();
872        assert!(!k.is_empty());
873    }
874}
bio_seq/kmer.rs

bio_seq/
kmer.rs
