Skip to main content

zoe/composition/
content.rs

1use crate::{data::types::nucleotides::NucleotidesReadable, private::Sealed, simd::SimdByteFunctions};
2use std::simd::prelude::*;
3
4/// Provides method for obtaining the GC content
5pub trait GcContent: NucleotidesReadable + Sealed {
6    /// Provides the count of G and C bases.
7    #[inline]
8    #[must_use]
9    fn gc_content(&self) -> usize {
10        gc_content_simd::<{ crate::DEFAULT_SIMD_LANES }>(self.nucleotide_bytes())
11    }
12}
13
14impl<T: NucleotidesReadable + Sealed> GcContent for T {}
15
16/// Calculates GC content using SIMD operations for better performance. Returns
17/// count of G and C nucleotides in the sequence.
18#[must_use]
19#[cfg_attr(feature = "multiversion", multiversion::multiversion(targets = "simd"))]
20pub fn gc_content_simd<const N: usize>(s: &[u8]) -> usize {
21    let g = Simd::splat(b'G');
22    let c = Simd::splat(b'C');
23
24    let (pre, mid, sfx) = s.as_simd::<N>();
25
26    let count: usize = mid
27        .iter()
28        .map(|&w| {
29            let v = w.to_ascii_uppercase();
30            (v.simd_eq(g) | v.simd_eq(c)).to_bitmask().count_ones() as usize
31        })
32        .sum();
33
34    count + gc_content(pre) + gc_content(sfx)
35}
36
37/// Calculates GC content using scalar operations. Returns count of G and C
38/// nucleotides in the sequence.
39#[must_use]
40pub fn gc_content(s: &[u8]) -> usize {
41    s.iter()
42        .map(|&b| b.to_ascii_uppercase())
43        .filter(|&b| b == b'G' || b == b'C')
44        .count()
45}
46
47#[cfg(all(test, feature = "rand"))]
48mod test {
49    use super::*;
50
51    #[test]
52    fn test_gc_content() {
53        for length in [16, 1200, 10000] {
54            let s = crate::generate::rand_sequence(b"ATCG", length, 42);
55            let (scalar, simd) = (gc_content(&s), gc_content_simd::<16>(&s));
56            assert_eq!(scalar, simd, "when testing for length {length}");
57        }
58    }
59}
60
61#[cfg(all(test, feature = "rand"))]
62mod bench {
63    use super::*;
64    use std::sync::LazyLock;
65    use test::Bencher;
66
67    extern crate test;
68
69    const N: usize = 1600;
70    static SEQ: LazyLock<Vec<u8>> = std::sync::LazyLock::new(|| crate::generate::rand_sequence(b"ATCG", N, 42));
71
72    #[bench]
73    fn gc_content_scalar(b: &mut Bencher) {
74        b.iter(|| gc_content(&SEQ));
75    }
76
77    #[bench]
78    fn gc_content_simd16(b: &mut Bencher) {
79        b.iter(|| gc_content_simd::<16>(&SEQ));
80    }
81}