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.inc"
32
33 /*
34 #################################
35 # PRIVATE FUNCTION DECLARATIONS #
36 #################################
37 */
38
39 PRIVATE int
40 eval_hp_loop_fake(vrna_fold_compound_t *fc,
41 int i,
42 int j);
43
44
45 /*
46 #################################
47 # BEGIN OF FUNCTION DEFINITIONS #
48 #################################
49 */
50
51 /**
52 * @brief Evaluate the free energy of a hairpin loop
53 * and consider possible hard constraints
54 *
55 * @note This function is polymorphic! The provided #vrna_fold_compound_t may be of type
56 * #VRNA_FC_TYPE_SINGLE or #VRNA_FC_TYPE_COMPARATIVE
57 *
58 */
59 PUBLIC int
vrna_E_hp_loop(vrna_fold_compound_t * fc,int i,int j)60 vrna_E_hp_loop(vrna_fold_compound_t *fc,
61 int i,
62 int j)
63 {
64 vrna_callback_hc_evaluate *evaluate;
65 struct hc_hp_def_dat hc_dat_local;
66
67 if (fc->hc->type == VRNA_HC_WINDOW)
68 evaluate = prepare_hc_hp_def_window(fc, &hc_dat_local);
69 else
70 evaluate = prepare_hc_hp_def(fc, &hc_dat_local);
71
72 if ((i > 0) && (j > 0)) {
73 /* is this base pair allowed to close a hairpin (like) loop ? */
74 if (evaluate(i, j, i, j, VRNA_DECOMP_PAIR_HP, &hc_dat_local)) {
75 if (j > i) /* linear case */
76 return vrna_eval_hp_loop(fc, i, j);
77 else /* circular case */
78 return vrna_eval_ext_hp_loop(fc, j, i);
79 }
80 }
81
82 return INF;
83 }
84
85
86 /**
87 * @brief Evaluate the free energy of an exterior hairpin loop
88 * and consider possible hard constraints
89 */
90 PUBLIC int
vrna_E_ext_hp_loop(vrna_fold_compound_t * fc,int i,int j)91 vrna_E_ext_hp_loop(vrna_fold_compound_t *fc,
92 int i,
93 int j)
94 {
95 return vrna_E_hp_loop(fc, j, i);
96 }
97
98
99 /**
100 * @brief Evaluate free energy of an exterior hairpin loop
101 *
102 * @ingroup eval
103 *
104 */
105 PUBLIC int
vrna_eval_ext_hp_loop(vrna_fold_compound_t * fc,int i,int j)106 vrna_eval_ext_hp_loop(vrna_fold_compound_t *fc,
107 int i,
108 int j)
109 {
110 char **Ss, loopseq[10] = {
111 0
112 };
113 unsigned int **a2s;
114 short *S, *S2, **SS, **S5, **S3;
115 int u1, u2, e, s, type, n_seq, length, noGUclosure;
116 vrna_param_t *P;
117 vrna_md_t *md;
118 struct sc_hp_dat sc_wrapper;
119
120 length = fc->length;
121 P = fc->params;
122 md = &(P->model_details);
123 noGUclosure = md->noGUclosure;
124 e = INF;
125
126 init_sc_hp(fc, &sc_wrapper);
127
128 u1 = length - j;
129 u2 = i - 1;
130
131 if ((u1 + u2) < 3)
132 return e;
133
134 switch (fc->type) {
135 /* single sequences and cofolding hybrids */
136 case VRNA_FC_TYPE_SINGLE:
137 S = fc->sequence_encoding;
138 S2 = fc->sequence_encoding2;
139 type = vrna_get_ptype_md(S2[j], S2[i], md);
140
141 if (noGUclosure && ((type == 3) || (type == 4)))
142 break;
143
144 /* maximum special hp loop size: 6 */
145 if ((u1 + u2) < 7) {
146 memcpy(loopseq, fc->sequence + j - 1, sizeof(char) * (u1 + 1));
147 memcpy(loopseq + u1 + 1, fc->sequence, sizeof(char) * (u2 + 1));
148 loopseq[u1 + u2 + 2] = '\0';
149 }
150
151 e = E_Hairpin(u1 + u2, type, S[j + 1], S[i - 1], loopseq, P);
152
153 break;
154
155 /* sequence alignments */
156 case VRNA_FC_TYPE_COMPARATIVE:
157 SS = fc->S;
158 S5 = fc->S5; /* S5[s][i] holds next base 5' of i in sequence s */
159 S3 = fc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */
160 Ss = fc->Ss;
161 a2s = fc->a2s;
162 n_seq = fc->n_seq;
163 e = 0;
164
165 for (s = 0; s < n_seq; s++) {
166 u1 = a2s[s][length] - a2s[s][j];
167 u2 = a2s[s][i - 1];
168 memset(loopseq, '\0', sizeof(loopseq));
169
170 if ((u1 + u2) < 7) {
171 memcpy(loopseq, Ss[s] + a2s[s][j] - 1, sizeof(char) * (u1 + 1));
172 memcpy(loopseq + u1 + 1, Ss[s], sizeof(char) * (u2 + 1));
173 loopseq[u1 + u2 + 2] = '\0';
174 }
175
176 if ((u1 + u2) < 3) {
177 e += 600;
178 } else {
179 type = vrna_get_ptype_md(SS[s][j], SS[s][i], md);
180 e += E_Hairpin(u1 + u2, type, S3[s][j], S5[s][i], loopseq, P);
181 }
182 }
183
184 break;
185
186 /* nothing */
187 default:
188 break;
189 }
190
191 if (e != INF)
192 if (sc_wrapper.pair_ext)
193 e += sc_wrapper.pair_ext(i, j, &sc_wrapper);
194
195 free_sc_hp(&sc_wrapper);
196
197 return e;
198 }
199
200
201 /**
202 * @brief Evaluate free energy of a hairpin loop
203 *
204 * @ingroup eval
205 *
206 * @note This function is polymorphic! The provided #vrna_fold_compound_t may be of type
207 * #VRNA_FC_TYPE_SINGLE or #VRNA_FC_TYPE_COMPARATIVE
208 *
209 * @param fc The #vrna_fold_compound_t for the particular energy evaluation
210 * @param i 5'-position of the base pair
211 * @param j 3'-position of the base pair
212 * @returns Free energy of the hairpin loop closed by @f$ (i,j) @f$ in deka-kal/mol
213 */
214 PUBLIC int
vrna_eval_hp_loop(vrna_fold_compound_t * fc,int i,int j)215 vrna_eval_hp_loop(vrna_fold_compound_t *fc,
216 int i,
217 int j)
218 {
219 char **Ss;
220 unsigned int **a2s;
221 short *S, *S2, **SS, **S5, **S3;
222 unsigned int *sn;
223 int u, e, s, type, n_seq, en, noGUclosure;
224 vrna_param_t *P;
225 vrna_md_t *md;
226 vrna_ud_t *domains_up;
227 struct sc_hp_dat sc_wrapper;
228
229 P = fc->params;
230 md = &(P->model_details);
231 noGUclosure = md->noGUclosure;
232 sn = fc->strand_number;
233 domains_up = fc->domains_up;
234 e = INF;
235
236 if (sn[j] != sn[i])
237 return eval_hp_loop_fake(fc, i, j);
238
239 init_sc_hp(fc, &sc_wrapper);
240
241 /* regular hairpin loop */
242 switch (fc->type) {
243 /* single sequences and cofolding hybrids */
244 case VRNA_FC_TYPE_SINGLE:
245 S = fc->sequence_encoding;
246 S2 = fc->sequence_encoding2;
247 u = j - i - 1;
248 type = vrna_get_ptype_md(S2[i], S2[j], md);
249
250 if (noGUclosure && ((type == 3) || (type == 4)))
251 break;
252
253 e = E_Hairpin(u, type, S[i + 1], S[j - 1], fc->sequence + i - 1, P);
254
255 break;
256
257 /* sequence alignments */
258 case VRNA_FC_TYPE_COMPARATIVE:
259 SS = fc->S;
260 S5 = fc->S5; /* S5[s][i] holds next base 5' of i in sequence s */
261 S3 = fc->S3; /* Sl[s][i] holds next base 3' of i in sequence s */
262 Ss = fc->Ss;
263 a2s = fc->a2s;
264 n_seq = fc->n_seq;
265
266 for (e = s = 0; s < n_seq; s++) {
267 u = a2s[s][j - 1] - a2s[s][i];
268 if (u < 3) {
269 e += 600; /* ??? really 600 ??? */
270 } else {
271 type = vrna_get_ptype_md(SS[s][i], SS[s][j], md);
272 e += E_Hairpin(u, type, S3[s][i], S5[s][j], Ss[s] + (a2s[s][i - 1]), P);
273 }
274 }
275
276 break;
277
278 /* nothing */
279 default:
280 break;
281 }
282
283 if (e != INF) {
284 if (sc_wrapper.pair)
285 e += sc_wrapper.pair(i, j, &sc_wrapper);
286
287 /* consider possible ligand binding */
288 if (domains_up && domains_up->energy_cb) {
289 en = domains_up->energy_cb(fc,
290 i + 1, j - 1,
291 VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP,
292 domains_up->data);
293 if (en != INF)
294 en += e;
295
296 e = MIN2(e, en);
297 }
298 }
299
300 free_sc_hp(&sc_wrapper);
301
302 return e;
303 }
304
305
306 PRIVATE int
eval_hp_loop_fake(vrna_fold_compound_t * fc,int i,int j)307 eval_hp_loop_fake(vrna_fold_compound_t *fc,
308 int i,
309 int j)
310 {
311 short *S, *S2;
312 unsigned int *sn;
313 int u, e, ij, type, *idx, en, noGUclosure;
314 vrna_param_t *P;
315 vrna_sc_t *sc;
316 vrna_md_t *md;
317 vrna_ud_t *domains_up;
318
319 idx = fc->jindx;
320 P = fc->params;
321 md = &(P->model_details);
322 noGUclosure = md->noGUclosure;
323 sn = fc->strand_number;
324 domains_up = fc->domains_up;
325 e = INF;
326
327 switch (fc->type) {
328 /* single sequences and cofolding hybrids */
329 case VRNA_FC_TYPE_SINGLE:
330 S = fc->sequence_encoding;
331 S2 = fc->sequence_encoding2;
332 sc = fc->sc;
333 u = j - i - 1;
334 ij = idx[j] + i;
335 type = vrna_get_ptype_md(S2[j], S2[i], md);
336
337 if (noGUclosure && ((type == 3) || (type == 4)))
338 break;
339
340 /* hairpin-like exterior loop (for cofolding) */
341 short si, sj;
342
343 si = (sn[i + 1] == sn[i]) ? S[i + 1] : -1;
344 sj = (sn[j] == sn[j - 1]) ? S[j - 1] : -1;
345
346 switch (md->dangles) {
347 case 0:
348 e = vrna_E_ext_stem(type, -1, -1, P);
349 break;
350
351 case 2:
352 e = vrna_E_ext_stem(type, sj, si, P);
353 break;
354
355 default:
356 e = vrna_E_ext_stem(type, -1, -1, P);
357 e = MIN2(e, vrna_E_ext_stem(type, sj, -1, P));
358 e = MIN2(e, vrna_E_ext_stem(type, -1, si, P));
359 e = MIN2(e, vrna_E_ext_stem(type, sj, si, P));
360 break;
361 }
362
363 /* add soft constraints */
364 if (sc) {
365 if (sc->energy_up)
366 e += sc->energy_up[i + 1][u];
367
368 if (sc->energy_bp)
369 e += sc->energy_bp[ij];
370
371 if (sc->f)
372 e += sc->f(i, j, i, j, VRNA_DECOMP_PAIR_HP, sc->data);
373 }
374
375 /* consider possible ligand binding */
376 if (domains_up && domains_up->energy_cb) {
377 en = domains_up->energy_cb(fc,
378 i + 1, j - 1,
379 VRNA_UNSTRUCTURED_DOMAIN_HP_LOOP,
380 domains_up->data);
381 if (en != INF)
382 en += e;
383
384 e = MIN2(e, en);
385 }
386
387 break;
388
389 /* nothing */
390 default:
391 break;
392 }
393
394 return e;
395 }
396