zoe 0.0.31

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

use crate::data::byte_types::ByteMappings;

/// Generates a 4x4 substitution matrix from two aligned sequences.
///
/// Each column refers to nucleotide bases in the first sequence (`A`, `C`, `G`,
/// and `T`, respectively), and each row refers to nucleotide bases in the
/// second sequence.
///
/// The substitution matrix can then be indexed to find counts of each type of
/// substitution between the sequences, or used to calculate evolutionary
/// distances.
///
/// ## Example
///
/// ```
/// # use zoe::distance::dna::dna_substitution_matrix;
/// let seq1: &[u8] = b"GATCAGATTTGCATTGGTT";
/// let seq2: &[u8] = b"GATCATATTAGCATTGCTT";
///
/// let sub_matrix: [[u32; 4]; 4] = dna_substitution_matrix(seq1, seq2);
/// # assert_eq!(sub_matrix, [[4, 0, 0, 1], [0, 2, 1, 0], [0, 0, 3, 0], [0, 0, 1, 7]])
/// ```
/// This should return the matrix (bases provided here for context):
///
/// $$ \begin{pmatrix}  & \text{A} & \text{C} & \text{G} & \text{T} \cr \text{A}
/// & 4 & 0 & 0 & 1 \cr \text{C} & 0 & 2 & 1 & 0 \cr \text{G} & 0 & 0 & 3 & 0
/// \cr \text{T} & 0 & 0 & 1 & 7 \end{pmatrix} $$
#[must_use]
pub fn dna_substitution_matrix(seq1: &[u8], seq2: &[u8]) -> [[u32; 4]; 4] {
    let mut sub_matrix = [[0u32; 4]; 4];
    std::iter::zip(
        seq1.iter().copied().map(ByteMappings::to_dna_index),
        seq2.iter().copied().map(ByteMappings::to_dna_index),
    )
    .filter(|(a, b)| *a < 4 && *b < 4)
    .for_each(|(a, b)| sub_matrix[b][a] += 1);
    sub_matrix
}

/// Calculates the Hamming distance from a nucleotide [substitution matrix](dna_substitution_matrix).
///
/// The Hamming distance is the number of substitutions between aligned
/// sequences, or sum of non-diagonal values in the substitution matrix.
#[allow(clippy::needless_range_loop)]
#[inline]
#[must_use]
pub(crate) fn hamming_dist_from_sub_matrix(sub_matrix: &[[u32; 4]; 4]) -> u32 {
    let mut hamming = 0;

    for i in 0..4 {
        for j in 0..4 {
            if i != j {
                hamming += sub_matrix[i][j];
            }
        }
    }
    hamming
}

/// Finds total number of bases and nucleotide base frequencies from a
/// nucleotide [substitution matrix](dna_substitution_matrix).
///
/// Returns a tuple of the total number of nucleotides in the two sequences, and
/// an array of base frequencies, corresponding to the frequencies of A, C, G,
/// and T bases in the sequences, respectively.
#[must_use]
pub(crate) fn total_and_frequencies(sub_matrix: &[[u32; 4]; 4]) -> (u32, [f64; 4]) {
    let mut bf = [0f64; 4];
    for i in 0..4 {
        // sum of ith column and row to get base frequencies
        bf[i] = f64::from(sub_matrix[i].iter().sum::<u32>() + sub_matrix.iter().map(|row| row[i]).sum::<u32>());
    }
    let total_bases = sub_matrix.iter().flatten().sum::<u32>();
    // counts of nucleotides are divided by 2 to remove double counting
    for freq in &mut bf {
        *freq /= 2.0 * f64::from(total_bases);
    }

    (total_bases, bf)
}