RNAlib-2.1.1: /mount/gensoft2/src/ViennaRNA/ViennaRNA-2.1.1/H/loop

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
/mount/gensoft2/src/ViennaRNA/ViennaRNA-2.1.1/H/loop_energies.h
