zoe 0.0.30

A nightly library for viral genomics
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

use crate::{
    DEFAULT_SIMD_LANES,
    data::types::amino_acids::AminoAcidsReadable,
    distance::{
        DistanceError::{self, NoData, NotComparable},
        hamming_simd,
    },
    private::Sealed,
};

/// An extension trait adding amino acid distance methods to types containing
/// amino acid data.
pub trait AminoAcidsDistance: AminoAcidsReadable + Sealed {
    /// Calculates hamming distance between `self` and another sequence.
    ///
    /// ## Example
    ///
    /// ```
    /// # use zoe::{
    /// #     data::types::amino_acids::{AminoAcids, AminoAcidsView},
    /// #     distance::aa::AminoAcidsDistance
    /// # };
    ///
    /// let s1: AminoAcids = b"MANATEEMANATEEMANATEE".into();
    /// let s2: AminoAcids = b"MANGAEEMANATEEMANGAEE".into();
    ///
    /// assert!(4 == s1.distance_hamming(&s2));
    ///
    /// let s3 = AminoAcidsView::from_bytes_unchecked(b"MANGAEEMANATEEMANGAEE");
    /// assert!(4 == s1.distance_hamming(&s3));
    /// ```
    #[inline]
    #[must_use]
    fn distance_hamming<T: AminoAcidsReadable>(&self, other_sequence: &T) -> usize {
        hamming_simd::<{ DEFAULT_SIMD_LANES }>(self.amino_acids_bytes(), other_sequence.amino_acids_bytes())
    }

    /// Calculates physiochemical distance between `self`and another protein
    /// sequence. See: [`crate::distance::aa::physiochemical`]
    ///
    /// ## Citation
    ///
    /// For factor analysis used by the function:
    ///
    /// > Atchley et al. 2008. "Solving the protein sequence metric problem."
    /// > Proc Natl Acad Sci U S A. 2005 May 3;102(18):6395-400. Published 2005
    /// > Apr 25.
    ///
    ///
    /// ## Errors
    ///
    /// If either argument is empty or the sequence characters are invalid, an
    /// error is thrown. See [`DistanceError`].
    ///
    /// ## Example
    ///
    /// ```
    /// # use zoe::{
    /// #     data::types::amino_acids::{AminoAcids, AminoAcidsView},
    /// #     distance::aa::AminoAcidsDistance
    /// # };
    ///
    /// let s1: AminoAcids = b"MANATEEMANATEEMANATEE".into();
    /// let s2: AminoAcids = b"MANGAEEMANATEEMANGAEE".into();
    ///
    /// assert!( (s1.distance_physiochemical(&s2).unwrap() - 0.7643302).abs() < 0.01 );
    ///
    /// let s3 = AminoAcidsView::from_bytes_unchecked(b"MANGAEEMANATEEMANGAEE");
    /// assert!( (s1.distance_physiochemical(&s3).unwrap() - 0.7643302).abs() < 0.01 );
    /// ```
    #[inline]
    fn distance_physiochemical<T: AminoAcidsReadable>(&self, other_sequence: &T) -> Result<f32, DistanceError> {
        physiochemical(self.amino_acids_bytes(), other_sequence.amino_acids_bytes())
    }
}

impl<T: AminoAcidsReadable + Sealed> AminoAcidsDistance for T {}

/// Calculates a "physiochemical" distance measure using Euclidean distances
/// over physiochemical factors, normalized by sequence length. Only valid amino
/// acids are counted in the denominator.
///
/// ## Example
///
/// ```
/// # use zoe::distance::aa::physiochemical;
///
/// let s1 = b"MANATEE";
/// let s2 = b"MANGAEE";
///
/// assert!( (physiochemical(s1,s2).unwrap()  - 1.1464952).abs() < 0.01 )
/// ```
///
/// ## Citation
///
/// For factor analysis used by the function:
///
/// 1. Atchley et al. 2008. "Solving the protein sequence metric problem." Proc
///    Natl Acad Sci U S A. 2005 May 3;102(18):6395-400. Published 2005 Apr 25.
///
///
/// ## Errors
///
/// If either argument is empty or the sequence characters are all invalid, a
/// [`DistanceError`] is thrown.
//
// TODO: Could make a generic function that depends on distance matrix.
#[allow(clippy::cast_precision_loss)]
pub fn physiochemical(seq1: &[u8], seq2: &[u8]) -> Result<f32, DistanceError> {
    use crate::data::matrices::PHYSIOCHEMICAL_FACTORS as dm;

    if seq1.is_empty() || seq2.is_empty() {
        return Err(NoData);
    }

    if seq1 == seq2 {
        return if seq1.iter().any(|&aa| dm[b'A' as usize][aa as usize].is_some()) {
            Ok(0.0)
        } else {
            Err(NotComparable)
        };
    }

    let (number_valid, mut pcd_distance) = seq1
        .iter()
        .zip(seq2)
        .filter_map(|(a1, a2)| dm[*a1 as usize][*a2 as usize])
        .map(|val| (1, val))
        .fold((0u32, 0.0), |(a_valid, a_value), (i_valid, i_value)| {
            (a_valid + i_valid, a_value + i_value)
        });

    if number_valid > 0 {
        pcd_distance /= number_valid as f32;
        Ok(pcd_distance)
    } else {
        Err(NotComparable)
    }
}

#[cfg(test)]
mod test {
    use super::*;
    use crate::{assert_fp_eq, distance::DistanceError::NotComparable};

    #[test]
    fn pcd_identical_sequences() {
        assert!(matches!(physiochemical(b"???", b"???"), Err(NotComparable)));
        assert_fp_eq!(physiochemical(b"MANATEE", b"MANATEE").unwrap(), 0.0);
        assert_fp_eq!(physiochemical(b"MA?A", b"MA?A").unwrap(), 0.0);
    }
}