binaries/src/ViennaRNA/lib/plex_functions.c

   1 /* Last changed Time-stamp: <2010-06-30 17:24:43 wolfgang> */
   2 /*
   3            compute potentially pseudoknotted structures of an RNA sequence
   4                              Ivo Hofacker
   5                           Vienna RNA package
   6 */
   7
   8 /*
   9   library containing the function used in PKplex
  10   it generates pseudoknotted structures by letting the sequence form a duplex structure with itself
  11 */
  12
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16 #include <ctype.h>
  17 #include <string.h>
  18 #include <time.h>
  19
  20 #include "energy_par.h"
  21 #include "fold_vars.h"
  22 #include "utils.h"
  23 #include "params.h"
  24 #include "loop_energies.h"
  25
  26 #include "pair_mat.h"
  27
  28 #include "fold.h"
  29 #include "PKplex.h"
  30
  31 #undef  MAXLOOP
  32 #define MAXLOOP 10
  33
  34 PRIVATE void  update_dfold_params(void);
  35 PRIVATE void  duplexfold_XS(const char *s1, int **access_s1, const int threshold, const int max_interaction_length);
  36 PRIVATE char  *backtrack_XS(int kk, int ll, const int ii, const int jj, const int max_interaction_length);
  37 PRIVATE void  make_ptypes(const char *structure);
  38
  39 PRIVATE paramT  *P = NULL;
  40 PRIVATE int     ***c3 = NULL;      /* energy array used in duplexfold */
  41 PRIVATE short   *S1 = NULL, *SS1 = NULL;
  42 PRIVATE int     n1;
  43 PRIVATE char    *ptype = NULL;     /* precomputed array of pair types */
  44 PRIVATE int     *indx = NULL;      /* index for moving in the triangle matrices ptype[] */
  45
  46 PUBLIC  dupVar  *PlexHits = NULL;
  47 PUBLIC  int     PlexHitsArrayLength = 100;
  48 PUBLIC  int     NumberOfHits        = 0;
  49 PUBLIC  int     verbose             = 0;
  50
  51 /*-----------------------------------------------------------------------duplexfold_XS---------------------------------------------------------------------------*/
  52
  53 PRIVATE void  duplexfold_XS(const char *s1,
  54                             int **access_s1,
  55                             const int threshold,
  56                             const int max_interaction_length){
  57
  58   int i, j, k, l, p, q, Emin=INF, l_min=0, k_min=0, j_min=0;
  59   int type, type2, type3, E, tempK;
  60   char *struc;
  61   int length = (int) strlen(s1);
  62   struc=NULL;
  63
  64   c3 = (int ***) space(sizeof(int **) * (length));
  65   for (i=0; i<length; i++){
  66     c3[i] = (int **) space(sizeof(int *) * max_interaction_length);
  67     for (j=0; j<max_interaction_length; j++) {
  68       c3[i][j]=(int *) space(sizeof(int) * max_interaction_length);
  69     }
  70   }
  71
  72   i = length - 9;
  73
  74   while( i-- > 11 ){
  75     Emin=INF;
  76     j_min=0;
  77     l_min=0;
  78     k_min=0;
  79
  80     /* init all matrix elements to INF */
  81     for (j=0; j < length; j++){
  82       for(k=0;k<max_interaction_length;k++){
  83         for(l=0;l<max_interaction_length;l++){
  84           c3[j][k][l]=INF;
  85         }
  86       }
  87     }
  88     char string[10] = {'\0'};
  89     /* matrix starting values for (i,j)-basepairs */
  90     for(j=i+4; j<n1-10; j++) {
  91       type=ptype[indx[j]+i];
  92       if (type) {
  93         c3[j-11][max_interaction_length-1][0] = P->DuplexInit;
  94         c3[j-11][max_interaction_length-1][0] += E_Hairpin(j-i-1, type,  SS1[i+1], SS1[j-1], string, P);
  95 /*        c3[j-11][max_interaction_length-1][0] += E_ExtLoop(type, SS1[i+1], SS1[j-1], P); */
  96 /*           c3[j-11][max_interaction_length-1][0] += E_ExtLoop(rtype[type], SS1[j-1], SS1[i+1], P); */
  97        }
  98     }
  99
 100     int i_pos_begin=MAX2(9, i-max_interaction_length); /* why 9 ??? */
 101
 102     /* fill matrix */
 103     for(k=i-1; k>i_pos_begin; k--) {
 104       tempK=max_interaction_length-i+k-1;
 105       for(l = i + 5; l < n1 - 9; l++) { /* again, why 9 less then the sequence length ? */
 106         type2 = ptype[indx[l] + k];
 107         if(!type2) continue;
 108         for(p = k + 1; (p <= i) && (p <= k + MAXLOOP + 1); p++){
 109           for(q = l - 1; (q>=i+4) && (q >= l - MAXLOOP - 1); q--){
 110             if (p - k + l - q - 2 > MAXLOOP) break;
 111             type3 = ptype[indx[q] + p];
 112             if(!type3) continue;
 113             E = E_IntLoop(p - k - 1, l - q - 1, type2, rtype[type3], SS1[k + 1], SS1[l - 1], SS1[p - 1], SS1[q + 1], P);
 114             for(j = MAX2(i + 4, l - max_interaction_length + 1); j <= q; j++){
 115               type = ptype[indx[j]+i];
 116               if (type){
 117                 c3[j-11][tempK][l-j] = MIN2(c3[j-11][tempK][l-j], c3[j-11][max_interaction_length-i+p-1][q-j]+E);
 118               }
 119             }/* next j */
 120           }/* next q */
 121         }/* next p */
 122       }/* next l */
 123     }/* next k */
 124
 125     /* read out matrix minimum */
 126     for(j=i+4; j<n1-10; j++) {
 127       type=ptype[indx[j]+i];
 128       if (!type) continue;
 129       int j_pos_end=MIN2(n1-9,j+max_interaction_length);
 130       for (k=i-1; k>i_pos_begin; k--) {
 131         for (l=j+1; l<j_pos_end; l++) {
 132           type2=ptype[indx[l]+k];
 133           if (!type2) continue;
 134           E = c3[j-11][max_interaction_length-i+k-1][l-j];
 135 /*           printf("[%d,%d][%d,%d]\t%6.2f\t%6.2f\t%6.2f\n", i, k, l, j, E/100., access_s1[i-k+1][i]/100., access_s1[l-j+1][l]/100.); */
 136           E+=access_s1[i-k+1][i]+access_s1[l-j+1][l];
 137           E+=E_ExtLoop(type2,((k>i_pos_begin+1)? SS1[k-1]:-1),((l<j_pos_end-1)? SS1[l+1]:-1),P);
 138           E+=E_ExtLoop(rtype[type], SS1[j-1], SS1[i+1], P);
 139           if (E<Emin) {
 140             Emin=E; k_min=k; l_min=l;
 141             j_min=j;
 142           }
 143         }
 144       }
 145     }
 146
 147     if(Emin  < threshold){
 148       struc = backtrack_XS(k_min, l_min, i, j_min, max_interaction_length);
 149
 150       /* lets take care of the dangles */
 151       /* find best combination */
 152       int dx_5, dx_3, dy_5, dy_3,dGx,dGy,bonus_x, bonus_y;
 153       dx_5 = dx_3 = dy_5 = dy_3 = dGx = dGy = bonus_x = bonus_y = 0;
 154       dGx = access_s1[i-k_min+1][i];
 155       dGy = access_s1[l_min-j_min+1][l_min];
 156       PlexHits[NumberOfHits].tb=k_min -10 -dx_5;
 157       PlexHits[NumberOfHits].te=i -10 + dx_3;
 158       PlexHits[NumberOfHits].qb=j_min -10 - dy_5;
 159       PlexHits[NumberOfHits].qe=l_min -10 + dy_3;
 160       PlexHits[NumberOfHits].ddG=(double) Emin * 0.01;
 161       PlexHits[NumberOfHits].dG1=(double) dGx*0.01 ;
 162       PlexHits[NumberOfHits].dG2=(double) dGy*0.01 ;
 163       PlexHits[NumberOfHits].energy = PlexHits[NumberOfHits].ddG - PlexHits[NumberOfHits].dG1 - PlexHits[NumberOfHits].dG2;
 164       PlexHits[NumberOfHits].structure = struc;
 165
 166       /* output: */
 167       if(PlexHits[NumberOfHits].energy * 100 < threshold){
 168         if (verbose) printf("%s %3d,%-3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f)\n", PlexHits[NumberOfHits].structure, PlexHits[NumberOfHits].tb, PlexHits[NumberOfHits].te, PlexHits[NumberOfHits].qb, PlexHits[NumberOfHits].qe, PlexHits[NumberOfHits].ddG, PlexHits[NumberOfHits].energy, PlexHits[NumberOfHits].dG1, PlexHits[NumberOfHits].dG2);
 169         NumberOfHits++;
 170         if(NumberOfHits==PlexHitsArrayLength-1){
 171           PlexHitsArrayLength*=2;
 172           PlexHits = (dupVar *) xrealloc(PlexHits,sizeof(dupVar)*PlexHitsArrayLength);
 173         }
 174       }
 175     }
 176   }
 177
 178   for (i=0; i<(n1-20); i++) {
 179     for (j=0; j<max_interaction_length; j++) {
 180       free(c3[i][j]);
 181     }
 182     free(c3[i]);
 183   }
 184   free(c3);
 185 }
 186
 187
 188 PRIVATE char *backtrack_XS(int k, int l, const int i, const int j, const int max_interaction_length) {
 189   /* backtrack structure going backwards from i, and forwards from j
 190      return structure in bracket notation with & as separator */
 191   int p, q, type, type2, E, traced, i0, j0;
 192   char *st1, *st2, *struc;
 193   st1 = (char *) space(sizeof(char)*(i-k+2));
 194   st1[i-k+1]='\0';
 195   st2 = (char *) space(sizeof(char)*(l-j+2));
 196   st2[l-j+1]='\0';
 197
 198   i0=k; j0=l;
 199   while (k<=i && l>=j) {
 200     E = c3[j-11][max_interaction_length-i+k-1][l-j]; traced=0;
 201     st1[k-i0] = '(';
 202     st2[l-j] = ')';
 203
 204     type=ptype[indx[l]+k];
 205     if (!type) nrerror("backtrack failed in fold duplex bli");
 206     for (p=k+1; p<=i; p++) {
 207       for (q=l-1; q>=j; q--) {
 208         int LE;
 209         if (p-k+l-q-2>MAXLOOP) break;
 210         type2=ptype[indx[q]+p];
 211         if (!type2) continue;
 212          LE = E_IntLoop(p-k-1, l-q-1, type, rtype[type2], SS1[k+1], SS1[l-1], SS1[p-1], SS1[q+1], P);
 213          if (E == c3[j-11][max_interaction_length-i+p-1][q-j]+LE) {
 214           traced=1;
 215            k=p; l=q;
 216           break;
 217         }
 218       }
 219       if (traced) break;
 220     }
 221     if (!traced) {
 222       E-=E_ExtLoop(type2, ((k<i)?SS1[k+1]:-1), ((l>j-1)? SS1[l-1]:-1), P);
 223       break;
 224       if (E != P->DuplexInit) {
 225         nrerror("backtrack failed in fold duplex bal");
 226       } else break;
 227     }
 228   }
 229   struc = (char *) space(k-i0+1+j0-l+1+2);
 230
 231   for (p=0; p<=i-i0; p++){
 232     if (!st1[p]) st1[p] = '.';
 233   }
 234
 235   for (p=0; p<=j0-j; p++) {
 236     if (!st2[p]) {
 237       st2[p] = '.';
 238     }
 239   }
 240
 241   strcpy(struc, st1);
 242   strcat(struc, "&");
 243   strcat(struc, st2);
 244   free(st1); free(st2);
 245   return struc;
 246 }
 247
 248 PUBLIC  dupVar  **PKLduplexfold_XS( const char *s1,
 249                                     int **access_s1,
 250                                     const int threshold,
 251                                     const int max_interaction_length,
 252                                     const int delta){
 253
 254   if ((!P) || (fabs(P->temperature - temperature)>1e-6))
 255     update_dfold_params();
 256
 257   n1 = (int) strlen(s1);
 258   S1 = encode_sequence(s1, 0);
 259   SS1 = encode_sequence(s1, 1);
 260
 261   indx  = get_indx(n1);
 262   ptype = (char *) space(sizeof(char)*((n1*(n1+1))/2+2));
 263   make_ptypes(s1);
 264
 265   P->DuplexInit=0;
 266   duplexfold_XS(s1,access_s1,threshold, max_interaction_length);
 267   free(S1); free(SS1);
 268   free(indx); free(ptype);
 269   return NULL;
 270 }
 271
 272 /*---------------------------------UTILS------------------------------------------*/
 273
 274 PRIVATE void update_dfold_params(void){
 275   if(P) free(P);
 276   P = scale_parameters();
 277   make_pair_matrix();
 278 }
 279
 280 /*---------------------------------------------------------------------------*/
 281
 282 PRIVATE void make_ptypes(const char *structure) {
 283   int n,i,j,k,l;
 284
 285   n=S1[0];
 286   for (k=1; k<n-TURN; k++)
 287     for (l=1; l<=2; l++) {
 288       int type,ntype=0,otype=0;
 289       i=k; j = i+TURN+l; if (j>n) continue;
 290       type = pair[S1[i]][S1[j]];
 291       while ((i>=1)&&(j<=n)) {
 292         if ((i>1)&&(j<n)) ntype = pair[S1[i-1]][S1[j+1]];
 293         if (noLonelyPairs && (!otype) && (!ntype))
 294           type = 0; /* i.j can only form isolated pairs */
 295         ptype[indx[j]+i] = (char) type;
 296         otype =  type;
 297         type  = ntype;
 298         i--; j++;
 299       }
 300     }
 301 }