symcurve/curve/
matrix.rs

1//! This module contains some constants/matrices for curvature calculation.
2use std::fmt;
3
4/// The number of nucleotides in a triplet, which is also the number of dimensions in the
5/// nucleotide matrices used for triplet -> value lookup.
6pub const TRIPLET_SIZE: usize = 3;
7
8/// A type alias for a 3D matrix sized 4x4x4 of f64 values. The first dimension is the
9/// first nucleotide in a triplet, the second dimension is the second nucleotide in a triplet,
10/// and the third dimension is the third nucleotide in a triplet.
11pub type NucMatrix = [[[f64; 4]; 4]; 4];
12
13/// The TWIST matrix is used to calculate the twist angle in three nucleotides of DNA.
14/// The values are all 0.598647428 for all combinations of nucleotide triplets.
15pub const TWIST: NucMatrix = [[[0.598647428; 4]; 4]; 4];
16
17/// The TILT matrix is not really used in the current implementation, but is included here
18/// for completeness. The values are all 0.0 for all combinations of nucleotide triplets.
19pub const TILT: NucMatrix = [[[0.0; 4]; 4]; 4];
20
21/// The "activated" version of the ROLL matrix is used to calculate the roll angle in three
22/// nucleotides of DNA. This matrix differs from the simple matrix in that the angle
23/// values represent a more activated state of the nucleosomes bound to the DNA.
24///
25/// Selected by [`RollType::Active`].
26pub const ROLL_ACTIVE: NucMatrix = [
27    [
28        [0.0633, 0.3500, 4.6709, 2.64115],
29        [6.2734, 0.3500, 7.7171, 4.44325],
30        [4.8884, 3.9232, 5.0523, 6.8829],
31        [5.4903, 3.9232, 5.3055, 5.3055],
32    ],
33    [
34        [4.6709, 6.2734, 5.00295, 5.0673],
35        [4.6709, 0.0633, 4.7618, 4.0633],
36        [7.7000, 5.4903, 3.05865, 6.75525],
37        [7.7000, 4.8884, 7.07195, 4.9907],
38    ],
39    [
40        [4.0633, 4.44325, 5.9806, 5.51645],
41        [5.0673, 2.64115, 6.62555, 5.51645],
42        [4.9907, 5.3055, 5.89135, 9.0823],
43        [6.75525, 6.8829, 5.89135, 9.0823],
44    ],
45    [
46        [4.7618, 7.7171, 6.8996, 6.62555],
47        [5.00295, 4.6709, 6.8996, 5.9806],
48        [7.07195, 5.3055, 3.869, 5.9000],
49        [3.05865, 5.0523, 3.869, 5.827],
50    ],
51];
52
53/// The simple version of the ROLL matrix is used to calculate the roll angle in three nucleotides
54/// of DNA.
55///
56/// Selected by [`RollType::Simple`].
57pub const ROLL_SIMPLE: NucMatrix = [
58    [
59        [0.1, 0.0, 4.2, 1.6],
60        [9.7, 0.0, 8.7, 3.6],
61        [6.5, 2.0, 4.7, 6.3],
62        [5.8, 2.0, 5.2, 5.2],
63    ],
64    [
65        [7.3, 9.7, 7.8, 6.4],
66        [7.3, 0.1, 6.2, 5.1],
67        [10.0, 5.8, 0.7, 7.5],
68        [10.0, 6.5, 5.8, 6.2],
69    ],
70    [
71        [5.1, 3.6, 6.6, 5.6],
72        [6.4, 1.6, 6.8, 5.6],
73        [6.2, 5.2, 5.7, 8.2],
74        [7.5, 6.3, 4.3, 8.2],
75    ],
76    [
77        [6.2, 8.7, 9.6, 6.8],
78        [7.8, 4.2, 9.6, 6.6],
79        [5.8, 5.2, 3.0, 4.3],
80        [0.7, 4.7, 3.0, 5.7],
81    ],
82];
83
84/// Why a matrix lookup could not be performed.
85///
86/// An enum rather than a message string: the two cases are distinct, callers can tell
87/// them apart, and neither needs an allocation to report.
88#[derive(Debug, Clone, Copy, PartialEq, Eq)]
89pub enum MatrixLookupError {
90    /// The slice handed in was not exactly three bases long.
91    WrongLength(usize),
92    /// A base was not A, C, G or T in either case.
93    UnknownBase(u8),
94}
95
96impl fmt::Display for MatrixLookupError {
97    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
98        match self {
99            Self::WrongLength(n) => write!(f, "triplet must be of length 3, got {n}"),
100            Self::UnknownBase(base) => {
101                write!(f, "unrecognized nucleotide {:?}", *base as char)
102            }
103        }
104    }
105}
106
107impl std::error::Error for MatrixLookupError {}
108
109/// Which ROLL matrix a curvature calculation should use.
110#[derive(Debug, Clone, Copy)]
111pub enum RollType {
112    Simple,
113    Active,
114}
115
116/// Sentinel stored in `NUC_TABLE` for any byte that is not a scoreable base.
117const NOT_A_BASE: u8 = u8::MAX;
118
119/// Direct lookup from an ASCII byte to its index in a `NucMatrix`.
120///
121/// Both cases are populated, so soft-masked sequence (RepeatMasker lowercases
122/// repetitive regions) is treated as ordinary sequence rather than as unknown bases.
123/// Everything else, including `N` and the IUPAC ambiguity codes, maps to `NOT_A_BASE`;
124/// callers are expected to have split those out already.
125///
126/// This runs three times per base over an entire genome, so it is a flat table rather
127/// than a case conversion followed by a match.
128const NUC_TABLE: [u8; 256] = {
129    let mut table = [NOT_A_BASE; 256];
130    table[b'A' as usize] = 0;
131    table[b'a' as usize] = 0;
132    table[b'T' as usize] = 1;
133    table[b't' as usize] = 1;
134    table[b'G' as usize] = 2;
135    table[b'g' as usize] = 2;
136    table[b'C' as usize] = 3;
137    table[b'c' as usize] = 3;
138    table
139};
140
141/// Maps a nucleotide to its index in a `NucMatrix`, in either case.
142///
143/// Returns `None` for anything that is not A, C, G, or T.
144#[inline]
145pub fn nuc_index(base: u8) -> Option<usize> {
146    match NUC_TABLE[base as usize] {
147        NOT_A_BASE => None,
148        ix => Some(ix as usize),
149    }
150}
151
152/// Decodes a triplet into its three `NucMatrix` indices.
153///
154/// Several matrices are consulted for the same triplet on every base, so the ASCII
155/// decode is done once here and the resulting indices reused, rather than re-decoding
156/// per matrix. The result is a fixed-size array, so no allocation is involved.
157///
158/// # Errors
159///
160/// Returns a `MatrixLookupError` naming the offending base if any of the three is not
161/// A, C, G, or T.
162pub fn triplet_indices(
163    triplet: &[u8; TRIPLET_SIZE],
164) -> Result<[usize; TRIPLET_SIZE], MatrixLookupError> {
165    let mut ixs = [0usize; TRIPLET_SIZE];
166    for (slot, &base) in ixs.iter_mut().zip(triplet.iter()) {
167        *slot = nuc_index(base).ok_or(MatrixLookupError::UnknownBase(base))?;
168    }
169    Ok(ixs)
170}
171
172/// Reads a value out of a matrix using indices already decoded by `triplet_indices`.
173#[inline]
174pub fn lookup_by_index(ixs: &[usize; TRIPLET_SIZE], matrix: &NucMatrix) -> f64 {
175    matrix[ixs[0]][ixs[1]][ixs[2]]
176}
177
178/// Looks up a value in a nucleotide matrix based on a triplet of nucleotides.
179///
180/// This is the convenience form taking an arbitrary slice. The hot path decodes once with
181/// `triplet_indices` and then calls `lookup_by_index` per matrix instead.
182///
183/// # Arguments
184///
185/// * `triplet` - A slice of u8 representing a triplet of nucleotides. Each u8 should be the ASCII
186///   value of 'A', 'C', 'G', or 'T', upper or lower case.
187/// * `matrix` - A reference to a `NucMatrix` to look up the value in.
188///
189/// # Returns
190///
191/// If the triplet is of length 3 and every base is recognized, this function returns a `Result`
192/// containing the value at the corresponding position in the matrix. Otherwise it returns a
193/// `Result` containing a `MatrixLookupError`.
194///
195/// # Errors
196///
197/// Returns a `MatrixLookupError` if the triplet is not of length 3, or if it contains a base
198/// that is not A, C, G, or T. These are reported as distinct errors rather than being conflated.
199pub fn matrix_lookup(triplet: &[u8], matrix: &NucMatrix) -> Result<f64, MatrixLookupError> {
200    let triplet: &[u8; TRIPLET_SIZE] = triplet
201        .try_into()
202        .map_err(|_| MatrixLookupError::WrongLength(triplet.len()))?;
203    Ok(lookup_by_index(&triplet_indices(triplet)?, matrix))
204}
205
206#[cfg(test)]
207mod tests {
208    use approx::assert_relative_eq;
209
210    use super::*;
211
212    #[test]
213    fn test_spot_check_indexing() {
214        assert_relative_eq!(TWIST[0][0][0], 0.598647428, epsilon = 1e-4);
215        assert_relative_eq!(TWIST[1][1][1], 0.598647428, epsilon = 1e-4);
216        assert_relative_eq!(ROLL_ACTIVE[1][2][0], 7.7, epsilon = 1e-4);
217        assert_relative_eq!(
218            matrix_lookup(b"AAA", &TWIST).unwrap(),
219            0.598647428,
220            epsilon = 1e-4
221        );
222        assert_relative_eq!(
223            matrix_lookup(b"CCC", &TWIST).unwrap(),
224            0.598647428,
225            epsilon = 1e-4
226        );
227        assert_relative_eq!(
228            matrix_lookup(b"CCA", &ROLL_SIMPLE).unwrap(),
229            0.7,
230            epsilon = 1e-4
231        );
232        assert!(matrix_lookup(b"AA", &ROLL_ACTIVE).is_err());
233        assert!(matrix_lookup(b"AAAA", &ROLL_ACTIVE).is_err());
234        assert!(matrix_lookup(b"AAN", &ROLL_ACTIVE).is_err());
235    }
236
237    #[test]
238    fn test_nuc_index_over_the_whole_byte_range() {
239        // The decode table is written out by hand, so check every possible byte
240        // rather than a sample: a wrong or missing entry is otherwise easy to miss.
241        for byte in 0u8..=255 {
242            let expected = match byte {
243                b'A' | b'a' => Some(0),
244                b'T' | b't' => Some(1),
245                b'G' | b'g' => Some(2),
246                b'C' | b'c' => Some(3),
247                _ => None,
248            };
249            assert_eq!(
250                nuc_index(byte),
251                expected,
252                "byte {byte:?} ({:?})",
253                byte as char
254            );
255        }
256    }
257
258    #[test]
259    fn test_triplet_indices_matches_per_base_decoding() {
260        let triplet = *b"CgA";
261        assert_eq!(triplet_indices(&triplet).unwrap(), [3, 2, 0]);
262        // Decoding once must agree with indexing the matrix the long way.
263        assert_relative_eq!(
264            lookup_by_index(&triplet_indices(&triplet).unwrap(), &ROLL_SIMPLE),
265            matrix_lookup(&triplet, &ROLL_SIMPLE).unwrap(),
266            epsilon = 1e-12
267        );
268        assert!(triplet_indices(b"CgN").is_err());
269    }
270
271    #[test]
272    fn test_matrix_lookup_error_display() {
273        assert_eq!(
274            MatrixLookupError::WrongLength(2).to_string(),
275            "triplet must be of length 3, got 2"
276        );
277        assert_eq!(
278            MatrixLookupError::UnknownBase(b'N').to_string(),
279            "unrecognized nucleotide 'N'"
280        );
281    }
282
283    #[test]
284    fn test_matrix_lookup_error_is_a_std_error() {
285        // So it composes with `?` into Box<dyn Error>, anyhow and the like.
286        fn takes_error<E: std::error::Error>(_: E) {}
287        takes_error(MatrixLookupError::UnknownBase(b'N'));
288        let boxed: Box<dyn std::error::Error> = MatrixLookupError::WrongLength(4).into();
289        assert!(boxed.to_string().contains("length 3"));
290    }
291}
symcurve/curve/matrix.rs

symcurve/curve/
matrix.rs
