1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <math.h>
8 #include <ctype.h>
9 #include <string.h>
10 #include "ViennaRNA/fold_vars.h"
11 #include "ViennaRNA/params/default.h"
12 #include "ViennaRNA/datastructures/basic.h"
13 #include "ViennaRNA/params/basic.h"
14 #include "ViennaRNA/utils/basic.h"
15 #include "ViennaRNA/constraints/hard.h"
16 #include "ViennaRNA/constraints/soft.h"
17 #include "ViennaRNA/loops/external.h"
18 #include "ViennaRNA/gquad.h"
19 #include "ViennaRNA/structured_domains.h"
20 #include "ViennaRNA/unstructured_domains.h"
21 #include "ViennaRNA/alphabet.h"
22 #include "ViennaRNA/loops/hairpin.h"
23
24 #ifdef __GNUC__
25 # define INLINE inline
26 #else
27 # define INLINE
28 #endif
29
30 #include "hairpin_hc.inc"
31 #include "hairpin_sc_pf.inc"
32
33 /*
34 #################################
35 # PRIVATE FUNCTION DECLARATIONS #
36 #################################
37 */
38
39 PRIVATE FLT_OR_DBL
40 exp_eval_hp_loop(vrna_fold_compound_t *fc,
41 int i,
42 int j);
43
44
45 PRIVATE FLT_OR_DBL
46 exp_eval_ext_hp_loop(vrna_fold_compound_t *fc,
47 int i,
48 int j);
49
50
51 PRIVATE FLT_OR_DBL
52 exp_eval_hp_loop_fake(vrna_fold_compound_t *fc,
53 int i,
54 int j);
55
56
57 /*
58 #################################
59 # BEGIN OF FUNCTION DEFINITIONS #
60 #################################
61 */
62
63 /**
64 * @brief High-Level function for hairpin loop energy evaluation (partition function variant)
65 *
66 * @see E_hp_loop() for it's free energy counterpart
67 */
68 PUBLIC FLT_OR_DBL
vrna_exp_E_hp_loop(vrna_fold_compound_t * fc,int i,int j)69 vrna_exp_E_hp_loop(vrna_fold_compound_t *fc,
70 int i,
71 int j)
72 {
73 vrna_callback_hc_evaluate *evaluate;
74 struct hc_hp_def_dat hc_dat_local;
75
76 if (fc->hc->type == VRNA_HC_WINDOW)
77 evaluate = prepare_hc_hp_def_window(fc, &hc_dat_local);
78 else
79 evaluate = prepare_hc_hp_def(fc, &hc_dat_local);
80
81 if ((i > 0) && (j > 0)) {
82 if (evaluate(i, j, i, j, VRNA_DECOMP_PAIR_HP, &hc_dat_local)) {
83 if (j > i) /* linear case */
84 return exp_eval_hp_loop(fc, i, j);
85 else /* circular case */
86 return exp_eval_ext_hp_loop(fc, j, i);
87 }
88 }
89
90 return 0.;
91 }
92
93
94 PRIVATE FLT_OR_DBL
exp_eval_hp_loop_fake(vrna_fold_compound_t * fc,int i,int j)95 exp_eval_hp_loop_fake(vrna_fold_compound_t *fc,
96 int i,
97 int j)
98 {
99 short *S, *S2, s5, s3;
100 unsigned int *sn, *ss, *se;
101 int u, type, *iidx, *jidx;
102 FLT_OR_DBL qq, temp, *q, *scale;
103 vrna_exp_param_t *pf_params;
104 vrna_sc_t *sc;
105 vrna_md_t *md;
106 vrna_ud_t *domains_up;
107
108 iidx = fc->iindx;
109 jidx = fc->jindx;
110 pf_params = fc->exp_params;
111 md = &(pf_params->model_details);
112 q = fc->exp_matrices->q;
113 scale = fc->exp_matrices->scale;
114 sn = fc->strand_number;
115 ss = fc->strand_start;
116 se = fc->strand_end;
117 domains_up = fc->domains_up;
118
119 qq = 0;
120
121 switch (fc->type) {
122 /* single sequences and cofolding hybrids */
123 case VRNA_FC_TYPE_SINGLE:
124 S = fc->sequence_encoding;
125 S2 = fc->sequence_encoding2;
126 sc = fc->sc;
127 u = j - i - 1;
128 type = vrna_get_ptype_md(S2[j], S2[i], md);
129
130 temp = scale[2];
131
132 if (u > 0) {
133 /* add contribution for [i + 1, end(strand(i))] */
134 if (sn[i] == sn[i + 1])
135 temp *= q[iidx[i + 1] - se[sn[i]]];
136
137 /* add contribution for [start(strand(j)), j - 1] */
138 if (sn[j - 1] == sn[j])
139 temp *= q[iidx[ss[sn[j]]] - (j - 1)];
140 }
141
142 s5 = (sn[j] == sn[j - 1]) ? S[j - 1] : -1;
143 s3 = (sn[i + 1] == sn[i]) ? S[i + 1] : -1;
144
145 temp *= vrna_exp_E_ext_stem(type, s5, s3, pf_params);
146
147 qq += temp;
148
149 /* add soft constraints */
150 if (sc) {
151 if (sc->exp_energy_up)
152 qq *= sc->exp_energy_up[i + 1][u];
153
154 if (sc->exp_energy_bp)
155 qq *= sc->exp_energy_bp[jidx[j] + i];
156
157 if (sc->exp_f)
158 qq *= sc->exp_f(i, j, i, j, VRNA_DECOMP_PAIR_HP, sc->data);
159 }
160
161 if (domains_up && domains_up->exp_energy_cb) {
162 /* we always consider both, bound and unbound state */
163 qq += qq * domains_up->exp_energy_cb(fc,
164 i + 1, j - 1,
165 VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP,
166 domains_up->data);
167 }
168
169 break;
170
171 /* nothing */
172 default:
173 break;
174 }
175
176 return qq;
177 }
178
179
180 PRIVATE FLT_OR_DBL
exp_eval_hp_loop(vrna_fold_compound_t * fc,int i,int j)181 exp_eval_hp_loop(vrna_fold_compound_t *fc,
182 int i,
183 int j)
184 {
185 char **Ss;
186 unsigned int **a2s;
187 short *S, *S2, **SS, **S5, **S3;
188 unsigned int *sn;
189 int u, type, n_seq, s;
190 FLT_OR_DBL q, qbt1, *scale;
191 vrna_exp_param_t *P;
192 vrna_md_t *md;
193 vrna_ud_t *domains_up;
194 struct sc_hp_exp_dat sc_wrapper;
195
196 P = fc->exp_params;
197 md = &(P->model_details);
198 sn = fc->strand_number;
199 scale = fc->exp_matrices->scale;
200 domains_up = fc->domains_up;
201
202 init_sc_hp_exp(fc, &sc_wrapper);
203 q = 0.;
204
205 if (sn[j] != sn[i])
206 return exp_eval_hp_loop_fake(fc, i, j);
207
208 switch (fc->type) {
209 case VRNA_FC_TYPE_SINGLE:
210 S = fc->sequence_encoding;
211 S2 = fc->sequence_encoding2;
212 u = j - i - 1;
213 type = vrna_get_ptype_md(S2[i], S2[j], md);
214
215 if (sn[j] == sn[i]) {
216 /* regular hairpin loop */
217 q = exp_E_Hairpin(u, type, S[i + 1], S[j - 1], fc->sequence + i - 1, P);
218 } else {
219 /*
220 * hairpin-like exterior loop (for cofolding)
221 * this is currently handle somewhere else
222 */
223 }
224
225 break;
226
227 case VRNA_FC_TYPE_COMPARATIVE:
228 SS = fc->S;
229 S5 = fc->S5; /* S5[s][i] holds next base 5' of i in sequence s */
230 S3 = fc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */
231 Ss = fc->Ss;
232 a2s = fc->a2s;
233 n_seq = fc->n_seq;
234 qbt1 = 1.;
235
236 for (s = 0; s < n_seq; s++) {
237 u = a2s[s][j - 1] - a2s[s][i];
238 if (a2s[s][i] < 1)
239 continue;
240
241 type = vrna_get_ptype_md(SS[s][i], SS[s][j], md);
242 qbt1 *= exp_E_Hairpin(u, type, S3[s][i], S5[s][j], Ss[s] + a2s[s][i] - 1, P);
243 }
244
245 q = qbt1;
246 break;
247
248 default:
249 break;
250 }
251
252 /* add soft constraints */
253 if (sc_wrapper.pair)
254 q *= sc_wrapper.pair(i, j, &sc_wrapper);
255
256 if (domains_up && domains_up->exp_energy_cb) {
257 /* we always consider both, bound and unbound state */
258 q += q * domains_up->exp_energy_cb(fc,
259 i + 1, j - 1,
260 VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP,
261 domains_up->data);
262 }
263
264 q *= scale[j - i + 1];
265
266 free_sc_hp_exp(&sc_wrapper);
267
268 return q;
269 }
270
271
272 PRIVATE FLT_OR_DBL
exp_eval_ext_hp_loop(vrna_fold_compound_t * fc,int i,int j)273 exp_eval_ext_hp_loop(vrna_fold_compound_t *fc,
274 int i,
275 int j)
276 {
277 char **Ss, *sequence, loopseq[10] = {
278 0
279 };
280 unsigned int **a2s;
281 short *S, *S2, **SS, **S5, **S3;
282 int u1, u2, n, type, n_seq, s, noGUclosure;
283 FLT_OR_DBL q, qbt1, *scale;
284 vrna_exp_param_t *P;
285 vrna_md_t *md;
286 vrna_ud_t *domains_up;
287 struct sc_hp_exp_dat sc_wrapper;
288
289 n = fc->length;
290 P = fc->exp_params;
291 md = &(P->model_details);
292 noGUclosure = md->noGUclosure;
293 scale = fc->exp_matrices->scale;
294 domains_up = fc->domains_up;
295
296 init_sc_hp_exp(fc, &sc_wrapper);
297
298 q = 0.;
299 u1 = n - j;
300 u2 = i - 1;
301
302 if ((u1 + u2) < 3)
303 return q;
304
305 switch (fc->type) {
306 case VRNA_FC_TYPE_SINGLE:
307 sequence = fc->sequence;
308 S = fc->sequence_encoding;
309 S2 = fc->sequence_encoding2;
310 type = vrna_get_ptype_md(S2[j], S2[i], md);
311
312 if (((type == 3) || (type == 4)) && noGUclosure)
313 return q;
314
315 /* get the loop sequence */
316 if ((u1 + u2) < 7) {
317 memcpy(loopseq, sequence + j - 1, sizeof(char) * (u1 + 1));
318 memcpy(loopseq + u1 + 1, sequence, sizeof(char) * (u2 + 1));
319 loopseq[u1 + u2 + 2] = '\0';
320 }
321
322 q = exp_E_Hairpin(u1 + u2, type, S[j + 1], S[i - 1], loopseq, P);
323
324 break;
325
326 case VRNA_FC_TYPE_COMPARATIVE:
327 SS = fc->S;
328 S5 = fc->S5; /* S5[s][i] holds next base 5' of i in sequence s */
329 S3 = fc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */
330 Ss = fc->Ss;
331 a2s = fc->a2s;
332 n_seq = fc->n_seq;
333 qbt1 = 1.;
334
335 for (s = 0; s < n_seq; s++) {
336 const int u1_local = a2s[s][n] - a2s[s][j];
337 const int u2_local = a2s[s][i - 1];
338 memset(loopseq, '\0', sizeof(loopseq));
339
340 if ((u1_local + u2_local) < 7) {
341 memcpy(loopseq, Ss[s] + a2s[s][j] - 1, sizeof(char) * (u1_local + 1));
342 memcpy(loopseq + u1_local + 1, Ss[s], sizeof(char) * (u2_local + 1));
343 loopseq[u1_local + u2_local + 2] = '\0';
344 }
345
346 type = vrna_get_ptype_md(SS[s][j], SS[s][i], md);
347 qbt1 *= exp_E_Hairpin(u1_local + u2_local, type, S3[s][j], S5[s][i], loopseq, P);
348 }
349
350 q = qbt1;
351
352 break;
353
354 default:
355 break;
356 }
357
358 /* add soft constraints */
359 if (sc_wrapper.pair_ext)
360 q *= sc_wrapper.pair_ext(i, j, &sc_wrapper);
361
362 if (domains_up && domains_up->exp_energy_cb) {
363 /* we always consider both, bound and unbound state */
364 q += q * domains_up->exp_energy_cb(fc,
365 j + 1, i - 1,
366 VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP,
367 domains_up->data);
368 }
369
370 q *= scale[u1 + u2];
371
372 free_sc_hp_exp(&sc_wrapper);
373
374 return q;
375 }
376