bio_seq/codec/
iupac.rs

1//! 4-bit IUPAC nucleotide ambiguity codes
2//!
3//! IUPAC nucleotide ambiguity codes are represented with 4 bits
4//!
5//! |   | A | C | G | T |
6//! | - | - | - | - | - |
7//! | A | 1 | 0 | 0 | 0 |
8//! | C | 0 | 1 | 0 | 0 |
9//! | G | 0 | 0 | 1 | 0 |
10//! | T | 0 | 0 | 0 | 1 |
11//! | Y | 0 | 1 | 0 | 1 |
12//! | R | 1 | 0 | 1 | 0 |
13//! | W | 1 | 0 | 0 | 1 |
14//! | S | 0 | 1 | 1 | 0 |
15//! | K | 0 | 0 | 1 | 1 |
16//! | M | 1 | 1 | 0 | 0 |
17//! | D | 1 | 0 | 1 | 1 |
18//! | V | 1 | 1 | 1 | 0 |
19//! | H | 1 | 1 | 0 | 1 |
20//! | B | 0 | 1 | 1 | 1 |
21//! | N | 1 | 1 | 1 | 1 |
22//! | X/- | 0 | 0 | 0 | 0 |
23//!
24//! The gap symbol [`Iupac::X`] (`-`) represents the empty set. In containment
25//! tests, every symbol contains a gap, while a gap contains only another gap.
26//!
27//! This means that we can treat each symbol as a set and we get meaningful bitwise operations:
28//!
29//! ```rust
30//! use bio_seq::prelude::*;
31//!
32//! // Set union:
33//! let union = iupac!("AS-GYTNAN") | iupac!("ANTGCAT-N");
34//! assert_eq!(union, iupac!("ANTGYWNAN"));
35//!
36//! // Set intersection:
37//! let intersection = iupac!("ACGTSWKMN") & iupac!("WKMSTNNAN");
38//! assert_eq!(intersection, iupac!("A----WKAN"));
39//! ```
40//!
41//! Which can be used to implement pattern matching:
42//!
43//! ```rust
44//! use bio_seq::prelude::*;
45//!
46//! let seq = iupac!("AGCTNNCAGTCGACGTATGTA");
47//! let pattern = iupac!("AYG");
48//!
49//! for slice in seq.windows(pattern.len()) {
50//!    if pattern.contains(slice) {
51//!        println!("{slice} matches pattern");
52//!    }
53//! }
54//!
55//! // ACG matches pattern
56//! // ATG matches pattern
57//! ```
58use crate::codec::{Codec, dna::Dna};
59use crate::seq::{Seq, SeqArray, SeqSlice};
60use crate::{Complement, ComplementMut};
61
62const IUPAC_COMPLEMENT_TABLE: [u8; 16] = {
63    let mut table = [0; 16];
64
65    table[Iupac::A as usize] = Iupac::T as u8;
66    table[Iupac::C as usize] = Iupac::G as u8;
67    table[Iupac::G as usize] = Iupac::C as u8;
68    table[Iupac::T as usize] = Iupac::A as u8;
69    table[Iupac::Y as usize] = Iupac::R as u8;
70    table[Iupac::R as usize] = Iupac::Y as u8;
71    table[Iupac::W as usize] = Iupac::W as u8;
72    table[Iupac::S as usize] = Iupac::S as u8;
73    table[Iupac::K as usize] = Iupac::M as u8;
74    table[Iupac::M as usize] = Iupac::K as u8;
75    table[Iupac::D as usize] = Iupac::H as u8;
76    table[Iupac::V as usize] = Iupac::B as u8;
77    table[Iupac::H as usize] = Iupac::D as u8;
78    table[Iupac::B as usize] = Iupac::V as u8;
79    table[Iupac::N as usize] = Iupac::N as u8;
80
81    table
82};
83
84#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, Codec)]
85#[bits(4)]
86#[repr(u8)]
87pub enum Iupac {
88    A = 0b1000,
89    C = 0b0100,
90    G = 0b0010,
91    T = 0b0001,
92    R = 0b1010,
93    Y = 0b0101,
94    S = 0b0110,
95    W = 0b1001,
96    K = 0b0011,
97    M = 0b1100,
98    B = 0b0111,
99    D = 0b1011,
100    H = 0b1101,
101    V = 0b1110,
102    N = 0b1111,
103    #[display('-')]
104    X = 0b0000,
105}
106
107impl From<Dna> for Iupac {
108    fn from(dna: Dna) -> Self {
109        match dna {
110            Dna::A => Iupac::A,
111            Dna::C => Iupac::C,
112            Dna::G => Iupac::G,
113            Dna::T => Iupac::T,
114        }
115    }
116}
117
118impl Seq<Iupac> {
119    #[must_use]
120    pub fn contains(&self, rhs: &SeqSlice<Iupac>) -> bool {
121        if rhs.len() != self.len() {
122            return false;
123        }
124
125        self.as_ref() & rhs == rhs
126    }
127}
128
129impl<const N: usize, const W: usize> SeqArray<Iupac, N, W> {
130    #[must_use]
131    pub fn contains(&self, rhs: &SeqSlice<Iupac>) -> bool {
132        if N != rhs.len() {
133            return false;
134        }
135        self.as_ref() & rhs == rhs
136    }
137}
138
139impl SeqSlice<Iupac> {
140    #[must_use]
141    pub fn contains(&self, rhs: &SeqSlice<Iupac>) -> bool {
142        if self.len() != rhs.len() {
143            return false;
144        }
145        self & rhs == rhs
146    }
147}
148
149/// The complement of an IUPAC base is the reverse of the bit-pattern
150impl ComplementMut for Iupac {
151    fn comp(&mut self) {
152        // Below are two methods:
153        // 1. Using a lookup table.
154        // 2. The 7 operation from https://graphics.stanford.edu/~seander/bithacks.html
155        // See: https://stackoverflow.com/questions/3587826/is-there-a-built-in-function-to-reverse-bit-order
156
157        // Use a lookup table
158        *self = Iupac::unsafe_from_bits(IUPAC_COMPLEMENT_TABLE[*self as usize]);
159
160        // Use the The 7 operation from
161        // let b = *self as u32;
162        // *self = Iupac::unsafe_from_bits(
163        //     ((((((b * 0x0802u32) & 0x22110u32) | ((b * 0x8020u32) & 0x88440u32)) * 0x10101u32)
164        //         >> 20) // 16 + 4
165        //         & 0x0f) as u8,
166        // );
167    }
168}
169
170impl Complement for Iupac {}
171
172/*
173impl Complement for Seq<Iupac> {
174    type Output = Self;
175
176    fn comp(&mut self) {
177        todo!()
178    }
179
180    fn to_comp(&self) -> Self::Output {
181        todo!()
182    }
183}
184
185/// Reverse complementing a sequence of IUPAC characters is simply a
186/// matter of reversing the entire bit sequence
187impl ReverseComplement for Seq<Iupac> {
188    type Output = Self;
189
190    fn revcomp(&mut self) {
191        // simply do bit-wise reversal
192        self.bv.reverse();
193    }
194
195    fn to_revcomp(&self) -> Self::Output {
196        todo!()
197    }
198}
199*/
200
201#[cfg(test)]
202mod tests {
203    use crate::prelude::*;
204
205    #[test]
206    fn iupac_ops() {
207        let seq = iupac!("AGCTNNCAGTCGACGTATGTASWAGG");
208
209        let pattern = iupac!("NAYGN");
210        let missing = iupac!("NATNG");
211
212        let matches: Vec<Seq<Iupac>> = seq
213            .windows(pattern.len())
214            .filter(|w| pattern.contains(w) && !missing.contains(w))
215            .collect();
216
217        assert_eq!(matches, vec![iupac!("GACGT"), iupac!("TATGT")]);
218    }
219
220    #[test]
221    fn iupac_gap_containment() {
222        let gap_array = SeqArray::<Iupac, 1, 1> {
223            _p: core::marker::PhantomData,
224            ba: crate::Ba::default(),
225        };
226        let gap: Seq<Iupac> = gap_array.as_ref().into();
227        assert_eq!(gap, iupac!("-"));
228
229        for symbol in Iupac::items() {
230            let seq: Seq<Iupac> = std::iter::once(symbol).collect();
231            assert!(seq.contains(&gap));
232            assert!(seq.as_ref().contains(&gap));
233            assert_eq!(gap.contains(&seq), symbol == Iupac::X);
234            assert_eq!(gap.as_ref().contains(&seq), symbol == Iupac::X);
235            assert_eq!(gap_array.contains(&seq), symbol == Iupac::X);
236        }
237    }
238
239    #[test]
240    fn iupac_complement() {
241        assert_eq!(
242            iupac!("AGCTYRWSKMDVHBN").to_comp(),
243            iupac!("TCGARYWSMKHBDVN")
244        );
245    }
246
247    #[test]
248    fn iupac_bits_sanity() {
249        assert_eq!(Iupac::A.to_bits(), 0b1000);
250        assert_eq!(Iupac::C.to_bits(), 0b0100);
251
252        assert_eq!(Iupac::unsafe_from_bits(0b1000), Iupac::A);
253        assert_eq!(Iupac::unsafe_from_bits(0b0100), Iupac::C);
254
255        assert_eq!(Iupac::unsafe_from_bits(Iupac::A.to_bits()), Iupac::A);
256        assert_eq!(Iupac::unsafe_from_bits(Iupac::C.to_bits()), Iupac::C);
257    }
258}
bio_seq/codec/iupac.rs

bio_seq/codec/
iupac.rs
