nale 0.1.2

Nale is a library that performs profile Hidden Markov Model (PHMM) biological sequence alignment.
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

use crate::structs::dp_matrix::DpMatrix;
use crate::structs::Profile;

pub fn posterior(
    profile: &Profile,
    forward_matrix: &impl DpMatrix,
    backward_matrix: &impl DpMatrix,
    posterior_matrix: &mut impl DpMatrix,
) {
    let target_length = forward_matrix.target_length();
    let overall_score: f32 = forward_matrix.get_special(target_length, Profile::SPECIAL_C_IDX)
        + profile.special_transition_score(Profile::SPECIAL_C_IDX, Profile::SPECIAL_MOVE_IDX);
    let mut denominator: f32;

    posterior_matrix.set_special(0, Profile::SPECIAL_E_IDX, 0.0);
    posterior_matrix.set_special(0, Profile::SPECIAL_N_IDX, 0.0);
    posterior_matrix.set_special(0, Profile::SPECIAL_J_IDX, 0.0);
    posterior_matrix.set_special(0, Profile::SPECIAL_B_IDX, 0.0);
    posterior_matrix.set_special(0, Profile::SPECIAL_C_IDX, 0.0);

    for k in 0..=profile.length {
        posterior_matrix.set_match(0, k, 0.0);
        posterior_matrix.set_insert(0, k, 0.0);
        posterior_matrix.set_delete(0, k, 0.0);
    }

    for i in 1..=target_length {
        denominator = 0.0;
        posterior_matrix.set_match(i, 0, 0.0);
        posterior_matrix.set_insert(i, 0, 0.0);
        posterior_matrix.set_delete(i, 0, 0.0);

        for k in 1..profile.length {
            posterior_matrix.set_match(
                i,
                k,
                (forward_matrix.get_match(i, k) + backward_matrix.get_match(i, k) - overall_score)
                    .exp(),
            );

            denominator += posterior_matrix.get_match(i, k);

            posterior_matrix.set_insert(
                i,
                k,
                (forward_matrix.get_insert(i, k) + backward_matrix.get_insert(i, k)
                    - overall_score)
                    .exp(),
            );
            denominator += posterior_matrix.get_insert(i, k);

            posterior_matrix.set_delete(i, k, 0.0);
        }

        posterior_matrix.set_match(
            i,
            profile.length,
            (forward_matrix.get_match(i, profile.length)
                + backward_matrix.get_match(i, profile.length)
                - overall_score)
                .exp(),
        );

        denominator += posterior_matrix.get_match(i, profile.length);
        posterior_matrix.set_insert(i, profile.length, 0.0);
        posterior_matrix.set_delete(i, profile.length, 0.0);

        posterior_matrix.set_special(i, Profile::SPECIAL_E_IDX, 0.0);
        posterior_matrix.set_special(
            i,
            Profile::SPECIAL_N_IDX,
            (forward_matrix.get_special(i - 1, Profile::SPECIAL_N_IDX)
                + backward_matrix.get_special(i, Profile::SPECIAL_N_IDX)
                + profile
                    .special_transition_score(Profile::SPECIAL_N_IDX, Profile::SPECIAL_LOOP_IDX)
                - overall_score)
                .exp(),
        );

        posterior_matrix.set_special(
            i,
            Profile::SPECIAL_J_IDX,
            (forward_matrix.get_special(i - 1, Profile::SPECIAL_J_IDX)
                + backward_matrix.get_special(i, Profile::SPECIAL_J_IDX)
                + profile
                    .special_transition_score(Profile::SPECIAL_J_IDX, Profile::SPECIAL_LOOP_IDX)
                - overall_score)
                .exp(),
        );

        posterior_matrix.set_special(i, Profile::SPECIAL_B_IDX, 0.0);

        posterior_matrix.set_special(
            i,
            Profile::SPECIAL_C_IDX,
            (forward_matrix.get_special(i - 1, Profile::SPECIAL_C_IDX)
                + backward_matrix.get_special(i, Profile::SPECIAL_C_IDX)
                + profile
                    .special_transition_score(Profile::SPECIAL_C_IDX, Profile::SPECIAL_LOOP_IDX)
                - overall_score)
                .exp(),
        );

        denominator += posterior_matrix.get_special(i, Profile::SPECIAL_N_IDX);
        denominator += posterior_matrix.get_special(i, Profile::SPECIAL_J_IDX);
        denominator += posterior_matrix.get_special(i, Profile::SPECIAL_C_IDX);

        denominator = 1.0 / denominator;

        for k in 1..profile.length {
            posterior_matrix.set_match(i, k, posterior_matrix.get_match(i, k) * denominator);
            posterior_matrix.set_insert(i, k, posterior_matrix.get_insert(i, k) * denominator);
        }
        posterior_matrix.set_match(
            i,
            profile.length,
            posterior_matrix.get_match(i, profile.length) * denominator,
        );
        posterior_matrix.set_special(
            i,
            Profile::SPECIAL_N_IDX,
            posterior_matrix.get_special(i, Profile::SPECIAL_N_IDX) * denominator,
        );
        posterior_matrix.set_special(
            i,
            Profile::SPECIAL_J_IDX,
            posterior_matrix.get_special(i, Profile::SPECIAL_J_IDX) * denominator,
        );
        posterior_matrix.set_special(
            i,
            Profile::SPECIAL_C_IDX,
            posterior_matrix.get_special(i, Profile::SPECIAL_C_IDX) * denominator,
        );
    }
}