1 /*
2 * partiton function for single RNA secondary structures
3 *
4 * Simplified interfaces and backward compatibility
5 * wrappers
6 *
7 * Ivo L Hofacker + Ronny Lorenz
8 * Vienna RNA package
9 */
10
11 #ifdef HAVE_CONFIG_H
12 #include "config.h"
13 #endif
14
15 /*###########################################*/
16 /*# deprecated functions below #*/
17 /*###########################################*/
18
19 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include <math.h>
25 #include <float.h> /* #defines FLT_MAX ... */
26 #include <limits.h>
27
28 #include "ViennaRNA/utils/basic.h"
29 #include "ViennaRNA/params/default.h"
30 #include "ViennaRNA/fold_vars.h"
31 #include "ViennaRNA/loops/all.h"
32 #include "ViennaRNA/gquad.h"
33 #include "ViennaRNA/constraints/hard.h"
34 #include "ViennaRNA/constraints/soft.h"
35 #include "ViennaRNA/mfe.h"
36 #include "ViennaRNA/part_func.h"
37
38 #ifdef _OPENMP
39 #include <omp.h>
40 #endif
41
42 /*
43 #################################
44 # GLOBAL VARIABLES #
45 #################################
46 */
47 PUBLIC int st_back = 0;
48
49 /*
50 #################################
51 # PRIVATE VARIABLES #
52 #################################
53 */
54
55 /* some backward compatibility stuff */
56 PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
57 PRIVATE int backward_compat = 0;
58
59 #ifdef _OPENMP
60
61 #pragma omp threadprivate(backward_compat_compound, backward_compat)
62
63 #endif
64
65 /*
66 #################################
67 # PRIVATE FUNCTION DECLARATIONS #
68 #################################
69 */
70
71 PRIVATE float
72 wrap_pf_fold(const char *sequence,
73 char *structure,
74 vrna_exp_param_t *parameters,
75 int calculate_bppm,
76 int is_constrained,
77 int is_circular);
78
79
80 PRIVATE double
81 wrap_mean_bp_distance(FLT_OR_DBL *p,
82 int length,
83 int *index,
84 int turn);
85
86 /*
87 #################################
88 # BEGIN OF FUNCTION DEFINITIONS #
89 #################################
90 */
91
92
93 PRIVATE double
wrap_mean_bp_distance(FLT_OR_DBL * p,int length,int * index,int turn)94 wrap_mean_bp_distance(FLT_OR_DBL *p,
95 int length,
96 int *index,
97 int turn)
98 {
99 int i, j;
100 double d = 0.;
101
102 /* compute the mean base pair distance in the thermodynamic ensemble */
103 /* <d> = \sum_{a,b} p_a p_b d(S_a,S_b)
104 * this can be computed from the pair probs p_ij as
105 * <d> = \sum_{ij} p_{ij}(1-p_{ij}) */
106
107 for (i = 1; i <= length; i++)
108 for (j = i + turn + 1; j <= length; j++)
109 d += p[index[i] - j] * (1 - p[index[i] - j]);
110
111 return 2 * d;
112 }
113
114
115 PRIVATE float
wrap_pf_fold(const char * sequence,char * structure,vrna_exp_param_t * parameters,int calculate_bppm,int is_constrained,int is_circular)116 wrap_pf_fold(const char *sequence,
117 char *structure,
118 vrna_exp_param_t *parameters,
119 int calculate_bppm,
120 int is_constrained,
121 int is_circular)
122 {
123 vrna_fold_compound_t *vc;
124 vrna_md_t md;
125
126 vc = NULL;
127
128 /* we need vrna_exp_param_t datastructure to correctly init default hard constraints */
129 if (parameters)
130 md = parameters->model_details;
131 else
132 set_model_details(&md); /* get global default parameters */
133
134 md.circ = is_circular;
135 md.compute_bpp = calculate_bppm;
136
137 vc = vrna_fold_compound(sequence, &md, VRNA_OPTION_DEFAULT);
138
139 /* prepare exp_params and set global pf_scale */
140 vc->exp_params = vrna_exp_params(&(vc->params->model_details));
141 vc->exp_params->pf_scale = pf_scale;
142
143 if (is_constrained && structure) {
144 unsigned int constraint_options = 0;
145 constraint_options |= VRNA_CONSTRAINT_DB
146 | VRNA_CONSTRAINT_DB_PIPE
147 | VRNA_CONSTRAINT_DB_DOT
148 | VRNA_CONSTRAINT_DB_X
149 | VRNA_CONSTRAINT_DB_ANG_BRACK
150 | VRNA_CONSTRAINT_DB_RND_BRACK;
151
152 vrna_constraints_add(vc, (const char *)structure, constraint_options);
153 }
154
155 if (backward_compat_compound && backward_compat)
156 vrna_fold_compound_free(backward_compat_compound);
157
158 backward_compat_compound = vc;
159 backward_compat = 1;
160 iindx = backward_compat_compound->iindx;
161
162 return vrna_pf(vc, structure);
163 }
164
165
166 PUBLIC vrna_ep_t *
stackProb(double cutoff)167 stackProb(double cutoff)
168 {
169 if (!(backward_compat_compound && backward_compat)) {
170 vrna_message_warning("stackProb: "
171 "run pf_fold() first!");
172 return NULL;
173 } else if (!backward_compat_compound->exp_matrices->probs) {
174 vrna_message_warning("stackProb: "
175 "probs == NULL!");
176 return NULL;
177 }
178
179 return vrna_stack_prob(backward_compat_compound, cutoff);
180 }
181
182
183 PUBLIC char *
centroid(int length,double * dist)184 centroid(int length,
185 double *dist)
186 {
187 if (pr == NULL) {
188 vrna_message_warning("centroid: "
189 "pr == NULL. You need to call pf_fold() before centroid()");
190 return NULL;
191 }
192
193 return vrna_centroid_from_probs(length, dist, pr);
194 }
195
196
197 PUBLIC double
mean_bp_dist(int length)198 mean_bp_dist(int length)
199 {
200 /* compute the mean base pair distance in the thermodynamic ensemble */
201 /* <d> = \sum_{a,b} p_a p_b d(S_a,S_b)
202 * this can be computed from the pair probs p_ij as
203 * <d> = \sum_{ij} p_{ij}(1-p_{ij}) */
204
205 int i, j, *my_iindx;
206 double d = 0;
207
208 if (pr == NULL) {
209 vrna_message_warning("mean_bp_dist: "
210 "pr == NULL. You need to call pf_fold() before mean_bp_dist()");
211 return d;
212 }
213
214 my_iindx = vrna_idx_row_wise(length);
215
216 for (i = 1; i <= length; i++)
217 for (j = i + TURN + 1; j <= length; j++)
218 d += pr[my_iindx[i] - j] * (1 - pr[my_iindx[i] - j]);
219
220 free(my_iindx);
221 return 2 * d;
222 }
223
224
225 /* get the free energy of a subsequence from the q[] array */
226 PUBLIC double
get_subseq_F(int i,int j)227 get_subseq_F(int i,
228 int j)
229 {
230 if (backward_compat_compound)
231 if (backward_compat_compound->exp_matrices)
232 if (backward_compat_compound->exp_matrices->q) {
233 int *my_iindx = backward_compat_compound->iindx;
234 vrna_exp_param_t *pf_params = backward_compat_compound->exp_params;
235 FLT_OR_DBL *q = backward_compat_compound->exp_matrices->q;
236 return (-log(q[my_iindx[i] - j]) - (j - i + 1) * log(pf_params->pf_scale)) * pf_params->kT /
237 1000.0;
238 }
239
240 vrna_message_warning("get_subseq_F: "
241 "call pf_fold() to fill q[] array before calling get_subseq_F()");
242
243 return 0.; /* we will never get to this point */
244 }
245
246
247 /*----------------------------------------------------------------------*/
248 PUBLIC double
expHairpinEnergy(int u,int type,short si1,short sj1,const char * string)249 expHairpinEnergy(int u,
250 int type,
251 short si1,
252 short sj1,
253 const char *string)
254 {
255 /* compute Boltzmann weight of a hairpin loop, multiply by scale[u+2] */
256
257 vrna_exp_param_t *pf_params = backward_compat_compound->exp_params;
258
259 double q, kT;
260
261 kT = pf_params->kT; /* kT in cal/mol */
262 if (u <= 30)
263 q = pf_params->exphairpin[u];
264 else
265 q = pf_params->exphairpin[30] * exp(-(pf_params->lxc * log(u / 30.)) * 10. / kT);
266
267 if ((tetra_loop) && (u == 4)) {
268 char tl[7] = {
269 0
270 }, *ts;
271 strncpy(tl, string, 6);
272 if ((ts = strstr(pf_params->Tetraloops, tl)))
273 return pf_params->exptetra[(ts - pf_params->Tetraloops) / 7];
274 }
275
276 if ((tetra_loop) && (u == 6)) {
277 char tl[9] = {
278 0
279 }, *ts;
280 strncpy(tl, string, 6);
281 if ((ts = strstr(pf_params->Hexaloops, tl)))
282 return pf_params->exphex[(ts - pf_params->Hexaloops) / 9];
283 }
284
285 if (u == 3) {
286 char tl[6] = {
287 0
288 }, *ts;
289 strncpy(tl, string, 5);
290 if ((ts = strstr(pf_params->Triloops, tl)))
291 return pf_params->exptri[(ts - pf_params->Triloops) / 6];
292
293 if (type > 2)
294 q *= pf_params->expTermAU;
295 } else {
296 /* no mismatches for tri-loops */
297 q *= pf_params->expmismatchH[type][si1][sj1];
298 }
299
300 return q;
301 }
302
303
304 PUBLIC double
expLoopEnergy(int u1,int u2,int type,int type2,short si1,short sj1,short sp1,short sq1)305 expLoopEnergy(int u1,
306 int u2,
307 int type,
308 int type2,
309 short si1,
310 short sj1,
311 short sp1,
312 short sq1)
313 {
314 /* compute Boltzmann weight of interior loop,
315 * multiply by scale[u1+u2+2] for scaling */
316 double z = 0;
317 int no_close = 0;
318 vrna_exp_param_t *pf_params = backward_compat_compound->exp_params;
319
320
321 if ((no_closingGU) && ((type2 == 3) || (type2 == 4) || (type == 2) || (type == 4)))
322 no_close = 1;
323
324 if ((u1 == 0) && (u2 == 0)) {
325 /* stack */
326 z = pf_params->expstack[type][type2];
327 } else if (no_close == 0) {
328 if ((u1 == 0) || (u2 == 0)) {
329 /* bulge */
330 int u;
331 u = (u1 == 0) ? u2 : u1;
332 z = pf_params->expbulge[u];
333 if (u2 + u1 == 1) {
334 z *= pf_params->expstack[type][type2];
335 } else {
336 if (type > 2)
337 z *= pf_params->expTermAU;
338
339 if (type2 > 2)
340 z *= pf_params->expTermAU;
341 }
342 } else {
343 /* interior loop */
344 if (u1 + u2 == 2) {
345 /* size 2 is special */
346 z = pf_params->expint11[type][type2][si1][sj1];
347 } else if ((u1 == 1) && (u2 == 2)) {
348 z = pf_params->expint21[type][type2][si1][sq1][sj1];
349 } else if ((u1 == 2) && (u2 == 1)) {
350 z = pf_params->expint21[type2][type][sq1][si1][sp1];
351 } else if ((u1 == 2) && (u2 == 2)) {
352 z = pf_params->expint22[type][type2][si1][sp1][sq1][sj1];
353 } else if (((u1 == 2) && (u2 == 3)) || ((u1 == 3) && (u2 == 2))) {
354 /*2-3 is special*/
355 z = pf_params->expinternal[5] *
356 pf_params->expmismatch23I[type][si1][sj1] *
357 pf_params->expmismatch23I[type2][sq1][sp1];
358 z *= pf_params->expninio[2][1];
359 } else if ((u1 == 1) || (u2 == 1)) {
360 /*1-n is special*/
361 z = pf_params->expinternal[u1 + u2] *
362 pf_params->expmismatch1nI[type][si1][sj1] *
363 pf_params->expmismatch1nI[type2][sq1][sp1];
364 z *= pf_params->expninio[2][abs(u1 - u2)];
365 } else {
366 z = pf_params->expinternal[u1 + u2] *
367 pf_params->expmismatchI[type][si1][sj1] *
368 pf_params->expmismatchI[type2][sq1][sp1];
369 z *= pf_params->expninio[2][abs(u1 - u2)];
370 }
371 }
372 }
373
374 return z;
375 }
376
377
378 PUBLIC void
init_pf_circ_fold(int length)379 init_pf_circ_fold(int length)
380 {
381 /* DO NOTHING */
382 }
383
384
385 PUBLIC void
init_pf_fold(int length)386 init_pf_fold(int length)
387 {
388 /* DO NOTHING */
389 }
390
391
392 /**
393 *** Allocate memory for all matrices and other stuff
394 **/
395 PUBLIC void
free_pf_arrays(void)396 free_pf_arrays(void)
397 {
398 if (backward_compat_compound && backward_compat) {
399 vrna_fold_compound_free(backward_compat_compound);
400 backward_compat_compound = NULL;
401 backward_compat = 0;
402 iindx = NULL;
403 }
404 }
405
406
407 PUBLIC FLT_OR_DBL *
export_bppm(void)408 export_bppm(void)
409 {
410 if (backward_compat_compound)
411 if (backward_compat_compound->exp_matrices)
412 if (backward_compat_compound->exp_matrices->probs)
413 return backward_compat_compound->exp_matrices->probs;
414
415 return NULL;
416 }
417
418
419 /*-------------------------------------------------------------------------*/
420 /* make arrays used for pf_fold available to other routines */
421 PUBLIC int
get_pf_arrays(short ** S_p,short ** S1_p,char ** ptype_p,FLT_OR_DBL ** qb_p,FLT_OR_DBL ** qm_p,FLT_OR_DBL ** q1k_p,FLT_OR_DBL ** qln_p)422 get_pf_arrays(short **S_p,
423 short **S1_p,
424 char **ptype_p,
425 FLT_OR_DBL **qb_p,
426 FLT_OR_DBL **qm_p,
427 FLT_OR_DBL **q1k_p,
428 FLT_OR_DBL **qln_p)
429 {
430 if (backward_compat_compound) {
431 if (backward_compat_compound->exp_matrices)
432 if (backward_compat_compound->exp_matrices->qb) {
433 *S_p = backward_compat_compound->sequence_encoding2;
434 *S1_p = backward_compat_compound->sequence_encoding;
435 *ptype_p = backward_compat_compound->ptype_pf_compat;
436 *qb_p = backward_compat_compound->exp_matrices->qb;
437 *qm_p = backward_compat_compound->exp_matrices->qm;
438 *q1k_p = backward_compat_compound->exp_matrices->q1k;
439 *qln_p = backward_compat_compound->exp_matrices->qln;
440 return 1;
441 }
442 }
443
444 return 0;
445 }
446
447
448 /*-----------------------------------------------------------------*/
449 PUBLIC float
pf_fold(const char * sequence,char * structure)450 pf_fold(const char *sequence,
451 char *structure)
452 {
453 return wrap_pf_fold(sequence, structure, NULL, do_backtrack, fold_constrained, 0);
454 }
455
456
457 PUBLIC float
pf_circ_fold(const char * sequence,char * structure)458 pf_circ_fold(const char *sequence,
459 char *structure)
460 {
461 return wrap_pf_fold(sequence, structure, NULL, do_backtrack, fold_constrained, 1);
462 }
463
464
465 PUBLIC float
pf_fold_par(const char * sequence,char * structure,vrna_exp_param_t * parameters,int calculate_bppm,int is_constrained,int is_circular)466 pf_fold_par(const char *sequence,
467 char *structure,
468 vrna_exp_param_t *parameters,
469 int calculate_bppm,
470 int is_constrained,
471 int is_circular)
472 {
473 return wrap_pf_fold(sequence, structure, parameters, calculate_bppm, is_constrained, is_circular);
474 }
475
476
477 PUBLIC char *
pbacktrack(char * seq)478 pbacktrack(char *seq)
479 {
480 int n = (int)strlen(seq);
481
482 return vrna_pbacktrack5(backward_compat_compound, n);
483 }
484
485
486 PUBLIC char *
pbacktrack5(char * seq,int length)487 pbacktrack5(char *seq,
488 int length)
489 {
490 /* the seq parameter must no differ to the one stored globally anyway, so we just ignore it */
491 return vrna_pbacktrack5(backward_compat_compound, length);
492 }
493
494
495 PUBLIC char *
pbacktrack_circ(char * seq)496 pbacktrack_circ(char *seq)
497 {
498 char *structure;
499 vrna_md_t *md;
500
501 structure = NULL;
502
503 if (backward_compat_compound) {
504 md = &(backward_compat_compound->exp_params->model_details);
505 if (md->circ && backward_compat_compound->exp_matrices->qm2)
506 structure = vrna_pbacktrack(backward_compat_compound);
507 }
508
509 return structure;
510 }
511
512
513 PUBLIC void
update_pf_params(int length)514 update_pf_params(int length)
515 {
516 if (backward_compat_compound && backward_compat) {
517 vrna_md_t md;
518 set_model_details(&md);
519 vrna_exp_params_reset(backward_compat_compound, &md);
520
521 /* compatibility with RNAup, may be removed sometime */
522 pf_scale = backward_compat_compound->exp_params->pf_scale;
523 }
524 }
525
526
527 PUBLIC void
update_pf_params_par(int length,vrna_exp_param_t * parameters)528 update_pf_params_par(int length,
529 vrna_exp_param_t *parameters)
530 {
531 if (backward_compat_compound && backward_compat) {
532 vrna_md_t md;
533 if (parameters) {
534 vrna_exp_params_subst(backward_compat_compound, parameters);
535 } else {
536 set_model_details(&md);
537 vrna_exp_params_reset(backward_compat_compound, &md);
538 }
539
540 /* compatibility with RNAup, may be removed sometime */
541 pf_scale = backward_compat_compound->exp_params->pf_scale;
542 }
543 }
544
545
546 PUBLIC char *
get_centroid_struct_gquad_pr(int length,double * dist)547 get_centroid_struct_gquad_pr(int length,
548 double *dist)
549 {
550 return vrna_centroid(backward_compat_compound, dist);
551 }
552
553
554 PUBLIC void
assign_plist_gquad_from_pr(vrna_ep_t ** pl,int length,double cut_off)555 assign_plist_gquad_from_pr(vrna_ep_t **pl,
556 int length, /* ignored */
557 double cut_off)
558 {
559 if (!backward_compat_compound)
560 *pl = NULL;
561 else if (!backward_compat_compound->exp_matrices->probs)
562 *pl = NULL;
563 else
564 *pl = vrna_plist_from_probs(backward_compat_compound, cut_off);
565 }
566
567
568 PUBLIC double
mean_bp_distance(int length)569 mean_bp_distance(int length)
570 {
571 if (backward_compat_compound)
572 if (backward_compat_compound->exp_matrices)
573 if (backward_compat_compound->exp_matrices->probs)
574 return vrna_mean_bp_distance(backward_compat_compound);
575
576 vrna_message_warning("mean_bp_distance: "
577 "you need to call vrna_pf_fold first");
578
579 return 0.; /* we will never get to this point */
580 }
581
582
583 PUBLIC double
mean_bp_distance_pr(int length,FLT_OR_DBL * p)584 mean_bp_distance_pr(int length,
585 FLT_OR_DBL *p)
586 {
587 double d = 0;
588 int *index = vrna_idx_row_wise((unsigned int)length);
589
590 if (p == NULL) {
591 vrna_message_warning("mean_bp_distance_pr: "
592 "p == NULL. You need to supply a valid probability matrix for mean_bp_distance_pr()");
593 return d;
594 }
595
596 d = wrap_mean_bp_distance(p, length, index, TURN);
597
598 free(index);
599 return d;
600 }
601
602
603 #endif
604