- Main Page
- Related Pages
- Modules
- Data Structures
- Files

- File List
- Globals

# /homes/fmmarquesmadeira/Projects/jabaws/binaries/src/ViennaRNA/H/loop_energies.h

Go to the documentation of this file.

    00001 #ifndef __VIENNA_RNA_PACKAGE_LOOP_ENERGIES_H__
    00002 #define __VIENNA_RNA_PACKAGE_LOOP_ENERGIES_H__
    00003
    00004 #include <stdio.h>
    00005 #include <stdlib.h>
    00006 #include <math.h>
    00007 #include <ctype.h>
    00008 #include <string.h>
    00009 #include "params.h"
    00010 #include "fold_vars.h"
    00011 #include "energy_par.h"
    00012
    00013 #ifdef __GNUC__
    00014 # define INLINE inline
    00015 #else
    00016 # define INLINE
    00017 #endif
    00018
    00054 INLINE  PRIVATE int E_MLstem( int type,
    00055                               int si1,
    00056                               int sj1,
    00057                               paramT *P);
    00058
    00065 INLINE  PRIVATE double exp_E_MLstem(int type,
    00066                                     int si1,
    00067                                     int sj1,
    00068                                     pf_paramT *P);
    00069
    00089 INLINE  PRIVATE int E_ExtLoop(int type,
    00090                               int si1,
    00091                               int sj1,
    00092                               paramT *P);
    00093
    00100 INLINE  PRIVATE double exp_E_ExtLoop( int type,
    00101                                       int si1,
    00102                                       int sj1,
    00103                                       pf_paramT *P);
    00104
    00149 INLINE  PRIVATE int E_IntLoop(int n1,
    00150                               int n2,
    00151                               int type,
    00152                               int type_2,
    00153                               int si1,
    00154                               int sj1,
    00155                               int sp1,
    00156                               int sq1,
    00157                               paramT *P);
    00158
    00159
    00191 INLINE  PRIVATE int E_Hairpin(int size,
    00192                               int type,
    00193                               int si1,
    00194                               int sj1,
    00195                               const char *string,
    00196                               paramT *P);
    00197
    00243 INLINE  PRIVATE int E_Stem( int type,
    00244                             int si1,
    00245                             int sj1,
    00246                             int extLoop,
    00247                             paramT *P);
    00248
    00257 INLINE  PRIVATE double exp_E_Stem(int type,
    00258                                   int si1,
    00259                                   int sj1,
    00260                                   int extLoop,
    00261                                   pf_paramT *P);
    00262
    00280 INLINE  PRIVATE double exp_E_Hairpin( int u,
    00281                                       int type,
    00282                                       short si1,
    00283                                       short sj1,
    00284                                       const char *string,
    00285                                       pf_paramT *P);
    00286
    00306 INLINE  PRIVATE double  exp_E_IntLoop(int u1,
    00307                                       int u2,
    00308                                       int type,
    00309                                       int type2,
    00310                                       short si1,
    00311                                       short sj1,
    00312                                       short sp1,
    00313                                       short sq1,
    00314                                       pf_paramT *P);
    00315
    00316
    00317 /*
    00318 #################################
    00319 # BEGIN OF FUNCTION DEFINITIONS #
    00320 #################################
    00321 */
    00322 INLINE  PRIVATE int E_Hairpin(int size, int type, int si1, int sj1, const char *string, paramT *P){
    00323   int energy;
    00324
    00325   energy = (size <= 30) ? P->hairpin[size] : P->hairpin[30]+(int)(P->lxc*log((size)/30.));
    00326   if (P->model_details.special_hp){
    00327     if (size == 4) { /* check for tetraloop bonus */
    00328       char tl[7]={0}, *ts;
    00329       strncpy(tl, string, 6);
    00330       if ((ts=strstr(P->Tetraloops, tl)))
    00331         return (P->Tetraloop_E[(ts - P->Tetraloops)/7]);
    00332     }
    00333     else if (size == 6) {
    00334       char tl[9]={0}, *ts;
    00335       strncpy(tl, string, 8);
    00336       if ((ts=strstr(P->Hexaloops, tl)))
    00337         return (energy = P->Hexaloop_E[(ts - P->Hexaloops)/9]);
    00338     }
    00339     else if (size == 3) {
    00340       char tl[6]={0,0,0,0,0,0}, *ts;
    00341       strncpy(tl, string, 5);
    00342       if ((ts=strstr(P->Triloops, tl))) {
    00343         return (P->Triloop_E[(ts - P->Triloops)/6]);
    00344       }
    00345       return (energy + (type>2 ? P->TerminalAU : 0));
    00346     }
    00347   }
    00348   energy += P->mismatchH[type][si1][sj1];
    00349
    00350   return energy;
    00351 }
    00352
    00353 INLINE  PRIVATE int E_IntLoop(int n1, int n2, int type, int type_2, int si1, int sj1, int sp1, int sq1, paramT *P){
    00354   /* compute energy of degree 2 loop (stack bulge or interior) */
    00355   int nl, ns, energy;
    00356   energy = INF;
    00357
    00358   if (n1>n2) { nl=n1; ns=n2;}
    00359   else {nl=n2; ns=n1;}
    00360
    00361   if (nl == 0)
    00362     return P->stack[type][type_2];  /* stack */
    00363
    00364   if (ns==0) {                      /* bulge */
    00365     energy = (nl<=MAXLOOP)?P->bulge[nl]:
    00366       (P->bulge[30]+(int)(P->lxc*log(nl/30.)));
    00367     if (nl==1) energy += P->stack[type][type_2];
    00368     else {
    00369       if (type>2) energy += P->TerminalAU;
    00370       if (type_2>2) energy += P->TerminalAU;
    00371     }
    00372     return energy;
    00373   }
    00374   else {                            /* interior loop */
    00375     if (ns==1) {
    00376       if (nl==1)                    /* 1x1 loop */
    00377         return P->int11[type][type_2][si1][sj1];
    00378       if (nl==2) {                  /* 2x1 loop */
    00379         if (n1==1)
    00380           energy = P->int21[type][type_2][si1][sq1][sj1];
    00381         else
    00382           energy = P->int21[type_2][type][sq1][si1][sp1];
    00383         return energy;
    00384       }
    00385       else {  /* 1xn loop */
    00386         energy = (nl+1<=MAXLOOP)?(P->internal_loop[nl+1]) : (P->internal_loop[30]+(int)(P->lxc*log((nl+1)/30.)));
    00387         energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
    00388         energy += P->mismatch1nI[type][si1][sj1] + P->mismatch1nI[type_2][sq1][sp1];
    00389         return energy;
    00390       }
    00391     }
    00392     else if (ns==2) {
    00393       if(nl==2)      {              /* 2x2 loop */
    00394         return P->int22[type][type_2][si1][sp1][sq1][sj1];}
    00395       else if (nl==3){              /* 2x3 loop */
    00396         energy = P->internal_loop[5]+P->ninio[2];
    00397         energy += P->mismatch23I[type][si1][sj1] + P->mismatch23I[type_2][sq1][sp1];
    00398         return energy;
    00399       }
    00400
    00401     }
    00402     { /* generic interior loop (no else here!)*/
    00403       energy = (n1+n2<=MAXLOOP)?(P->internal_loop[n1+n2]) : (P->internal_loop[30]+(int)(P->lxc*log((n1+n2)/30.)));
    00404
    00405       energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
    00406
    00407       energy += P->mismatchI[type][si1][sj1] + P->mismatchI[type_2][sq1][sp1];
    00408     }
    00409   }
    00410   return energy;
    00411 }
    00412
    00413 INLINE  PRIVATE int E_Stem(int type, int si1, int sj1, int extLoop, paramT *P){
    00414   int energy = 0;
    00415   int d5 = (si1 >= 0) ? P->dangle5[type][si1] : 0;
    00416   int d3 = (sj1 >= 0) ? P->dangle3[type][sj1] : 0;
    00417
    00418   if(type > 2)
    00419     energy += P->TerminalAU;
    00420
    00421   if(si1 >= 0 && sj1 >= 0)
    00422     energy += (extLoop) ? P->mismatchExt[type][si1][sj1] : P->mismatchM[type][si1][sj1];
    00423   else
    00424     energy += d5 + d3;
    00425
    00426   if(!extLoop) energy += P->MLintern[type];
    00427   return energy;
    00428 }
    00429
    00430 INLINE  PRIVATE int E_ExtLoop(int type, int si1, int sj1, paramT *P){
    00431   int energy = 0;
    00432   if(si1 >= 0 && sj1 >= 0){
    00433     energy += P->mismatchExt[type][si1][sj1];
    00434   }
    00435   else if (si1 >= 0){
    00436     energy += P->dangle5[type][si1];
    00437   }
    00438   else if (sj1 >= 0){
    00439     energy += P->dangle3[type][sj1];
    00440   }
    00441
    00442   if(type > 2)
    00443     energy += P->TerminalAU;
    00444
    00445   return energy;
    00446 }
    00447
    00448 INLINE  PRIVATE int E_MLstem(int type, int si1, int sj1, paramT *P){
    00449   int energy = 0;
    00450   if(si1 >= 0 && sj1 >= 0){
    00451     energy += P->mismatchM[type][si1][sj1];
    00452   }
    00453   else if (si1 >= 0){
    00454     energy += P->dangle5[type][si1];
    00455   }
    00456   else if (sj1 >= 0){
    00457     energy += P->dangle3[type][sj1];
    00458   }
    00459
    00460   if(type > 2)
    00461     energy += P->TerminalAU;
    00462
    00463   energy += P->MLintern[type];
    00464
    00465   return energy;
    00466 }
    00467
    00468 INLINE  PRIVATE double exp_E_Hairpin(int u, int type, short si1, short sj1, const char *string, pf_paramT *P){
    00469   double q, kT;
    00470   kT = P->kT;   /* kT in cal/mol  */
    00471
    00472   if(u <= 30)
    00473     q = P->exphairpin[u];
    00474   else
    00475     q = P->exphairpin[30] * exp( -(P->lxc*log( u/30.))*10./kT);
    00476
    00477   if(u < 3) return q; /* should only be the case when folding alignments */
    00478
    00479   if ((P->model_details.special_hp)&&(u==4)) {
    00480     char tl[7]={0,0,0,0,0,0,0}, *ts;
    00481     strncpy(tl, string, 6);
    00482     if ((ts=strstr(P->Tetraloops, tl))){
    00483       if(type != 7)
    00484         return (P->exptetra[(ts-P->Tetraloops)/7]);
    00485       else
    00486         q *= P->exptetra[(ts-P->Tetraloops)/7];
    00487     }
    00488   }
    00489   if ((tetra_loop)&&(u==6)) {
    00490     char tl[9]={0,0,0,0,0,0,0,0,0}, *ts;
    00491     strncpy(tl, string, 6);
    00492     if ((ts=strstr(P->Hexaloops, tl)))
    00493       return  (P->exphex[(ts-P->Hexaloops)/9]);
    00494   }
    00495   if (u==3) {
    00496     char tl[6]={0,0,0,0,0,0}, *ts;
    00497     strncpy(tl, string, 5);
    00498     if ((ts=strstr(P->Triloops, tl)))
    00499       return (P->exptri[(ts-P->Triloops)/6]);
    00500     if (type>2)
    00501       q *= P->expTermAU;
    00502   }
    00503   else /* no mismatches for tri-loops */
    00504     q *= P->expmismatchH[type][si1][sj1];
    00505
    00506   return q;
    00507 }
    00508
    00509 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){
    00510   int ul, us, no_close = 0;
    00511   double z = 0.;
    00512
    00513   if ((no_closingGU) && ((type2==3)||(type2==4)||(type==3)||(type==4)))
    00514     no_close = 1;
    00515
    00516   if (u1>u2) { ul=u1; us=u2;}
    00517   else {ul=u2; us=u1;}
    00518
    00519   if (ul==0) /* stack */
    00520     z = P->expstack[type][type2];
    00521   else if(!no_close){
    00522     if (us==0) {                      /* bulge */
    00523       z = P->expbulge[ul];
    00524       if (ul==1) z *= P->expstack[type][type2];
    00525       else {
    00526         if (type>2) z *= P->expTermAU;
    00527         if (type2>2) z *= P->expTermAU;
    00528       }
    00529       return z;
    00530     }
    00531     else if (us==1) {
    00532       if (ul==1){                    /* 1x1 loop */
    00533         return P->expint11[type][type2][si1][sj1];
    00534       }
    00535       if (ul==2) {                  /* 2x1 loop */
    00536         if (u1==1)
    00537           return P->expint21[type][type2][si1][sq1][sj1];
    00538         else
    00539           return P->expint21[type2][type][sq1][si1][sp1];
    00540       }
    00541       else {  /* 1xn loop */
    00542         z = P->expinternal[ul+us] * P->expmismatch1nI[type][si1][sj1] * P->expmismatch1nI[type2][sq1][sp1];
    00543         return z * P->expninio[2][ul-us];
    00544       }
    00545     }
    00546     else if (us==2) {
    00547       if(ul==2) /* 2x2 loop */
    00548         return P->expint22[type][type2][si1][sp1][sq1][sj1];
    00549       else if(ul==3){              /* 2x3 loop */
    00550         z = P->expinternal[5]*P->expmismatch23I[type][si1][sj1]*P->expmismatch23I[type2][sq1][sp1];
    00551         return z * P->expninio[2][1];
    00552       }
    00553     }
    00554     /* generic interior loop (no else here!)*/
    00555     z = P->expinternal[ul+us] * P->expmismatchI[type][si1][sj1] * P->expmismatchI[type2][sq1][sp1];
    00556     return z * P->expninio[2][ul-us];
    00557
    00558   }
    00559   return z;
    00560 }
    00561
    00562 INLINE  PRIVATE double exp_E_Stem(int type, int si1, int sj1, int extLoop, pf_paramT *P){
    00563   double energy = 1.0;
    00564   double d5 = (si1 >= 0) ? P->expdangle5[type][si1] : 1.;
    00565   double d3 = (sj1 >= 0) ? P->expdangle3[type][sj1] : 1.;
    00566
    00567   if(type > 2)
    00568     energy *= P->expTermAU;
    00569
    00570   if(si1 >= 0 && sj1 >= 0)
    00571     energy *= (extLoop) ? P->expmismatchExt[type][si1][sj1] : P->expmismatchM[type][si1][sj1];
    00572   else
    00573     energy *= d5 * d3;
    00574
    00575   if(!extLoop) energy *= P->expMLintern[type];
    00576   return energy;
    00577 }
    00578
    00579 INLINE  PRIVATE double exp_E_MLstem(int type, int si1, int sj1, pf_paramT *P){
    00580   double energy = 1.0;
    00581   if(si1 >= 0 && sj1 >= 0){
    00582     energy *= P->expmismatchM[type][si1][sj1];
    00583   }
    00584   else if(si1 >= 0){
    00585     energy *= P->expdangle5[type][si1];
    00586   }
    00587   else if(sj1 >= 0){
    00588     energy *= P->expdangle3[type][sj1];
    00589   }
    00590
    00591   if(type > 2)
    00592     energy *= P->expTermAU;
    00593
    00594   energy *= P->expMLintern[type];
    00595   return energy;
    00596 }
    00597
    00598 INLINE  PRIVATE double exp_E_ExtLoop(int type, int si1, int sj1, pf_paramT *P){
    00599   double energy = 1.0;
    00600   if(si1 >= 0 && sj1 >= 0){
    00601     energy *= P->expmismatchExt[type][si1][sj1];
    00602   }
    00603   else if(si1 >= 0){
    00604     energy *= P->expdangle5[type][si1];
    00605   }
    00606   else if(sj1 >= 0){
    00607     energy *= P->expdangle3[type][sj1];
    00608   }
    00609
    00610   if(type > 2)
    00611     energy *= P->expTermAU;
    00612
    00613   return energy;
    00614 }
    00615
    00616 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){
    00617   int energy = 0;
    00618   if(type > 2)   energy += P->TerminalAU;
    00619   if(type_2 > 2) energy += P->TerminalAU;
    00620
    00621   if(!dangles) return energy;
    00622
    00623   int ci = (i>=cutpoint)||((i+1)<cutpoint);
    00624   int cj = ((j-1)>=cutpoint)||(j<cutpoint);
    00625   int cp = ((p-1)>=cutpoint)||(p<cutpoint);
    00626   int cq = (q>=cutpoint)||((q+1)<cutpoint);
    00627
    00628   int d3    = ci  ? P->dangle3[type][si1]   : 0;
    00629   int d5    = cj  ? P->dangle5[type][sj1]   : 0;
    00630   int d5_2  = cp  ? P->dangle5[type_2][sp1] : 0;
    00631   int d3_2  = cq  ? P->dangle3[type_2][sq1] : 0;
    00632
    00633   int tmm   = (cj && ci) ? P->mismatchExt[type][sj1][si1]   : d5 + d3;
    00634   int tmm_2 = (cp && cq) ? P->mismatchExt[type_2][sp1][sq1] : d5_2 + d3_2;
    00635
    00636   if(dangles == 2) return energy + tmm + tmm_2;
    00637
    00638   /* now we may have non-double dangles only */
    00639   if(i+2 < p){
    00640     if(q+2 < j){ energy += tmm + tmm_2;}
    00641     else if(q+2 == j){ energy += (cj && cq) ? MIN2(tmm + d5_2, tmm_2 + d3) : tmm + tmm_2;}
    00642     else energy += d3 + d5_2;
    00643   }
    00644   else if(i+2 == p){
    00645     if(q+2 < j){ energy += (ci && cp) ? MIN2(tmm + d3_2, tmm_2 + d5) : tmm + tmm_2;}
    00646     else if(q+2 == j){
    00647       energy += MIN2(tmm, MIN2(tmm_2, MIN2(d5 + d5_2, d3 + d3_2)));
    00648     }
    00649     else energy += MIN2(d3, d5_2);
    00650   }
    00651   else{
    00652     if(q+2 < j){ energy += d5 + d3_2;}
    00653     else if(q+2 == j){ energy += MIN2(d5, d3_2);}
    00654   }
    00655   return energy;
    00656 }
    00657
    00658 #endif