RNAlib-2.1.9
loop_energies.h
Go to the documentation of this file.
1 #ifndef __VIENNA_RNA_PACKAGE_LOOP_ENERGIES_H__
2 #define __VIENNA_RNA_PACKAGE_LOOP_ENERGIES_H__
3 
4 #include <stdio.h>
5 #include <stdlib.h>
6 #include <math.h>
7 #include <ctype.h>
8 #include <string.h>
9 #include "params.h"
10 #include "fold_vars.h"
11 #include "energy_par.h"
12 
13 #ifdef __GNUC__
14 # define INLINE inline
15 #else
16 # define INLINE
17 #endif
18 
54 INLINE PRIVATE int E_MLstem( int type,
55  int si1,
56  int sj1,
57  paramT *P);
58 
65 INLINE PRIVATE double exp_E_MLstem(int type,
66  int si1,
67  int sj1,
68  pf_paramT *P);
69 
89 INLINE PRIVATE int E_ExtLoop(int type,
90  int si1,
91  int sj1,
92  paramT *P);
93 
100 INLINE PRIVATE double exp_E_ExtLoop( int type,
101  int si1,
102  int sj1,
103  pf_paramT *P);
104 
149 INLINE PRIVATE int E_IntLoop(int n1,
150  int n2,
151  int type,
152  int type_2,
153  int si1,
154  int sj1,
155  int sp1,
156  int sq1,
157  paramT *P);
158 
159 
191 INLINE PRIVATE int E_Hairpin(int size,
192  int type,
193  int si1,
194  int sj1,
195  const char *string,
196  paramT *P);
197 
243 INLINE PRIVATE int E_Stem( int type,
244  int si1,
245  int sj1,
246  int extLoop,
247  paramT *P);
248 
257 INLINE PRIVATE double exp_E_Stem(int type,
258  int si1,
259  int sj1,
260  int extLoop,
261  pf_paramT *P);
262 
280 INLINE PRIVATE double exp_E_Hairpin( int u,
281  int type,
282  short si1,
283  short sj1,
284  const char *string,
285  pf_paramT *P);
286 
306 INLINE PRIVATE double exp_E_IntLoop(int u1,
307  int u2,
308  int type,
309  int type2,
310  short si1,
311  short sj1,
312  short sp1,
313  short sq1,
314  pf_paramT *P);
315 
316 
317 /*
318 #################################
319 # BEGIN OF FUNCTION DEFINITIONS #
320 #################################
321 */
322 INLINE PRIVATE int E_Hairpin(int size, int type, int si1, int sj1, const char *string, paramT *P){
323  int energy;
324 
325  energy = (size <= 30) ? P->hairpin[size] : P->hairpin[30]+(int)(P->lxc*log((size)/30.));
326  if (P->model_details.special_hp){
327  if (size == 4) { /* check for tetraloop bonus */
328  char tl[7]={0}, *ts;
329  strncpy(tl, string, 6);
330  if ((ts=strstr(P->Tetraloops, tl)))
331  return (P->Tetraloop_E[(ts - P->Tetraloops)/7]);
332  }
333  else if (size == 6) {
334  char tl[9]={0}, *ts;
335  strncpy(tl, string, 8);
336  if ((ts=strstr(P->Hexaloops, tl)))
337  return (energy = P->Hexaloop_E[(ts - P->Hexaloops)/9]);
338  }
339  else if (size == 3) {
340  char tl[6]={0,0,0,0,0,0}, *ts;
341  strncpy(tl, string, 5);
342  if ((ts=strstr(P->Triloops, tl))) {
343  return (P->Triloop_E[(ts - P->Triloops)/6]);
344  }
345  return (energy + (type>2 ? P->TerminalAU : 0));
346  }
347  }
348  energy += P->mismatchH[type][si1][sj1];
349 
350  return energy;
351 }
352 
353 INLINE PRIVATE int E_IntLoop(int n1, int n2, int type, int type_2, int si1, int sj1, int sp1, int sq1, paramT *P){
354  /* compute energy of degree 2 loop (stack bulge or interior) */
355  int nl, ns, energy;
356  energy = INF;
357 
358  if (n1>n2) { nl=n1; ns=n2;}
359  else {nl=n2; ns=n1;}
360 
361  if (nl == 0)
362  return P->stack[type][type_2]; /* stack */
363 
364  if (ns==0) { /* bulge */
365  energy = (nl<=MAXLOOP)?P->bulge[nl]:
366  (P->bulge[30]+(int)(P->lxc*log(nl/30.)));
367  if (nl==1) energy += P->stack[type][type_2];
368  else {
369  if (type>2) energy += P->TerminalAU;
370  if (type_2>2) energy += P->TerminalAU;
371  }
372  return energy;
373  }
374  else { /* interior loop */
375  if (ns==1) {
376  if (nl==1) /* 1x1 loop */
377  return P->int11[type][type_2][si1][sj1];
378  if (nl==2) { /* 2x1 loop */
379  if (n1==1)
380  energy = P->int21[type][type_2][si1][sq1][sj1];
381  else
382  energy = P->int21[type_2][type][sq1][si1][sp1];
383  return energy;
384  }
385  else { /* 1xn loop */
386  energy = (nl+1<=MAXLOOP)?(P->internal_loop[nl+1]) : (P->internal_loop[30]+(int)(P->lxc*log((nl+1)/30.)));
387  energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
388  energy += P->mismatch1nI[type][si1][sj1] + P->mismatch1nI[type_2][sq1][sp1];
389  return energy;
390  }
391  }
392  else if (ns==2) {
393  if(nl==2) { /* 2x2 loop */
394  return P->int22[type][type_2][si1][sp1][sq1][sj1];}
395  else if (nl==3){ /* 2x3 loop */
396  energy = P->internal_loop[5]+P->ninio[2];
397  energy += P->mismatch23I[type][si1][sj1] + P->mismatch23I[type_2][sq1][sp1];
398  return energy;
399  }
400 
401  }
402  { /* generic interior loop (no else here!)*/
403  energy = (n1+n2<=MAXLOOP)?(P->internal_loop[n1+n2]) : (P->internal_loop[30]+(int)(P->lxc*log((n1+n2)/30.)));
404 
405  energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
406 
407  energy += P->mismatchI[type][si1][sj1] + P->mismatchI[type_2][sq1][sp1];
408  }
409  }
410  return energy;
411 }
412 
413 INLINE PRIVATE int E_Stem(int type, int si1, int sj1, int extLoop, paramT *P){
414  int energy = 0;
415  int d5 = (si1 >= 0) ? P->dangle5[type][si1] : 0;
416  int d3 = (sj1 >= 0) ? P->dangle3[type][sj1] : 0;
417 
418  if(type > 2)
419  energy += P->TerminalAU;
420 
421  if(si1 >= 0 && sj1 >= 0)
422  energy += (extLoop) ? P->mismatchExt[type][si1][sj1] : P->mismatchM[type][si1][sj1];
423  else
424  energy += d5 + d3;
425 
426  if(!extLoop) energy += P->MLintern[type];
427  return energy;
428 }
429 
430 INLINE PRIVATE int E_ExtLoop(int type, int si1, int sj1, paramT *P){
431  int energy = 0;
432  if(si1 >= 0 && sj1 >= 0){
433  energy += P->mismatchExt[type][si1][sj1];
434  }
435  else if (si1 >= 0){
436  energy += P->dangle5[type][si1];
437  }
438  else if (sj1 >= 0){
439  energy += P->dangle3[type][sj1];
440  }
441 
442  if(type > 2)
443  energy += P->TerminalAU;
444 
445  return energy;
446 }
447 
448 INLINE PRIVATE int E_MLstem(int type, int si1, int sj1, paramT *P){
449  int energy = 0;
450  if(si1 >= 0 && sj1 >= 0){
451  energy += P->mismatchM[type][si1][sj1];
452  }
453  else if (si1 >= 0){
454  energy += P->dangle5[type][si1];
455  }
456  else if (sj1 >= 0){
457  energy += P->dangle3[type][sj1];
458  }
459 
460  if(type > 2)
461  energy += P->TerminalAU;
462 
463  energy += P->MLintern[type];
464 
465  return energy;
466 }
467 
468 INLINE PRIVATE double exp_E_Hairpin(int u, int type, short si1, short sj1, const char *string, pf_paramT *P){
469  double q, kT;
470  kT = P->kT; /* kT in cal/mol */
471 
472  if(u <= 30)
473  q = P->exphairpin[u];
474  else
475  q = P->exphairpin[30] * exp( -(P->lxc*log( u/30.))*10./kT);
476 
477  if(u < 3) return q; /* should only be the case when folding alignments */
478 
479  if(P->model_details.special_hp){
480  if(u==4) {
481  char tl[7]={0,0,0,0,0,0,0}, *ts;
482  strncpy(tl, string, 6);
483  if ((ts=strstr(P->Tetraloops, tl))){
484  if(type != 7)
485  return (P->exptetra[(ts-P->Tetraloops)/7]);
486  else
487  q *= P->exptetra[(ts-P->Tetraloops)/7];
488  }
489  }
490  if (u==6) {
491  char tl[9]={0,0,0,0,0,0,0,0,0}, *ts;
492  strncpy(tl, string, 8);
493  if ((ts=strstr(P->Hexaloops, tl)))
494  return (P->exphex[(ts-P->Hexaloops)/9]);
495  }
496  if (u==3) {
497  char tl[6]={0,0,0,0,0,0}, *ts;
498  strncpy(tl, string, 5);
499  if ((ts=strstr(P->Triloops, tl)))
500  return (P->exptri[(ts-P->Triloops)/6]);
501  if (type>2)
502  return q * P->expTermAU;
503  return q;
504  }
505  }
506  /* no mismatches for tri-loops */
507  q *= P->expmismatchH[type][si1][sj1];
508 
509  return q;
510 }
511 
512 INLINE PRIVATE double exp_E_IntLoop(int u1, int u2, int type, int type2, short si1, short sj1, short sp1, short sq1, pf_paramT *P){
513  int ul, us, no_close = 0;
514  double z = 0.;
515 
516  if ((no_closingGU) && ((type2==3)||(type2==4)||(type==3)||(type==4)))
517  no_close = 1;
518 
519  if (u1>u2) { ul=u1; us=u2;}
520  else {ul=u2; us=u1;}
521 
522  if (ul==0) /* stack */
523  z = P->expstack[type][type2];
524  else if(!no_close){
525  if (us==0) { /* bulge */
526  z = P->expbulge[ul];
527  if (ul==1) z *= P->expstack[type][type2];
528  else {
529  if (type>2) z *= P->expTermAU;
530  if (type2>2) z *= P->expTermAU;
531  }
532  return z;
533  }
534  else if (us==1) {
535  if (ul==1){ /* 1x1 loop */
536  return P->expint11[type][type2][si1][sj1];
537  }
538  if (ul==2) { /* 2x1 loop */
539  if (u1==1)
540  return P->expint21[type][type2][si1][sq1][sj1];
541  else
542  return P->expint21[type2][type][sq1][si1][sp1];
543  }
544  else { /* 1xn loop */
545  z = P->expinternal[ul+us] * P->expmismatch1nI[type][si1][sj1] * P->expmismatch1nI[type2][sq1][sp1];
546  return z * P->expninio[2][ul-us];
547  }
548  }
549  else if (us==2) {
550  if(ul==2) /* 2x2 loop */
551  return P->expint22[type][type2][si1][sp1][sq1][sj1];
552  else if(ul==3){ /* 2x3 loop */
553  z = P->expinternal[5]*P->expmismatch23I[type][si1][sj1]*P->expmismatch23I[type2][sq1][sp1];
554  return z * P->expninio[2][1];
555  }
556  }
557  /* generic interior loop (no else here!)*/
558  z = P->expinternal[ul+us] * P->expmismatchI[type][si1][sj1] * P->expmismatchI[type2][sq1][sp1];
559  return z * P->expninio[2][ul-us];
560 
561  }
562  return z;
563 }
564 
565 INLINE PRIVATE double exp_E_Stem(int type, int si1, int sj1, int extLoop, pf_paramT *P){
566  double energy = 1.0;
567  double d5 = (si1 >= 0) ? P->expdangle5[type][si1] : 1.;
568  double d3 = (sj1 >= 0) ? P->expdangle3[type][sj1] : 1.;
569 
570  if(type > 2)
571  energy *= P->expTermAU;
572 
573  if(si1 >= 0 && sj1 >= 0)
574  energy *= (extLoop) ? P->expmismatchExt[type][si1][sj1] : P->expmismatchM[type][si1][sj1];
575  else
576  energy *= d5 * d3;
577 
578  if(!extLoop) energy *= P->expMLintern[type];
579  return energy;
580 }
581 
582 INLINE PRIVATE double exp_E_MLstem(int type, int si1, int sj1, pf_paramT *P){
583  double energy = 1.0;
584  if(si1 >= 0 && sj1 >= 0){
585  energy *= P->expmismatchM[type][si1][sj1];
586  }
587  else if(si1 >= 0){
588  energy *= P->expdangle5[type][si1];
589  }
590  else if(sj1 >= 0){
591  energy *= P->expdangle3[type][sj1];
592  }
593 
594  if(type > 2)
595  energy *= P->expTermAU;
596 
597  energy *= P->expMLintern[type];
598  return energy;
599 }
600 
601 INLINE PRIVATE double exp_E_ExtLoop(int type, int si1, int sj1, pf_paramT *P){
602  double energy = 1.0;
603  if(si1 >= 0 && sj1 >= 0){
604  energy *= P->expmismatchExt[type][si1][sj1];
605  }
606  else if(si1 >= 0){
607  energy *= P->expdangle5[type][si1];
608  }
609  else if(sj1 >= 0){
610  energy *= P->expdangle3[type][sj1];
611  }
612 
613  if(type > 2)
614  energy *= P->expTermAU;
615 
616  return energy;
617 }
618 
619 INLINE PRIVATE int E_IntLoop_Co(int type, int type_2, int i, int j, int p, int q, int cutpoint, short si1, short sj1, short sp1, short sq1, int dangles, paramT *P){
620  int energy = 0;
621  if(type > 2) energy += P->TerminalAU;
622  if(type_2 > 2) energy += P->TerminalAU;
623 
624  if(!dangles) return energy;
625 
626  int ci = (i>=cutpoint)||((i+1)<cutpoint);
627  int cj = ((j-1)>=cutpoint)||(j<cutpoint);
628  int cp = ((p-1)>=cutpoint)||(p<cutpoint);
629  int cq = (q>=cutpoint)||((q+1)<cutpoint);
630 
631  int d3 = ci ? P->dangle3[type][si1] : 0;
632  int d5 = cj ? P->dangle5[type][sj1] : 0;
633  int d5_2 = cp ? P->dangle5[type_2][sp1] : 0;
634  int d3_2 = cq ? P->dangle3[type_2][sq1] : 0;
635 
636  int tmm = (cj && ci) ? P->mismatchExt[type][sj1][si1] : d5 + d3;
637  int tmm_2 = (cp && cq) ? P->mismatchExt[type_2][sp1][sq1] : d5_2 + d3_2;
638 
639  if(dangles == 2) return energy + tmm + tmm_2;
640 
641  /* now we may have non-double dangles only */
642  if(i+2 < p){
643  if(q+2 < j){ energy += tmm + tmm_2;}
644  else if(q+2 == j){ energy += (cj && cq) ? MIN2(tmm + d5_2, tmm_2 + d3) : tmm + tmm_2;}
645  else energy += d3 + d5_2;
646  }
647  else if(i+2 == p){
648  if(q+2 < j){ energy += (ci && cp) ? MIN2(tmm + d3_2, tmm_2 + d5) : tmm + tmm_2;}
649  else if(q+2 == j){
650  energy += MIN2(tmm, MIN2(tmm_2, MIN2(d5 + d5_2, d3 + d3_2)));
651  }
652  else energy += MIN2(d3, d5_2);
653  }
654  else{
655  if(q+2 < j){ energy += d5 + d3_2;}
656  else if(q+2 == j){ energy += MIN2(d5, d3_2);}
657  }
658  return energy;
659 }
660 
661 #endif