1 /* p7 model configuration for defining envelopes for truncated hits:
2 *
3 * Contents:
4 * 1. Routines in the exposed API.
5 *
6 * EPN, Mon Nov 21 10:02:03 2011
7 */
8 #include "esl_config.h"
9 #include "p7_config.h"
10 #include "config.h"
11
12 #include <math.h>
13 #include <float.h>
14 #include <string.h>
15 #include <ctype.h>
16
17 #include "easel.h"
18 #include "esl_vectorops.h"
19
20 #include "hmmer.h"
21
22 #include "infernal.h"
23
24 /*****************************************************************
25 * 1. Routines in the exposed API.
26 *****************************************************************/
27
28 /* Function: p7_ProfileConfig5PrimeTrunc()
29 */
30 int
p7_ProfileConfig5PrimeTrunc(P7_PROFILE * gm,int L)31 p7_ProfileConfig5PrimeTrunc(P7_PROFILE *gm, int L)
32 {
33 int status;
34 int k;
35
36 /* we should be in glocal mode */
37 assert(! p7_IsLocal(gm->mode));
38
39 /* Local mode entry, uniform: 1/M */
40 for (k = 1; k <= gm->M; k++)
41 p7P_TSC(gm, k-1, p7P_BM) = log(1. / gm->M); /* note off-by-one: entry at Mk stored as [k-1][BM] */
42
43 /* set profile's mode to UNIGLOCAL, even though we're really a hybrid of local/glocal
44 * We do this so p7_GForward, p7_GBackward set end scores to -eslINFINITY
45 * forcing us to only exit from last node of the model
46 */
47 gm->mode = p7_UNIGLOCAL;
48
49 /* borrowed and modified from p7_ReconfigUnihit(),
50 * make J unreachable (unihit) and N->N transitions impossible (special to this func)
51 */
52 gm->xsc[p7P_N][p7P_MOVE] = gm->xsc[p7P_E][p7P_MOVE] = 0.0f;
53 gm->xsc[p7P_N][p7P_LOOP] = gm->xsc[p7P_E][p7P_LOOP] = -eslINFINITY;
54
55 /* Remaining specials, [C][MOVE | LOOP] are set by ReconfigLength5PrimeTrunc()
56 */
57 gm->L = 0; /* force ReconfigLength to reconfig */
58 if ((status = p7_ReconfigLength5PrimeTrunc(gm, L)) != eslOK) return status;
59
60 return eslOK;
61 }
62
63 /* Function: p7_ProfileConfig3PrimeTrunc()
64 */
65 int
p7_ProfileConfig3PrimeTrunc(const P7_HMM * hmm,P7_PROFILE * gm,int L)66 p7_ProfileConfig3PrimeTrunc(const P7_HMM *hmm, P7_PROFILE *gm, int L)
67 {
68 int status;
69 int k;
70 float Z;
71
72 /* we should be in glocal mode */
73 assert(! p7_IsLocal(gm->mode));
74
75 /* glocal mode entry: left wing retraction, must be in log space for precision
76 * (copied from HMMER's modelconfig.c:p7_ProfileConfig()) */
77 Z = log(hmm->t[0][p7H_MD]);
78 p7P_TSC(gm, 0, p7P_BM) = log(1.0 - hmm->t[0][p7H_MD]);
79 for (k = 1; k < hmm->M; k++)
80 {
81 p7P_TSC(gm, k, p7P_BM) = Z + log(hmm->t[k][p7H_DM]);
82 Z += log(hmm->t[k][p7H_DD]);
83 }
84 /* set profile's mode to UNILOCAL, even though we're really a hybrid of local/glocal
85 * We do this so p7_GForward, p7_GBackward do NOT set end scores to -eslINFINITY
86 * this allows us to exit from any model node (like local ends in a CP9)
87 */
88 gm->mode = p7_UNILOCAL;
89
90 /* borrowed and modified from p7_ReconfigUnihit(),
91 * make J unreachable (unihit) and C->C transitions impossible (special to this func)
92 */
93 gm->xsc[p7P_C][p7P_MOVE] = gm->xsc[p7P_E][p7P_MOVE] = 0.0f;
94 gm->xsc[p7P_C][p7P_LOOP] = gm->xsc[p7P_E][p7P_LOOP] = -eslINFINITY;
95
96 /* Remaining specials, [N][MOVE | LOOP] are set by ReconfigLength5PrimeTrunc()
97 */
98 gm->L = 0; /* force ReconfigLength to reconfig */
99 if ((status = p7_ReconfigLength3PrimeTrunc(gm, L)) != eslOK) return status;
100
101 return eslOK;
102 }
103
104 /* Function: p7_ProfileConfig5PrimeAnd3PrimeTrunc()
105 */
106 int
p7_ProfileConfig5PrimeAnd3PrimeTrunc(P7_PROFILE * gm,int L)107 p7_ProfileConfig5PrimeAnd3PrimeTrunc(P7_PROFILE *gm, int L)
108 {
109 assert(gm->mode == p7_LOCAL);
110
111 /* gm should already be set up as local with equiprobable begins and ends,
112 * we just need to make it unihit, and make N->N and C->C transitions impossible
113 */
114 gm->xsc[p7P_N][p7P_MOVE] = gm->xsc[p7P_C][p7P_MOVE] = gm->xsc[p7P_E][p7P_MOVE] = 0.0f;
115 gm->xsc[p7P_N][p7P_LOOP] = gm->xsc[p7P_C][p7P_LOOP] = gm->xsc[p7P_E][p7P_LOOP] = -eslINFINITY;
116 gm->nj = 0.0f;
117 gm->L = L;
118 gm->mode = p7_UNILOCAL;
119 /* don't call ReconfigLength() */
120
121 return eslOK;
122 }
123
124 /* Function: p7_ReconfigLength5PrimeTrunc()
125 */
126 int
p7_ReconfigLength5PrimeTrunc(P7_PROFILE * gm,int L)127 p7_ReconfigLength5PrimeTrunc(P7_PROFILE *gm, int L)
128 {
129 float ploop, pmove;
130
131 /* mode should be set as p7_UNIGLOCAL for 5' trunc case, even though its
132 * actually in a hybrid local/glocal mode. It's impt it is p7_UNIGLOCAL
133 * so it is behaves appropriately in other HMMER functions.
134 */
135 assert(gm->mode == p7_UNIGLOCAL);
136
137 /* Configure C transitions so it bears the total
138 * unannotated sequence length L.
139 */
140 pmove = (1.0f + gm->nj) / ((float) L + 1.0f + gm->nj); /* 1/(L+1) */
141 ploop = 1.0f - pmove;
142 gm->xsc[p7P_C][p7P_LOOP] = gm->xsc[p7P_J][p7P_LOOP] = log(ploop); /* J should be irrelevant, b/c we've set unihit */
143 gm->xsc[p7P_C][p7P_MOVE] = gm->xsc[p7P_J][p7P_MOVE] = log(pmove);
144 gm->L = L;
145 return eslOK;
146 }
147
148
149 /* Function: p7_ReconfigLength3PrimeTrunc()
150 */
151 int
p7_ReconfigLength3PrimeTrunc(P7_PROFILE * gm,int L)152 p7_ReconfigLength3PrimeTrunc(P7_PROFILE *gm, int L)
153 {
154 float ploop, pmove;
155
156 /* mode should be set as p7_UNILOCAL for 5' trunc case, even though its
157 * actually in a hybrid local/glocal mode. It's impt it is p7_UNILOCAL
158 * so it is behaves appropriately in other HMMER functions.
159 */
160 assert(gm->mode == p7_UNILOCAL);
161
162 /* Configure N transitions so it bears the total
163 * unannotated sequence length L.
164 */
165 pmove = (1.0f + gm->nj) / ((float) L + 1.0f + gm->nj); /* 1/(L+1) */
166 ploop = 1.0f - pmove;
167 gm->xsc[p7P_N][p7P_LOOP] = gm->xsc[p7P_J][p7P_LOOP] = log(ploop); /* J should be irrelevant, b/c we've set unihit */
168 gm->xsc[p7P_N][p7P_MOVE] = gm->xsc[p7P_J][p7P_MOVE] = log(pmove);
169 gm->L = L;
170 return eslOK;
171 }
172