binaries/src/ViennaRNA/lib/circfold.inc

   1 /* -*-C-*- */
   2 /* this file contains code for folding circular RNAs */
   3 /* it's #include'd into fold.c */
   4
   5 PRIVATE void fill_arrays_circ(const char *string, int *bt){
   6   /* variant of fold() for circular RNAs */
   7   /* auxiliarry arrays:
   8      fM2 = multiloop region with exactly two stems, extending to 3' end
   9      for stupid dangles=1 case we also need:
  10      fM_d3 = multiloop region with >= 2 stems, starting at pos 2
  11              (a pair (k,n) will form 3' dangle with pos 1)
  12      fM_d5 = multiloop region with >= 2 stems, extending to pos n-1
  13              (a pair (1,k) will form a 5' dangle with pos n)
  14   */
  15   int Hi, Hj, Ii, Ij, Ip, Iq, Mi;
  16   int *fM_d3, *fM_d5, Md3i, Md5i, FcMd3, FcMd5;
  17   int i,j, p,q,length, energy;
  18   int dangle_model = P->model_details.dangles;
  19
  20   length = (int) strlen(string);
  21
  22   FcH = FcI= FcM = FcMd3= FcMd5= Fc = INF;
  23   for (i=1; i<length; i++)
  24     for (j=i+TURN+1; j <= length; j++) {
  25       int ij, bonus=0, type, u, new_c, no_close;
  26       u = length-j + i-1;
  27       if (u<TURN) continue;
  28
  29       ij = indx[j]+i;
  30       type = ptype[ij];
  31
  32       /* enforcing structure constraints */
  33       if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
  34       if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
  35       if ((BP[i]==-4)||(BP[j]==-4)) type=0;
  36
  37       no_close = (((type==3)||(type==4))&&no_closingGU&&(bonus==0));
  38
  39       /* if (j-i-1 > max_separation) type = 0; */  /* forces locality degree */
  40
  41       type=rtype[type];
  42       if (!type) continue;
  43       if (no_close) new_c = FORBIDDEN;
  44       else {
  45         char loopseq[10];
  46         /*int si1, sj1;*/
  47         if (u<7) {
  48           strcpy(loopseq , string+j-1);
  49           strncat(loopseq, string, i);
  50         }
  51         new_c = E_Hairpin(u, type, S1[j+1], S1[i-1],  loopseq, P)+bonus+c[ij];
  52       }
  53       if (new_c<FcH) {
  54         FcH = new_c; Hi=i; Hj=j;
  55       }
  56
  57       for (p = j+1; p < length ; p++) {
  58         int u1, qmin;
  59         u1 = p-j-1;
  60         if (u1+i-1>MAXLOOP) break;
  61         qmin = u1+i-1+length-MAXLOOP;
  62         if (qmin<p+TURN+1) qmin = p+TURN+1;
  63         for (q = qmin; q <=length; q++) {
  64           int u2, type_2/*, si1, sq1*/;
  65           type_2 = rtype[ptype[indx[q]+p]];
  66           if (type_2==0) continue;
  67           u2 = i-1 + length-q;
  68           if (u1+u2>MAXLOOP) continue;
  69           energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P);
  70           new_c = c[ij] + c[indx[q]+p] + energy;
  71           if (new_c<FcI) {
  72             FcI = new_c; Ii=i; Ij=j; Ip=p; Iq=q;
  73           }
  74         }
  75       }
  76     }
  77   Fc = MIN2(FcI, FcH);
  78
  79   /* compute the fM2 array (multi loops with exactly 2 helices) */
  80   /* to get a unique ML decomposition, just use fM1 instead of fML
  81      below. However, that will not work with dangle_model==1  */
  82   for (i=1; i<length-TURN; i++) {
  83     int u;
  84     fM2[i] = INF;
  85     for (u=i+TURN; u<length-TURN; u++)
  86       fM2[i] = MIN2(fM2[i], fML[indx[u]+i] + fML[indx[length]+u+1]);
  87   }
  88
  89   for (i=TURN+1; i<length-2*TURN; i++) {
  90     int fm;
  91     fm = fML[indx[i]+1]+fM2[i+1]+P->MLclosing;
  92     if (fm<FcM) {
  93       FcM=fm; Mi=i;
  94     }
  95   }
  96   Fc = MIN2(Fc, FcM);
  97
  98   if (dangle_model==1) {
  99     int u;
 100     fM_d3 =  (int *) space(sizeof(int)*(length+2));
 101     fM_d5 =  (int *) space(sizeof(int)*(length+2));
 102     for (i=TURN+1; i<length-TURN; i++) {
 103       fM_d3[i] = INF;
 104       for (u=2+TURN; u<i-TURN; u++)
 105         fM_d3[i] = MIN2(fM_d3[i], fML[indx[u]+2] + fML[indx[i]+u+1]);
 106     }
 107     for (i=2*TURN+1; i<length-TURN; i++) {
 108       int fm, type;
 109       type = ptype[indx[length]+i+1];
 110       if (type==0) continue;
 111       fm = fM_d3[i]+c[indx[length]+i+1]+E_MLstem(type, -1, S1[1], P) + P->MLclosing;
 112       if (fm<FcMd3) {
 113         FcMd3=fm; Md3i=i;
 114       }
 115       fm = fM_d3[i-1]+c[indx[length]+i+1]+E_MLstem(type, S1[i], S1[1], P) + P->MLclosing;
 116       if (fm<FcMd3) {
 117         FcMd3=fm; Md3i=-i;
 118       }
 119     }
 120
 121     for (i=TURN+1; i<length-TURN; i++) {
 122       fM_d5[i] = INF;
 123       for (u=i+TURN; u<length-TURN; u++)
 124         fM_d5[i] = MIN2(fM_d5[i], fML[indx[u]+i] + fML[indx[length-1]+u+1]);
 125     }
 126     for (i=TURN+1; i<length-2*TURN; i++) {
 127       int fm, type;
 128       type = ptype[indx[i]+1];
 129       if (type==0) continue;
 130       fm = E_MLstem(type, S1[length], -1, P) + c[indx[i]+1] + fM_d5[i+1] + P->MLclosing;
 131       if (fm<FcMd5) {
 132         FcMd5=fm; Md5i=i;
 133       }
 134       fm = E_MLstem(type, S1[length], S1[i+1], P) + c[indx[i]+1] + fM_d5[i+2] + P->MLclosing;
 135       if (fm<FcMd5) {
 136         FcMd5=fm; Md5i=-i;
 137       }
 138     }
 139     if (FcMd5<MIN2(Fc,FcMd3)) {
 140       /* looks like we have to do this ... */
 141       sector[++(*bt)].i = 1;
 142       sector[(*bt)].j = (Md5i>0)?Md5i:-Md5i;
 143       sector[(*bt)].ml = 2;
 144       i = (Md5i>0)?Md5i+1 : -Md5i+2; /* let's backtrack fm_d5[Md5i+1] */
 145       for (u=i+TURN; u<length-TURN; u++)
 146         if (fM_d5[i] == fML[indx[u]+i] + fML[indx[length-1]+u+1]) {
 147           sector[++(*bt)].i = i;
 148           sector[(*bt)].j = u;
 149           sector[(*bt)].ml = 1;
 150           sector[++(*bt)].i =u+1;
 151           sector[(*bt)].j = length-1;
 152           sector[(*bt)].ml = 1;
 153           break;
 154         }
 155       Fc = FcMd5;
 156     } else if (FcMd3<Fc) {
 157       /* here we go again... */
 158       sector[++(*bt)].i = (Md3i>0)?Md3i+1:-Md3i+1;
 159       sector[(*bt)].j = length;
 160       sector[(*bt)].ml = 2;
 161       i = (Md3i>0)? Md3i : -Md3i-1; /* let's backtrack fm_d3[Md3i] */
 162       for (u=2+TURN; u<i-TURN; u++)
 163         if (fM_d3[i] == fML[indx[u]+2] + fML[indx[i]+u+1]) {
 164           sector[++(*bt)].i = 2;
 165           sector[(*bt)].j = u;
 166           sector[(*bt)].ml = 1;
 167           sector[++(*bt)].i =u+1;
 168           sector[(*bt)].j = i;
 169           sector[(*bt)].ml = 1;
 170           break;
 171         }
 172       Fc = FcMd3;
 173     }
 174     free(fM_d3);
 175     free(fM_d5);
 176   }
 177   else if(Fc < INF){
 178     if (FcH==Fc) {
 179       sector[++(*bt)].i = Hi;
 180       sector[(*bt)].j = Hj;
 181       sector[(*bt)].ml = 2;
 182     }
 183     else if (FcI==Fc) {
 184       sector[++(*bt)].i = Ii;
 185       sector[(*bt)].j = Ij;
 186       sector[(*bt)].ml = 2;
 187       sector[++(*bt)].i = Ip;
 188       sector[(*bt)].j = Iq;
 189       sector[(*bt)].ml = 2;
 190     }
 191     else if (FcM==Fc) { /* grumpf we found a Multiloop */
 192       int fm, u;
 193       /* backtrack in fM2 */
 194       fm = fM2[Mi+1];
 195       for (u=Mi+TURN+1; u<length-TURN; u++)
 196         if (fm == fML[indx[u]+Mi+1] + fML[indx[length]+u+1]) {
 197                 sector[++(*bt)].i=Mi+1;
 198                 sector[(*bt)].j=u;
 199                 sector[(*bt)].ml = 1;
 200                 sector[++(*bt)].i=u+1;
 201                 sector[(*bt)].j=length;
 202                 sector[(*bt)].ml = 1;
 203                 break;
 204         }
 205       sector[++(*bt)].i = 1;
 206       sector[(*bt)].j = Mi;
 207       sector[(*bt)].ml = 1;
 208     }
 209   }
 210 }
 211