1 /*
2 * partiton function for RNA secondary structures
3 *
4 * Ivo L Hofacker
5 * Stephan Bernhart
6 * Ronny Lorenz
7 * ViennaRNA package
8 */
9
10 #ifdef HAVE_CONFIG_H
11 #include "config.h"
12 #endif
13
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include <math.h>
18 #include <float.h> /* #defines FLT_MAX ... */
19 #include <limits.h>
20
21 #include "ViennaRNA/utils/basic.h"
22 #include "ViennaRNA/utils/structures.h"
23 #include "ViennaRNA/params/default.h"
24 #include "ViennaRNA/fold_vars.h"
25 #include "ViennaRNA/plotting/probabilities.h"
26 #include "ViennaRNA/params/basic.h"
27 #include "ViennaRNA/loops/all.h"
28 #include "ViennaRNA/mfe.h"
29 #include "ViennaRNA/part_func.h"
30 #include "ViennaRNA/part_func_co.h"
31
32 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
33 #ifdef _OPENMP
34 #include <omp.h>
35 #endif
36 #endif
37
38 /*
39 #################################
40 # GLOBAL VARIABLES #
41 #################################
42 */
43
44 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
45
46 int mirnatog = 0;
47 double F_monomer[2] = {
48 0, 0
49 }; /* free energies of the two monomers */
50
51 #endif
52
53 /*
54 #################################
55 # PRIVATE VARIABLES #
56 #################################
57 */
58
59 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
60
61 /* some backward compatibility stuff */
62 PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
63 PRIVATE int backward_compat = 0;
64
65 #ifdef _OPENMP
66
67 #pragma omp threadprivate(backward_compat_compound, backward_compat)
68
69 #endif
70
71 #endif
72
73 /*
74 #################################
75 # PRIVATE FUNCTION DECLARATIONS #
76 #################################
77 */
78
79 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
80
81 PRIVATE vrna_dimer_pf_t
82 wrap_co_pf_fold(char *sequence,
83 char *structure,
84 vrna_exp_param_t *parameters,
85 int calculate_bppm,
86 int is_constrained);
87
88
89 #endif
90
91 /*
92 #################################
93 # BEGIN OF FUNCTION DEFINITIONS #
94 #################################
95 */
96 PUBLIC vrna_dimer_pf_t
vrna_pf_co_fold(const char * seq,char * structure,vrna_ep_t ** pl)97 vrna_pf_co_fold(const char *seq,
98 char *structure,
99 vrna_ep_t **pl)
100 {
101 double mfe;
102 vrna_dimer_pf_t X;
103 vrna_fold_compound_t *vc;
104 vrna_md_t md;
105
106 vrna_md_set_default(&md);
107
108 /* no need to backtrack MFE structure */
109 md.backtrack = 0;
110 if (pl)
111 md.compute_bpp = 1;
112 else
113 md.compute_bpp = 0;
114
115 vc = vrna_fold_compound(seq, &md, VRNA_OPTION_DEFAULT);
116
117 mfe = (double)vrna_mfe_dimer(vc, NULL);
118 vrna_exp_params_rescale(vc, &mfe);
119
120 X = vrna_pf_dimer(vc, structure);
121
122 if (pl)
123 *pl = vrna_plist_from_probs(vc, /*cut_off:*/ 1e-6);
124
125 vrna_fold_compound_free(vc);
126
127 return X;
128 }
129
130
131 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
132
133 /*
134 *****************************************
135 * BEGIN backward compatibility wrappers *
136 *****************************************
137 */
138 PRIVATE vrna_dimer_pf_t
wrap_co_pf_fold(char * sequence,char * structure,vrna_exp_param_t * parameters,int calculate_bppm,int is_constrained)139 wrap_co_pf_fold(char *sequence,
140 char *structure,
141 vrna_exp_param_t *parameters,
142 int calculate_bppm,
143 int is_constrained)
144 {
145 int length;
146 char *seq;
147 vrna_fold_compound_t *vc;
148 vrna_md_t md;
149
150 vc = NULL;
151 length = strlen(sequence);
152
153 seq = (char *)vrna_alloc(sizeof(char) * (length + 2));
154 if (cut_point > -1) {
155 int i;
156 for (i = 0; i < cut_point - 1; i++)
157 seq[i] = sequence[i];
158 seq[i] = '&';
159 for (; i < (int)length; i++)
160 seq[i + 1] = sequence[i];
161 } else {
162 /* this ensures the allocation of all cofold matrices via vrna_fold_compound_t */
163 free(seq);
164 seq = strdup(sequence);
165 }
166
167 /*
168 * if present, extract model details from provided parameters variable,
169 * to properly initialize the fold compound. Otherwise use default
170 * settings taken from deprecated global variables
171 */
172 if (parameters)
173 vrna_md_copy(&md, &(parameters->model_details));
174 else
175 set_model_details(&md);
176
177 /* set min_loop_size and backtracing options */
178 md.compute_bpp = calculate_bppm;
179 md.min_loop_size = 0;
180
181 vc = vrna_fold_compound(seq, &md, VRNA_OPTION_DEFAULT);
182
183 /*
184 * if present, attach a copy of the parameters structure instead of the
185 * default parameters but take care of re-setting it to (initialized)
186 * model details
187 */
188 free(vc->exp_params);
189 if (parameters) {
190 vrna_md_copy(&(parameters->model_details), &(vc->params->model_details));
191 vc->exp_params = vrna_exp_params_copy(parameters);
192 } else {
193 vc->exp_params = vrna_exp_params(&(vc->params->model_details));
194 }
195
196 /* propagate global pf_scale into vc->exp_params */
197 vc->exp_params->pf_scale = pf_scale;
198
199 if (is_constrained && structure) {
200 unsigned int constraint_options = 0;
201 constraint_options |= VRNA_CONSTRAINT_DB
202 | VRNA_CONSTRAINT_DB_PIPE
203 | VRNA_CONSTRAINT_DB_DOT
204 | VRNA_CONSTRAINT_DB_X
205 | VRNA_CONSTRAINT_DB_ANG_BRACK
206 | VRNA_CONSTRAINT_DB_RND_BRACK;
207
208 vrna_constraints_add(vc, (const char *)structure, constraint_options);
209 }
210
211 if (backward_compat_compound)
212 vrna_fold_compound_free(backward_compat_compound);
213
214 backward_compat_compound = vc;
215 backward_compat = 1;
216 iindx = backward_compat_compound->iindx;
217
218 free(seq);
219 return vrna_pf_dimer(vc, structure);
220 }
221
222
223 /*
224 *****************************************
225 * END backward compatibility wrappers *
226 *****************************************
227 */
228
229
230 /*###########################################*/
231 /*# deprecated functions below #*/
232 /*###########################################*/
233
234 PUBLIC vrna_dimer_pf_t
co_pf_fold(char * sequence,char * structure)235 co_pf_fold(char *sequence,
236 char *structure)
237 {
238 return wrap_co_pf_fold(sequence, structure, NULL, do_backtrack, fold_constrained);
239 }
240
241
242 PUBLIC vrna_dimer_pf_t
co_pf_fold_par(char * sequence,char * structure,vrna_exp_param_t * parameters,int calculate_bppm,int is_constrained)243 co_pf_fold_par(char *sequence,
244 char *structure,
245 vrna_exp_param_t *parameters,
246 int calculate_bppm,
247 int is_constrained)
248 {
249 return wrap_co_pf_fold(sequence, structure, parameters, calculate_bppm, is_constrained);
250 }
251
252
253 PUBLIC vrna_ep_t *
get_plist(vrna_ep_t * pl,int length,double cut_off)254 get_plist(vrna_ep_t *pl,
255 int length,
256 double cut_off)
257 {
258 int i, j, n, count, *my_iindx;
259
260 my_iindx = backward_compat_compound->iindx;
261 /* get pair probibilities out of pr array */
262 count = 0;
263 n = 2;
264 for (i = 1; i < length; i++) {
265 for (j = i + 1; j <= length; j++) {
266 if (pr[my_iindx[i] - j] < cut_off)
267 continue;
268
269 if (count == n * length - 1) {
270 n *= 2;
271 pl = (vrna_ep_t *)vrna_realloc(pl, n * length * sizeof(vrna_ep_t));
272 }
273
274 pl[count].i = i;
275 pl[count].j = j;
276 pl[count++].p = pr[my_iindx[i] - j];
277 /* printf("gpl: %2d %2d %.9f\n",i,j,pr[my_iindx[i]-j]); */
278 }
279 }
280 pl[count].i = 0;
281 pl[count].j = 0; /* ->?? */
282 pl[count++].p = 0.;
283 pl = (vrna_ep_t *)vrna_realloc(pl, (count) * sizeof(vrna_ep_t));
284 return pl;
285 }
286
287
288 PUBLIC void
compute_probabilities(double FAB,double FA,double FB,vrna_ep_t * prAB,vrna_ep_t * prA,vrna_ep_t * prB,int Alength)289 compute_probabilities(double FAB,
290 double FA,
291 double FB,
292 vrna_ep_t *prAB,
293 vrna_ep_t *prA,
294 vrna_ep_t *prB,
295 int Alength)
296 {
297 if (backward_compat_compound && backward_compat) {
298 vrna_pf_dimer_probs(FAB,
299 FA,
300 FB,
301 prAB,
302 (const vrna_ep_t *)prA,
303 (const vrna_ep_t *)prB,
304 Alength,
305 (const vrna_exp_param_t *)backward_compat_compound->exp_params);
306 }
307 }
308
309
310 PUBLIC vrna_dimer_conc_t *
get_concentrations(double FcAB,double FcAA,double FcBB,double FEA,double FEB,double * startconc)311 get_concentrations(double FcAB,
312 double FcAA,
313 double FcBB,
314 double FEA,
315 double FEB,
316 double *startconc)
317 {
318 return vrna_pf_dimer_concentrations(FcAB,
319 FcAA,
320 FcBB,
321 FEA,
322 FEB,
323 (const double *)startconc,
324 (const vrna_exp_param_t *)backward_compat_compound->exp_params);
325 }
326
327
328 PUBLIC void
init_co_pf_fold(int length)329 init_co_pf_fold(int length)
330 {
331 /* DO NOTHING */
332 }
333
334
335 PUBLIC void
free_co_pf_arrays(void)336 free_co_pf_arrays(void)
337 {
338 if (backward_compat_compound && backward_compat) {
339 vrna_fold_compound_free(backward_compat_compound);
340 backward_compat_compound = NULL;
341 backward_compat = 0;
342 }
343 }
344
345
346 PUBLIC FLT_OR_DBL *
export_co_bppm(void)347 export_co_bppm(void)
348 {
349 if (backward_compat_compound)
350 return backward_compat_compound->exp_matrices->probs;
351 else
352 return NULL;
353 }
354
355
356 PUBLIC void
update_co_pf_params(int length)357 update_co_pf_params(int length)
358 {
359 if (backward_compat_compound && backward_compat) {
360 vrna_md_t md;
361 set_model_details(&md);
362 vrna_exp_params_reset(backward_compat_compound, &md);
363
364 /* compatibility with RNAup, may be removed sometime */
365 pf_scale = backward_compat_compound->exp_params->pf_scale;
366 }
367 }
368
369
370 PUBLIC void
update_co_pf_params_par(int length,vrna_exp_param_t * parameters)371 update_co_pf_params_par(int length,
372 vrna_exp_param_t *parameters)
373 {
374 if (backward_compat_compound && backward_compat) {
375 vrna_md_t md;
376 if (parameters) {
377 vrna_exp_params_subst(backward_compat_compound, parameters);
378 } else {
379 set_model_details(&md);
380 vrna_exp_params_reset(backward_compat_compound, &md);
381 }
382
383 /* compatibility with RNAup, may be removed sometime */
384 pf_scale = backward_compat_compound->exp_params->pf_scale;
385 }
386 }
387
388
389 #endif
390