1 /*
2 *
3 * c Ivo Hofacker
4 *
5 * Vienna RNA package
6 */
7
8 #ifdef HAVE_CONFIG_H
9 #include "config.h"
10 #endif
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <math.h>
15 #include <string.h>
16 #include "ViennaRNA/params/default.h"
17 #include "ViennaRNA/fold_vars.h"
18 #include "ViennaRNA/utils/basic.h"
19 #include "ViennaRNA/params/io.h"
20 #include "ViennaRNA/params/basic.h"
21
22 /**
23 *** \file ViennaRNA/params/basic.c
24 *** <P>
25 *** This file provides functions that return temperature scaled energy parameters and
26 *** Boltzmann weights packed in datastructures
27 *** </P>
28 ***/
29
30 /*------------------------------------------------------------------------*/
31 #define SCALE 10
32 /**
33 *** dangling ends should never be destabilizing, i.e. expdangle>=1<BR>
34 *** specific heat needs smooth function (2nd derivative)<BR>
35 *** we use a*(sin(x+b)+1)^2, with a=2/(3*sqrt(3)), b=Pi/6-sqrt(3)/2,
36 *** in the interval b<x<sqrt(3)/2
37 */
38 #define CLIP_NEGATIVE(X) ((X) < 0 ? 0 : (X))
39 #define SMOOTH(X) ((!pf_smooth) ? CLIP_NEGATIVE(X) : \
40 ((X) / SCALE < -1.2283697) ? 0 : \
41 ((X) / SCALE > 0.8660254) ? (X) : \
42 SCALE *0.38490018 \
43 * (sin((X) / SCALE - 0.34242663) + 1) \
44 * (sin((X) / SCALE - 0.34242663) + 1) \
45 )
46
47 /* #define SMOOTH(X) ((X)<0 ? 0 : (X)) */
48
49 /*
50 * If the global use_mfelike_energies flag is set, truncate doubles to int
51 * values and cast back to double. This makes the energy parameters of the
52 * partition (folding get_scaled_exp_params()) compatible with the mfe folding
53 * parameters (get_scaled_exp_params()), e.g. for explicit partition function
54 * computations.
55 */
56 #define TRUNC_MAYBE(X) ((!pf_smooth) ? (double)((int)(X)) : (X))
57
58
59 /* Rescale Free energy contribution according to deviation of temperature from measurement conditions */
60 #define RESCALE_dG(dG, dH, dT) ((dH) - ((dH) - (dG)) * dT)
61
62 /*
63 * Rescale Free energy contribution according to deviation of temperature from measurement conditions
64 * and convert it to Boltzmann Factor for specific kT
65 */
66 #define RESCALE_BF(dG, dH, dT, kT) ( \
67 exp( \
68 -TRUNC_MAYBE((double)RESCALE_dG((dG), (dH), (dT))) \
69 * 10. \
70 / kT \
71 ) \
72 )
73
74 #define RESCALE_BF_SMOOTH(dG, dH, dT, kT) ( \
75 exp( \
76 SMOOTH( \
77 -TRUNC_MAYBE((double)RESCALE_dG((dG), (dH), (dT))) \
78 ) \
79 * 10. \
80 / kT \
81 ) \
82 )
83
84 /*
85 #################################
86 # PRIVATE VARIABLES #
87 #################################
88 */
89 PRIVATE vrna_param_t p;
90 PRIVATE int id = -1;
91 /* variables for partition function */
92 PRIVATE vrna_exp_param_t pf;
93 PRIVATE int pf_id = -1;
94
95 #ifdef _OPENMP
96 #pragma omp threadprivate(id, pf_id)
97 #endif
98
99 /*
100 #################################
101 # PRIVATE FUNCTION DECLARATIONS #
102 #################################
103 */
104
105 PRIVATE vrna_param_t *
106 get_scaled_params(vrna_md_t *md);
107
108
109 PRIVATE vrna_exp_param_t *
110 get_scaled_exp_params(vrna_md_t *md,
111 double pfs);
112
113
114 PRIVATE vrna_exp_param_t *
115 get_exp_params_ali(vrna_md_t *md,
116 unsigned int n_seq,
117 double pfs);
118
119
120 PRIVATE void
121 rescale_params(vrna_fold_compound_t *vc);
122
123
124 /*
125 #################################
126 # BEGIN OF FUNCTION DEFINITIONS #
127 #################################
128 */
129 PUBLIC vrna_param_t *
vrna_params(vrna_md_t * md)130 vrna_params(vrna_md_t *md)
131 {
132 if (md) {
133 return get_scaled_params(md);
134 } else {
135 vrna_md_t md;
136 vrna_md_set_default(&md);
137 return get_scaled_params(&md);
138 }
139 }
140
141
142 PUBLIC vrna_exp_param_t *
vrna_exp_params(vrna_md_t * md)143 vrna_exp_params(vrna_md_t *md)
144 {
145 if (md) {
146 return get_scaled_exp_params(md, -1.);
147 } else {
148 vrna_md_t md;
149 vrna_md_set_default(&md);
150 return get_scaled_exp_params(&md, -1.);
151 }
152 }
153
154
155 PUBLIC vrna_exp_param_t *
vrna_exp_params_comparative(unsigned int n_seq,vrna_md_t * md)156 vrna_exp_params_comparative(unsigned int n_seq,
157 vrna_md_t *md)
158 {
159 if (md) {
160 return get_exp_params_ali(md, n_seq, -1.);
161 } else {
162 vrna_md_t md;
163 vrna_md_set_default(&md);
164 return get_exp_params_ali(&md, n_seq, -1.);
165 }
166 }
167
168
169 PUBLIC vrna_param_t *
vrna_params_copy(vrna_param_t * par)170 vrna_params_copy(vrna_param_t *par)
171 {
172 vrna_param_t *copy = NULL;
173
174 if (par) {
175 copy = (vrna_param_t *)vrna_alloc(sizeof(vrna_param_t));
176 memcpy(copy, par, sizeof(vrna_param_t));
177 }
178
179 return copy;
180 }
181
182
183 PUBLIC vrna_exp_param_t *
vrna_exp_params_copy(vrna_exp_param_t * par)184 vrna_exp_params_copy(vrna_exp_param_t *par)
185 {
186 vrna_exp_param_t *copy = NULL;
187
188 if (par) {
189 copy = (vrna_exp_param_t *)vrna_alloc(sizeof(vrna_exp_param_t));
190 memcpy(copy, par, sizeof(vrna_exp_param_t));
191 }
192
193 return copy;
194 }
195
196
197 PUBLIC void
vrna_params_subst(vrna_fold_compound_t * vc,vrna_param_t * parameters)198 vrna_params_subst(vrna_fold_compound_t *vc,
199 vrna_param_t *parameters)
200 {
201 if (vc) {
202 if (vc->params)
203 free(vc->params);
204
205 if (parameters) {
206 vc->params = vrna_params_copy(parameters);
207 } else {
208 switch (vc->type) {
209 case VRNA_FC_TYPE_SINGLE: /* fall through */
210
211 case VRNA_FC_TYPE_COMPARATIVE:
212 vc->params = vrna_params(NULL);
213 break;
214
215 default:
216 break;
217 }
218 }
219 }
220 }
221
222
223 PUBLIC void
vrna_params_reset(vrna_fold_compound_t * vc,vrna_md_t * md_p)224 vrna_params_reset(vrna_fold_compound_t *vc,
225 vrna_md_t *md_p)
226 {
227 if (vc) {
228 switch (vc->type) {
229 case VRNA_FC_TYPE_SINGLE: /* fall through */
230
231 case VRNA_FC_TYPE_COMPARATIVE:
232 if (vc->params)
233 free(vc->params);
234
235 vc->params = vrna_params(md_p);
236
237 if (vc->exp_params) {
238 free(vc->exp_params);
239
240 vc->exp_params = vrna_exp_params(md_p);
241 }
242
243 break;
244
245 default:
246 break;
247 }
248 }
249 }
250
251
252 PUBLIC void
vrna_exp_params_reset(vrna_fold_compound_t * vc,vrna_md_t * md_p)253 vrna_exp_params_reset(vrna_fold_compound_t *vc,
254 vrna_md_t *md_p)
255 {
256 if (vc) {
257 switch (vc->type) {
258 case VRNA_FC_TYPE_SINGLE: /* fall through */
259
260 case VRNA_FC_TYPE_COMPARATIVE:
261 if (vc->exp_params)
262 free(vc->exp_params);
263
264 vc->exp_params = vrna_exp_params(md_p);
265 break;
266
267 default:
268 break;
269 }
270 }
271 }
272
273
274 PUBLIC void
vrna_exp_params_subst(vrna_fold_compound_t * vc,vrna_exp_param_t * params)275 vrna_exp_params_subst(vrna_fold_compound_t *vc,
276 vrna_exp_param_t *params)
277 {
278 if (vc) {
279 if (vc->exp_params)
280 free(vc->exp_params);
281
282 if (params) {
283 vc->exp_params = vrna_exp_params_copy(params);
284 } else {
285 switch (vc->type) {
286 case VRNA_FC_TYPE_SINGLE:
287 vc->exp_params = vrna_exp_params(NULL);
288 if (vc->strands > 1)
289 vc->exp_params->model_details.min_loop_size = 0;
290
291 break;
292
293 case VRNA_FC_TYPE_COMPARATIVE:
294 vc->exp_params = vrna_exp_params_comparative(vc->n_seq, NULL);
295 break;
296
297 default:
298 break;
299 }
300 }
301
302 /* fill additional helper arrays for scaling etc. */
303 vrna_exp_params_rescale(vc, NULL);
304 }
305 }
306
307
308 PUBLIC void
vrna_exp_params_rescale(vrna_fold_compound_t * vc,double * mfe)309 vrna_exp_params_rescale(vrna_fold_compound_t *vc,
310 double *mfe)
311 {
312 vrna_exp_param_t *pf;
313 double e_per_nt, kT;
314 vrna_md_t *md;
315
316 if (vc) {
317 if (!vc->exp_params) {
318 switch (vc->type) {
319 case VRNA_FC_TYPE_SINGLE:
320 vc->exp_params = vrna_exp_params(&(vc->params->model_details));
321 break;
322 case VRNA_FC_TYPE_COMPARATIVE:
323 vc->exp_params = vrna_exp_params_comparative(vc->n_seq, &(vc->params->model_details));
324 break;
325 }
326 } else if (memcmp(&(vc->params->model_details),
327 &(vc->exp_params->model_details),
328 sizeof(vrna_md_t)) != 0) {
329 /* make sure that model details are matching */
330 (void)vrna_md_copy(&(vc->exp_params->model_details), &(vc->params->model_details));
331 /* we probably need some mechanism to check whether DP matrices still match the new model settings! */
332 }
333
334 pf = vc->exp_params;
335 if (pf) {
336 kT = pf->kT;
337 md = &(pf->model_details);
338
339 if (vc->type == VRNA_FC_TYPE_COMPARATIVE)
340 kT /= vc->n_seq;
341
342 /* re-compute scaling factor if necessary */
343 if ((mfe) || (pf->pf_scale < 1.)) {
344 if (mfe) /* use largest known Boltzmann factor for scaling */
345 e_per_nt = *mfe * 1000. / vc->length;
346 else /* use mean energy for random sequences: 184.3*length cal for scaling */
347 e_per_nt = -185 + (pf->temperature - 37.) * 7.27;
348
349 /* apply user-defined scaling factor to allow scaling for unusually stable/unstable structure enembles */
350 pf->pf_scale = exp(-(md->sfact * e_per_nt) / kT);
351 }
352
353 if (pf->pf_scale < 1.)
354 pf->pf_scale = 1.;
355
356 rescale_params(vc);
357 }
358 }
359 }
360
361
362 PUBLIC void
vrna_params_prepare(vrna_fold_compound_t * fc,unsigned int options)363 vrna_params_prepare(vrna_fold_compound_t *fc,
364 unsigned int options)
365 {
366 if (fc) {
367 vrna_md_t *md_p;
368
369 /*
370 * every vrna_fold_compound_t must have a vrna_paramt_t structure attached
371 * to it that holds the current model details. So we just use this here as
372 * the reference model
373 */
374 md_p = &(fc->params->model_details);
375
376 if (options & VRNA_OPTION_PF) {
377 /* remove previous parameters if present and they differ from reference model */
378 if (fc->exp_params) {
379 if (memcmp(md_p, &(fc->exp_params->model_details), sizeof(vrna_md_t)) != 0) {
380 free(fc->exp_params);
381 fc->exp_params = NULL;
382 }
383 }
384
385 if (!fc->exp_params)
386 fc->exp_params = (fc->type == VRNA_FC_TYPE_SINGLE) ? \
387 vrna_exp_params(md_p) : \
388 vrna_exp_params_comparative(fc->n_seq, md_p);
389 }
390 }
391 }
392
393
394 /*
395 #####################################
396 # BEGIN OF STATIC HELPER FUNCTIONS #
397 #####################################
398 */
399 PRIVATE vrna_param_t *
get_scaled_params(vrna_md_t * md)400 get_scaled_params(vrna_md_t *md)
401 {
402 unsigned int i, j, k, l;
403 double tempf;
404 vrna_param_t *params;
405
406 params = (vrna_param_t *)vrna_alloc(sizeof(vrna_param_t));
407
408 memset(params->param_file, '\0', 256);
409 if (last_parameter_file() != NULL)
410 strncpy(params->param_file, last_parameter_file(), 255);
411
412 params->model_details = *md; /* copy over the model details */
413 params->temperature = md->temperature;
414 tempf = ((params->temperature + K0) / Tmeasure);
415
416 params->ninio[2] = RESCALE_dG(ninio37, niniodH, tempf);
417 params->lxc = lxc37 * tempf;
418 params->TripleC = RESCALE_dG(TripleC37, TripleCdH, tempf);
419 params->MultipleCA = RESCALE_dG(MultipleCA37, MultipleCAdH, tempf);
420 params->MultipleCB = RESCALE_dG(MultipleCB37, MultipleCBdH, tempf);
421 params->TerminalAU = RESCALE_dG(TerminalAU37, TerminalAUdH, tempf);
422 params->DuplexInit = RESCALE_dG(DuplexInit37, DuplexInitdH, tempf);
423 params->MLbase = RESCALE_dG(ML_BASE37, ML_BASEdH, tempf);
424 params->MLclosing = RESCALE_dG(ML_closing37, ML_closingdH, tempf);
425 params->gquadLayerMismatch = RESCALE_dG(GQuadLayerMismatch37, GQuadLayerMismatchH, tempf);
426 params->gquadLayerMismatchMax = GQuadLayerMismatchMax;
427
428 for (i = VRNA_GQUAD_MIN_STACK_SIZE; i <= VRNA_GQUAD_MAX_STACK_SIZE; i++)
429 for (j = 3 * VRNA_GQUAD_MIN_LINKER_LENGTH; j <= 3 * VRNA_GQUAD_MAX_LINKER_LENGTH; j++) {
430 double GQuadAlpha_T = RESCALE_dG(GQuadAlpha37, GQuadAlphadH, tempf);
431 double GQuadBeta_T = RESCALE_dG(GQuadBeta37, GQuadBetadH, tempf);
432 params->gquad[i][j] = (int)GQuadAlpha_T * (i - 1) + (int)(((double)GQuadBeta_T) * log(j - 2));
433 }
434
435 for (i = 0; i < 31; i++)
436 params->hairpin[i] = RESCALE_dG(hairpin37[i], hairpindH[i], tempf);
437
438 for (i = 0; i <= MIN2(30, MAXLOOP); i++) {
439 params->bulge[i] = RESCALE_dG(bulge37[i], bulgedH[i], tempf);
440 params->internal_loop[i] = RESCALE_dG(internal_loop37[i], internal_loopdH[i], tempf);
441 }
442
443 for (; i <= MAXLOOP; i++) {
444 params->bulge[i] = params->bulge[30] +
445 (int)(params->lxc * log((double)(i) / 30.));
446 params->internal_loop[i] = params->internal_loop[30] +
447 (int)(params->lxc * log((double)(i) / 30.));
448 }
449
450 for (i = 0; (i * 7) < strlen(Tetraloops); i++)
451 params->Tetraloop_E[i] = RESCALE_dG(Tetraloop37[i], TetraloopdH[i], tempf);
452
453 for (i = 0; (i * 5) < strlen(Triloops); i++)
454 params->Triloop_E[i] = RESCALE_dG(Triloop37[i], TriloopdH[i], tempf);
455
456 for (i = 0; (i * 9) < strlen(Hexaloops); i++)
457 params->Hexaloop_E[i] = RESCALE_dG(Hexaloop37[i], HexaloopdH[i], tempf);
458
459 for (i = 0; i <= NBPAIRS; i++)
460 params->MLintern[i] = RESCALE_dG(ML_intern37, ML_interndH, tempf);
461
462 /* stacks G(T) = H - [H - G(T0)]*T/T0 */
463 for (i = 0; i <= NBPAIRS; i++)
464 for (j = 0; j <= NBPAIRS; j++)
465 params->stack[i][j] = RESCALE_dG(stack37[i][j],
466 stackdH[i][j],
467 tempf);
468
469 /* mismatches */
470 for (i = 0; i <= NBPAIRS; i++)
471 for (j = 0; j < 5; j++)
472 for (k = 0; k < 5; k++) {
473 int mm;
474 params->mismatchI[i][j][k] = RESCALE_dG(mismatchI37[i][j][k],
475 mismatchIdH[i][j][k],
476 tempf);
477 params->mismatchH[i][j][k] = RESCALE_dG(mismatchH37[i][j][k],
478 mismatchHdH[i][j][k],
479 tempf);
480 params->mismatch1nI[i][j][k] = RESCALE_dG(mismatch1nI37[i][j][k],
481 mismatch1nIdH[i][j][k],
482 tempf);
483 params->mismatch23I[i][j][k] = RESCALE_dG(mismatch23I37[i][j][k],
484 mismatch23IdH[i][j][k],
485 tempf);
486 if (md->dangles) {
487 mm = RESCALE_dG(mismatchM37[i][j][k],
488 mismatchMdH[i][j][k],
489 tempf);
490 params->mismatchM[i][j][k] = (mm > 0) ? 0 : mm;
491 mm = RESCALE_dG(mismatchExt37[i][j][k],
492 mismatchExtdH[i][j][k],
493 tempf);
494 params->mismatchExt[i][j][k] = (mm > 0) ? 0 : mm;
495 } else {
496 params->mismatchM[i][j][k] = params->mismatchExt[i][j][k] = 0;
497 }
498 }
499
500 /* dangles */
501 for (i = 0; i <= NBPAIRS; i++)
502 for (j = 0; j < 5; j++) {
503 int dd;
504 dd = RESCALE_dG(dangle5_37[i][j],
505 dangle5_dH[i][j],
506 tempf);
507 params->dangle5[i][j] = (dd > 0) ? 0 : dd; /* must be <= 0 */
508 dd = RESCALE_dG(dangle3_37[i][j],
509 dangle3_dH[i][j],
510 tempf);
511 params->dangle3[i][j] = (dd > 0) ? 0 : dd; /* must be <= 0 */
512 }
513
514 /* interior 1x1 loops */
515 for (i = 0; i <= NBPAIRS; i++)
516 for (j = 0; j <= NBPAIRS; j++)
517 for (k = 0; k < 5; k++)
518 for (l = 0; l < 5; l++)
519 params->int11[i][j][k][l] = RESCALE_dG(int11_37[i][j][k][l],
520 int11_dH[i][j][k][l],
521 tempf);
522
523 /* interior 2x1 loops */
524 for (i = 0; i <= NBPAIRS; i++)
525 for (j = 0; j <= NBPAIRS; j++)
526 for (k = 0; k < 5; k++)
527 for (l = 0; l < 5; l++) {
528 int m;
529 for (m = 0; m < 5; m++)
530 params->int21[i][j][k][l][m] = RESCALE_dG(int21_37[i][j][k][l][m],
531 int21_dH[i][j][k][l][m],
532 tempf);
533 }
534
535 /* interior 2x2 loops */
536 for (i = 0; i <= NBPAIRS; i++)
537 for (j = 0; j <= NBPAIRS; j++)
538 for (k = 0; k < 5; k++)
539 for (l = 0; l < 5; l++) {
540 int m, n;
541 for (m = 0; m < 5; m++)
542 for (n = 0; n < 5; n++)
543 params->int22[i][j][k][l][m][n] = RESCALE_dG(int22_37[i][j][k][l][m][n],
544 int22_dH[i][j][k][l][m][n],
545 tempf);
546 }
547
548 strncpy(params->Tetraloops, Tetraloops, 281);
549 strncpy(params->Triloops, Triloops, 241);
550 strncpy(params->Hexaloops, Hexaloops, 361);
551
552 params->id = ++id;
553 return params;
554 }
555
556
557 PRIVATE vrna_exp_param_t *
get_scaled_exp_params(vrna_md_t * md,double pfs)558 get_scaled_exp_params(vrna_md_t *md,
559 double pfs)
560 {
561 unsigned int i, j, k, l;
562 int pf_smooth;
563 double kT, TT;
564 double GT;
565 vrna_exp_param_t *pf;
566
567 pf = (vrna_exp_param_t *)vrna_alloc(sizeof(vrna_exp_param_t));
568
569 memset(pf->param_file, '\0', 256);
570 if (last_parameter_file() != NULL)
571 strncpy(pf->param_file, last_parameter_file(), 255);
572
573 pf->model_details = *md;
574 pf->temperature = md->temperature;
575 pf->alpha = md->betaScale;
576 pf->kT = kT = md->betaScale * (md->temperature + K0) * GASCONST; /* kT in cal/mol */
577 pf->pf_scale = pfs;
578 pf_smooth = md->pf_smooth;
579 TT = (md->temperature + K0) / (Tmeasure);
580
581 pf->lxc = lxc37 * TT;
582 pf->expDuplexInit = RESCALE_BF(DuplexInit37, DuplexInitdH, TT, kT);
583 pf->expTermAU = RESCALE_BF(TerminalAU37, TerminalAUdH, TT, kT);
584 pf->expMLbase = RESCALE_BF(ML_BASE37, ML_BASEdH, TT, kT);
585 pf->expMLclosing = RESCALE_BF(ML_closing37, ML_closingdH, TT, kT);
586 pf->expgquadLayerMismatch = RESCALE_BF(GQuadLayerMismatch37, GQuadLayerMismatchH, TT, kT);
587 pf->gquadLayerMismatchMax = GQuadLayerMismatchMax;
588
589 for (i = VRNA_GQUAD_MIN_STACK_SIZE; i <= VRNA_GQUAD_MAX_STACK_SIZE; i++)
590 for (j = 3 * VRNA_GQUAD_MIN_LINKER_LENGTH; j <= 3 * VRNA_GQUAD_MAX_LINKER_LENGTH; j++) {
591 double GQuadAlpha_T = RESCALE_dG(GQuadAlpha37, GQuadAlphadH, TT);
592 double GQuadBeta_T = RESCALE_dG(GQuadBeta37, GQuadBetadH, TT);
593 GT = ((double)GQuadAlpha_T) * ((double)(i - 1)) + ((double)GQuadBeta_T) *
594 log(((double)j) - 2.);
595 pf->expgquad[i][j] = exp(-TRUNC_MAYBE(GT) * 10. / kT);
596 }
597
598 /* loop energies: hairpins, bulges, interior, mulit-loops */
599 for (i = 0; i < 31; i++)
600 pf->exphairpin[i] = RESCALE_BF(hairpin37[i], hairpindH[i], TT, kT);
601
602 for (i = 0; i <= MIN2(30, MAXLOOP); i++) {
603 pf->expbulge[i] = RESCALE_BF(bulge37[i], bulgedH[i], TT, kT);
604 pf->expinternal[i] = RESCALE_BF(internal_loop37[i], internal_loopdH[i], TT, kT);
605 }
606
607 /* special case of size 2 interior loops (single mismatch) */
608 if (james_rule)
609 pf->expinternal[2] = exp(-80 * 10. / kT);
610
611 GT = RESCALE_dG(bulge37[30],
612 bulgedH[30],
613 TT);
614 for (i = 31; i <= MAXLOOP; i++)
615 pf->expbulge[i] = exp(-TRUNC_MAYBE(GT + (pf->lxc * log(i / 30.))) * 10. / kT);
616
617 GT = RESCALE_dG(internal_loop37[30],
618 internal_loopdH[30],
619 TT);
620 for (i = 31; i <= MAXLOOP; i++)
621 pf->expinternal[i] = exp(-TRUNC_MAYBE(GT + (pf->lxc * log(i / 30.))) * 10. / kT);
622
623 GT = RESCALE_dG(ninio37, niniodH, TT);
624 for (j = 0; j <= MAXLOOP; j++)
625 pf->expninio[2][j] = exp(-MIN2(MAX_NINIO, j * TRUNC_MAYBE(GT)) * 10. / kT);
626
627 for (i = 0; (i * 7) < strlen(Tetraloops); i++)
628 pf->exptetra[i] = RESCALE_BF(Tetraloop37[i], TetraloopdH[i], TT, kT);
629
630 for (i = 0; (i * 5) < strlen(Triloops); i++)
631 pf->exptri[i] = RESCALE_BF(Triloop37[i], TriloopdH[i], TT, kT);
632
633 for (i = 0; (i * 9) < strlen(Hexaloops); i++)
634 pf->exphex[i] = RESCALE_BF(Hexaloop37[i], HexaloopdH[i], TT, kT);
635
636 for (i = 0; i <= NBPAIRS; i++)
637 pf->expMLintern[i] = RESCALE_BF(ML_intern37, ML_interndH, TT, kT);
638
639 /* if dangles==0 just set their energy to 0,
640 * don't let dangle energies become > 0 (at large temps),
641 * but make sure go smoothly to 0 */
642 for (i = 0; i <= NBPAIRS; i++)
643 for (j = 0; j <= 4; j++) {
644 if (md->dangles) {
645 pf->expdangle5[i][j] = RESCALE_BF_SMOOTH(dangle5_37[i][j], dangle5_dH[i][j], TT, kT);
646 pf->expdangle3[i][j] = RESCALE_BF_SMOOTH(dangle3_37[i][j], dangle3_dH[i][j], TT, kT);
647 } else {
648 pf->expdangle3[i][j] = pf->expdangle5[i][j] = 1;
649 }
650 }
651
652 /* stacking energies */
653 for (i = 0; i <= NBPAIRS; i++)
654 for (j = 0; j <= NBPAIRS; j++)
655 pf->expstack[i][j] = RESCALE_BF(stack37[i][j], stackdH[i][j], TT, kT);
656
657 /* mismatch energies */
658 for (i = 0; i <= NBPAIRS; i++)
659 for (j = 0; j < 5; j++)
660 for (k = 0; k < 5; k++) {
661 pf->expmismatchI[i][j][k] = RESCALE_BF(mismatchI37[i][j][k],
662 mismatchIdH[i][j][k],
663 TT,
664 kT);
665 pf->expmismatch1nI[i][j][k] = RESCALE_BF(mismatch1nI37[i][j][k],
666 mismatch1nIdH[i][j][k],
667 TT,
668 kT);
669 pf->expmismatchH[i][j][k] = RESCALE_BF(mismatchH37[i][j][k],
670 mismatchHdH[i][j][k],
671 TT,
672 kT);
673 pf->expmismatch23I[i][j][k] = RESCALE_BF(mismatch23I37[i][j][k],
674 mismatch23IdH[i][j][k],
675 TT,
676 kT);
677
678 if (md->dangles) {
679 pf->expmismatchM[i][j][k] = RESCALE_BF_SMOOTH(mismatchM37[i][j][k],
680 mismatchMdH[i][j][k],
681 TT,
682 kT);
683 pf->expmismatchExt[i][j][k] = RESCALE_BF_SMOOTH(mismatchExt37[i][j][k],
684 mismatchExtdH[i][j][k],
685 TT,
686 kT);
687 } else {
688 pf->expmismatchM[i][j][k] = pf->expmismatchExt[i][j][k] = 1.;
689 }
690 }
691
692 /* interior lops of length 2 */
693 for (i = 0; i <= NBPAIRS; i++)
694 for (j = 0; j <= NBPAIRS; j++)
695 for (k = 0; k < 5; k++)
696 for (l = 0; l < 5; l++) {
697 pf->expint11[i][j][k][l] = RESCALE_BF(int11_37[i][j][k][l],
698 int11_dH[i][j][k][l],
699 TT,
700 kT);
701 }
702
703 /* interior 2x1 loops */
704 for (i = 0; i <= NBPAIRS; i++)
705 for (j = 0; j <= NBPAIRS; j++)
706 for (k = 0; k < 5; k++)
707 for (l = 0; l < 5; l++) {
708 int m;
709 for (m = 0; m < 5; m++) {
710 pf->expint21[i][j][k][l][m] = RESCALE_BF(int21_37[i][j][k][l][m],
711 int21_dH[i][j][k][l][m],
712 TT,
713 kT);
714 }
715 }
716
717 /* interior 2x2 loops */
718 for (i = 0; i <= NBPAIRS; i++)
719 for (j = 0; j <= NBPAIRS; j++)
720 for (k = 0; k < 5; k++)
721 for (l = 0; l < 5; l++) {
722 int m, n;
723 for (m = 0; m < 5; m++)
724 for (n = 0; n < 5; n++) {
725 pf->expint22[i][j][k][l][m][n] = RESCALE_BF(int22_37[i][j][k][l][m][n],
726 int22_dH[i][j][k][l][m][n],
727 TT,
728 kT);
729 }
730 }
731
732 strncpy(pf->Tetraloops, Tetraloops, 281);
733 strncpy(pf->Triloops, Triloops, 241);
734 strncpy(pf->Hexaloops, Hexaloops, 361);
735
736 return pf;
737 }
738
739
740 PRIVATE vrna_exp_param_t *
get_exp_params_ali(vrna_md_t * md,unsigned int n_seq,double pfs)741 get_exp_params_ali(vrna_md_t *md,
742 unsigned int n_seq,
743 double pfs)
744 {
745 /* scale energy parameters and pre-calculate Boltzmann weights */
746 unsigned int i, j, k, l;
747 int pf_smooth;
748 double kTn, TT;
749 double GT;
750 vrna_exp_param_t *pf;
751
752 pf = (vrna_exp_param_t *)vrna_alloc(sizeof(vrna_exp_param_t));
753 pf->model_details = *md;
754 pf->alpha = md->betaScale;
755 pf->temperature = md->temperature;
756 pf->pf_scale = pfs;
757 pf->kT = kTn = ((double)n_seq) * md->betaScale * (md->temperature + K0) * GASCONST; /* kT in cal/mol */
758 pf_smooth = md->pf_smooth;
759 TT = (md->temperature + K0) / (Tmeasure);
760
761 pf->lxc = lxc37 * TT;
762 pf->expDuplexInit = RESCALE_BF(DuplexInit37, DuplexInitdH, TT, kTn);
763 pf->expTermAU = RESCALE_BF(TerminalAU37, TerminalAUdH, TT, kTn);
764 pf->expMLbase = RESCALE_BF(ML_BASE37, ML_BASEdH, TT, kTn / n_seq);
765 pf->expMLclosing = RESCALE_BF(ML_closing37, ML_closingdH, TT, kTn);
766 pf->expgquadLayerMismatch = RESCALE_BF(GQuadLayerMismatch37, GQuadLayerMismatchH, TT, kTn);
767 pf->gquadLayerMismatchMax = GQuadLayerMismatchMax;
768
769 for (i = VRNA_GQUAD_MIN_STACK_SIZE; i <= VRNA_GQUAD_MAX_STACK_SIZE; i++)
770 for (j = 3 * VRNA_GQUAD_MIN_LINKER_LENGTH; j <= 3 * VRNA_GQUAD_MAX_LINKER_LENGTH; j++) {
771 double GQuadAlpha_T = RESCALE_dG(GQuadAlpha37, GQuadAlphadH, TT);
772 double GQuadBeta_T = RESCALE_dG(GQuadBeta37, GQuadBetadH, TT);
773 GT = ((double)GQuadAlpha_T) * ((double)(i - 1)) + ((double)GQuadBeta_T) *
774 log(((double)j) - 2.);
775 pf->expgquad[i][j] = exp(-TRUNC_MAYBE(GT) * 10. / kTn);
776 }
777
778 /* loop energies: hairpins, bulges, interior, mulit-loops */
779 for (i = 0; i < 31; i++)
780 pf->exphairpin[i] = RESCALE_BF(hairpin37[i], hairpindH[i], TT, kTn);
781 /*add penalty for too short hairpins*/
782 for (i = 0; i < 3; i++) {
783 GT = 600 /*Penalty*/ * TT;
784 pf->exphairpin[i] = exp(-GT * 10. / kTn);
785 }
786
787 for (i = 0; i <= MIN2(30, MAXLOOP); i++) {
788 pf->expbulge[i] = RESCALE_BF(bulge37[i], bulgedH[i], TT, kTn);
789 pf->expinternal[i] = RESCALE_BF(internal_loop37[i], internal_loopdH[i], TT, kTn);
790 }
791
792 /* special case of size 2 interior loops (single mismatch) */
793 if (james_rule)
794 pf->expinternal[2] = exp(-80 * 10. / kTn);
795
796 GT = RESCALE_dG(bulge37[30], bulgedH[30], TT);
797 for (i = 31; i <= MAXLOOP; i++)
798 pf->expbulge[i] = exp(-(GT + (pf->lxc * log(i / 30.))) * 10. / kTn);
799
800 GT = RESCALE_dG(internal_loop37[30], internal_loopdH[30], TT);
801 for (i = 31; i <= MAXLOOP; i++)
802 pf->expinternal[i] = exp(-(GT + (pf->lxc * log(i / 30.))) * 10. / kTn);
803
804 GT = RESCALE_dG(ninio37, niniodH, TT);
805 for (j = 0; j <= MAXLOOP; j++)
806 pf->expninio[2][j] = exp(-MIN2(MAX_NINIO, j * GT) * 10. / kTn);
807
808 for (i = 0; (i * 7) < strlen(Tetraloops); i++)
809 pf->exptetra[i] = RESCALE_BF(Tetraloop37[i], TetraloopdH[i], TT, kTn);
810
811 for (i = 0; (i * 5) < strlen(Triloops); i++)
812 pf->exptri[i] = RESCALE_BF(Triloop37[i], TriloopdH[i], TT, kTn);
813
814 for (i = 0; (i * 9) < strlen(Hexaloops); i++)
815 pf->exphex[i] = RESCALE_BF(Hexaloop37[i], HexaloopdH[i], TT, kTn);
816
817 for (i = 0; i <= NBPAIRS; i++)
818 /* includes AU penalty */
819 pf->expMLintern[i] = RESCALE_BF(ML_intern37, ML_interndH, TT, kTn);
820
821 /* if dangle_model==0 just set their energy to 0,
822 * don't let dangle energies become > 0 (at large temps),
823 * but make sure go smoothly to 0 */
824 for (i = 0; i <= NBPAIRS; i++)
825 for (j = 0; j <= 4; j++) {
826 if (md->dangles) {
827 pf->expdangle5[i][j] = RESCALE_BF_SMOOTH(dangle5_37[i][j],
828 dangle5_dH[i][j],
829 TT,
830 kTn);
831 pf->expdangle3[i][j] = RESCALE_BF_SMOOTH(dangle3_37[i][j],
832 dangle3_dH[i][j],
833 TT,
834 kTn);
835 } else {
836 pf->expdangle3[i][j] = pf->expdangle5[i][j] = 1;
837 }
838 }
839
840 /* stacking energies */
841 for (i = 0; i <= NBPAIRS; i++)
842 for (j = 0; j <= NBPAIRS; j++) {
843 pf->expstack[i][j] = RESCALE_BF(stack37[i][j],
844 stackdH[i][j],
845 TT,
846 kTn);
847 }
848
849 /* mismatch energies */
850 for (i = 0; i <= NBPAIRS; i++)
851 for (j = 0; j < 5; j++)
852 for (k = 0; k < 5; k++) {
853 pf->expmismatchI[i][j][k] = RESCALE_BF(mismatchI37[i][j][k],
854 mismatchIdH[i][j][k],
855 TT,
856 kTn);
857 pf->expmismatch1nI[i][j][k] = RESCALE_BF(mismatch1nI37[i][j][k],
858 mismatch1nIdH[i][j][k],
859 TT,
860 kTn);
861 pf->expmismatchH[i][j][k] = RESCALE_BF(mismatchH37[i][j][k],
862 mismatchHdH[i][j][k],
863 TT,
864 kTn);
865 pf->expmismatch23I[i][j][k] = RESCALE_BF(mismatch23I37[i][j][k],
866 mismatch23IdH[i][j][k],
867 TT,
868 kTn);
869
870 if (md->dangles) {
871 pf->expmismatchM[i][j][k] = RESCALE_BF_SMOOTH(mismatchM37[i][j][k],
872 mismatchMdH[i][j][k],
873 TT,
874 kTn);
875 pf->expmismatchExt[i][j][k] = RESCALE_BF_SMOOTH(mismatchExt37[i][j][k],
876 mismatchExtdH[i][j][k],
877 TT,
878 kTn);
879 } else {
880 pf->expmismatchM[i][j][k] = pf->expmismatchExt[i][j][k] = 1.;
881 }
882 }
883
884 /* interior lops of length 2 */
885 for (i = 0; i <= NBPAIRS; i++)
886 for (j = 0; j <= NBPAIRS; j++)
887 for (k = 0; k < 5; k++)
888 for (l = 0; l < 5; l++) {
889 pf->expint11[i][j][k][l] = RESCALE_BF(int11_37[i][j][k][l],
890 int11_dH[i][j][k][l],
891 TT,
892 kTn);
893 }
894
895 /* interior 2x1 loops */
896 for (i = 0; i <= NBPAIRS; i++)
897 for (j = 0; j <= NBPAIRS; j++)
898 for (k = 0; k < 5; k++)
899 for (l = 0; l < 5; l++) {
900 int m;
901 for (m = 0; m < 5; m++) {
902 pf->expint21[i][j][k][l][m] = RESCALE_BF(int21_37[i][j][k][l][m],
903 int21_dH[i][j][k][l][m],
904 TT,
905 kTn);
906 }
907 }
908
909 /* interior 2x2 loops */
910 for (i = 0; i <= NBPAIRS; i++)
911 for (j = 0; j <= NBPAIRS; j++)
912 for (k = 0; k < 5; k++)
913 for (l = 0; l < 5; l++) {
914 int m, n;
915 for (m = 0; m < 5; m++)
916 for (n = 0; n < 5; n++) {
917 pf->expint22[i][j][k][l][m][n] = RESCALE_BF(int22_37[i][j][k][l][m][n],
918 int22_dH[i][j][k][l][m][n],
919 TT,
920 kTn);
921 }
922 }
923
924 strncpy(pf->Tetraloops, Tetraloops, 281);
925 strncpy(pf->Triloops, Triloops, 241);
926 strncpy(pf->Hexaloops, Hexaloops, 361);
927
928 return pf;
929 }
930
931
932 PRIVATE void
rescale_params(vrna_fold_compound_t * vc)933 rescale_params(vrna_fold_compound_t *vc)
934 {
935 int i;
936 vrna_exp_param_t *pf = vc->exp_params;
937 vrna_mx_pf_t *m = vc->exp_matrices;
938
939 if (m && pf) {
940 m->scale[0] = 1.;
941 m->scale[1] = (FLT_OR_DBL)(1. / pf->pf_scale);
942 m->expMLbase[0] = 1;
943 m->expMLbase[1] = (FLT_OR_DBL)(pf->expMLbase / pf->pf_scale);
944 for (i = 2; i <= vc->length; i++) {
945 m->scale[i] = m->scale[i / 2] * m->scale[i - (i / 2)];
946 m->expMLbase[i] = (FLT_OR_DBL)pow(pf->expMLbase, (double)i) * m->scale[i];
947 }
948 }
949 }
950
951
952 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
953
954 /*
955 *###########################################
956 *# deprecated functions below #
957 *###########################################
958 */
959 PUBLIC vrna_param_t *
scale_parameters(void)960 scale_parameters(void)
961 {
962 vrna_md_t md;
963
964 set_model_details(&md);
965 return vrna_params(&md);
966 }
967
968
969 PUBLIC vrna_param_t *
get_scaled_parameters(double temp,vrna_md_t md)970 get_scaled_parameters(double temp,
971 vrna_md_t md)
972 {
973 md.temperature = temp;
974 return get_scaled_params(&md);
975 }
976
977
978 PUBLIC vrna_exp_param_t *
get_boltzmann_factors(double temp,double betaScale,vrna_md_t md,double pfs)979 get_boltzmann_factors(double temp,
980 double betaScale,
981 vrna_md_t md,
982 double pfs)
983 {
984 md.temperature = temp;
985 md.betaScale = betaScale;
986 pf_scale = pfs;
987
988 return get_scaled_exp_params(&md, pfs);
989 }
990
991
992 PUBLIC vrna_exp_param_t *
get_scaled_pf_parameters(void)993 get_scaled_pf_parameters(void)
994 {
995 vrna_md_t md;
996 vrna_exp_param_t *pf;
997
998 set_model_details(&md);
999
1000 pf = vrna_exp_params(&md);
1001 pf->pf_scale = pf_scale;
1002
1003 return pf;
1004 }
1005
1006
1007 PUBLIC vrna_exp_param_t *
get_boltzmann_factors_ali(unsigned int n_seq,double temp,double betaScale,vrna_md_t md,double pfs)1008 get_boltzmann_factors_ali(unsigned int n_seq,
1009 double temp,
1010 double betaScale,
1011 vrna_md_t md,
1012 double pfs)
1013 {
1014 md.temperature = temp;
1015 md.betaScale = betaScale;
1016 pf_scale = pfs;
1017
1018 return get_exp_params_ali(&md, n_seq, pfs);
1019 }
1020
1021
1022 PUBLIC vrna_exp_param_t *
get_scaled_alipf_parameters(unsigned int n_seq)1023 get_scaled_alipf_parameters(unsigned int n_seq)
1024 {
1025 vrna_md_t md;
1026
1027 set_model_details(&md);
1028
1029 return get_exp_params_ali(&md, n_seq, pf_scale);
1030 }
1031
1032
1033 PUBLIC vrna_exp_param_t *
get_boltzmann_factor_copy(vrna_exp_param_t * par)1034 get_boltzmann_factor_copy(vrna_exp_param_t *par)
1035 {
1036 return vrna_exp_params_copy(par);
1037 }
1038
1039
1040 PUBLIC vrna_param_t *
get_parameter_copy(vrna_param_t * par)1041 get_parameter_copy(vrna_param_t *par)
1042 {
1043 return vrna_params_copy(par);
1044 }
1045
1046
1047 PUBLIC vrna_param_t *
copy_parameters(void)1048 copy_parameters(void)
1049 {
1050 vrna_param_t *copy;
1051
1052 if (p.id != id) {
1053 vrna_md_t md;
1054 set_model_details(&md);
1055 return vrna_params(&md);
1056 } else {
1057 copy = (vrna_param_t *)vrna_alloc(sizeof(vrna_param_t));
1058 memcpy(copy, &p, sizeof(vrna_param_t));
1059 }
1060
1061 return copy;
1062 }
1063
1064
1065 PUBLIC vrna_param_t *
set_parameters(vrna_param_t * dest)1066 set_parameters(vrna_param_t *dest)
1067 {
1068 memcpy(&p, dest, sizeof(vrna_param_t));
1069 return &p;
1070 }
1071
1072
1073 PUBLIC vrna_exp_param_t *
copy_pf_param(void)1074 copy_pf_param(void)
1075 {
1076 vrna_exp_param_t *copy;
1077
1078 if (pf.id != pf_id) {
1079 vrna_md_t md;
1080 set_model_details(&md);
1081 copy = vrna_exp_params(&md);
1082 copy->pf_scale = pf_scale;
1083 return copy;
1084 } else {
1085 copy = (vrna_exp_param_t *)vrna_alloc(sizeof(vrna_exp_param_t));
1086 memcpy(copy, &pf, sizeof(vrna_exp_param_t));
1087 }
1088
1089 return copy;
1090 }
1091
1092
1093 PUBLIC vrna_exp_param_t *
set_pf_param(vrna_param_t * dest)1094 set_pf_param(vrna_param_t *dest)
1095 {
1096 memcpy(&pf, dest, sizeof(vrna_exp_param_t));
1097 return &pf;
1098 }
1099
1100
1101 PUBLIC vrna_exp_param_t *
scale_pf_parameters(void)1102 scale_pf_parameters(void)
1103 {
1104 vrna_md_t md;
1105 vrna_exp_param_t *pf;
1106
1107 set_model_details(&md);
1108
1109 pf = vrna_exp_params(&md);
1110 pf->pf_scale = pf_scale;
1111
1112 return pf;
1113 }
1114
1115
1116 #endif
1117