binaries/src/ViennaRNA/lib/alicircfold.inc

   1 /* -*-C-*- */
   2 /* this file contains code for alifolding circular RNAs */
   3 /* it's #include'd into alifold.c */
   4
   5 float circalifold(const char **strings, char *structure) {
   6   /* variant of alifold() 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 *fM2, Fc, FcH, FcI, FcM, Hi, Hj, Ii, Ij, Ip, Iq, Mi;
  16   int FcMd3, FcMd5;
  17   int i,j, p,q,length, energy;
  18   int n_seq, psc, s=0;
  19   int *type;
  20
  21   circular = 1;
  22   /* if (!uniq_ML) {uniq_ML=1; init_length=-1;} */
  23   length = (int) strlen(strings[0]);
  24
  25 #ifdef USE_OPENMP
  26   /* always init everything since all global static variables are uninitialized when entering a thread */
  27   init_alifold(length);
  28 #else
  29   if (length > init_length) init_alifold(length);
  30 #endif
  31   if (fabs(P->temperature - temperature)>1e-6) update_alifold_params();
  32
  33   oldAliEn=1; /* TODO allow new gap treatment in alicircfold() */
  34   for (s=0; strings[s]!=NULL; s++);
  35   n_seq = s;
  36   type = (int *) space(n_seq*sizeof(int));
  37
  38   alloc_sequence_arrays(strings, &S, &S5, &S3, &a2s, &Ss, circular);
  39   make_pscores((const short **) S, strings, n_seq, structure);
  40
  41   /* compute all arrays used by fold() for linear RNA */
  42   energy = fill_arrays((const char **) strings);
  43
  44   /* extra arrays for circfold() */
  45   fM2 =  (int *) space(sizeof(int)*(length+2));
  46   FcH = FcI= FcM = FcMd3= FcMd5= Fc = INF;
  47   for (i=1; i<length; i++)
  48     for (j=i+TURN+1; j <= length; j++) {
  49       int ij, u, new_c;
  50       u = length-j + i-1;
  51       if (u<TURN) continue;
  52
  53       ij = indx[j]+i;
  54
  55       psc = pscore[indx[j]+i];
  56
  57       if (psc<MINPSCORE) continue;
  58
  59       new_c = 0;
  60       for (s=0; s<n_seq; s++) {
  61         char loopseq[10];
  62         int rt, si1, sj1;
  63
  64         type[s] = pair[S[s][i]][S[s][j]];
  65         if (type[s]==0) type[s]=7;
  66
  67         rt=rtype[type[s]];
  68         if (u<7) {
  69           strcpy(loopseq , strings[s]+j-1);
  70           strncat(loopseq, strings[s], i);
  71         }
  72         si1 = (i==1)?S[s][length] : S[s][i-1];
  73         sj1 = (j==length)?S[s][1] : S[s][j+1];
  74         new_c += E_Hairpin(u, rt, sj1, si1,  loopseq, P);
  75       }
  76       new_c += +c[ij];
  77       if (new_c<FcH) {
  78         FcH = new_c; Hi=i; Hj=j;
  79       }
  80
  81       for (p = j+1; p < length ; p++) {
  82         int u1, qmin;
  83         u1 = p-j-1;
  84         if (u1+i-1>MAXLOOP) break;
  85         qmin = u1+i-1+length-MAXLOOP;
  86         if (qmin<p+TURN+1) qmin = p+TURN+1;
  87         for (q = qmin; q <=length; q++) {
  88           int u2, si1, sq1, type_2;
  89
  90           if (pscore[indx[q]+p]<MINPSCORE) continue;
  91
  92           u2 = i-1 + length-q;
  93           if (u1+u2>MAXLOOP) continue;
  94           for (energy = s=0; s<n_seq; s++) {
  95             int rt;
  96             rt=rtype[type[s]];
  97             type_2 = pair[S[s][q]][S[s][p]]; /* q,p not p,q! */
  98             if (type_2 ==0) type_2=7;
  99             si1 = (i==1)? S[s][length] : S[s][i-1];
 100             sq1 = (q==length)? S[s][1] : S[s][q+1];
 101             energy += E_IntLoop(u1, u2, rt, type_2,
 102                                  S[s][j+1], si1, S[s][p-1], sq1, P);
 103           }
 104           new_c = c[ij] + c[indx[q]+p] + energy;
 105           if (new_c<FcI) {
 106             FcI = new_c; Ii=i; Ij=j; Ip=p; Iq=q;
 107           }
 108         }
 109       }
 110     }
 111   Fc = MIN2(FcI, FcH);
 112
 113   /* compute the fM2 array (multi loops with exactly 2 helices) */
 114   /* to get a unique ML decomposition, just use fM1 instead of fML
 115      below. However, that will not work with dangles==1  */
 116   for (i=1; i<length-TURN; i++) {
 117     int u;
 118     fM2[i] = INF;
 119     for (u=i+TURN; u<length-TURN; u++)
 120       fM2[i] = MIN2(fM2[i], fML[indx[u]+i] + fML[indx[length]+u+1]);
 121   }
 122
 123   for (i=TURN+1; i<length-2*TURN; i++) {
 124     int fm;
 125     fm = fML[indx[i]+1]+fM2[i+1]+n_seq*P->MLclosing;
 126     if (fm<FcM) {
 127       FcM=fm; Mi=i;
 128     }
 129   }
 130   Fc = MIN2(Fc, FcM);
 131
 132   if (FcH==Fc) {
 133     sector[++s].i = Hi;
 134     sector[s].j = Hj;
 135     sector[s].ml = 2;
 136   }
 137   else if (FcI==Fc) {
 138     sector[++s].i = Ii;
 139     sector[s].j = Ij;
 140     sector[s].ml = 2;
 141     sector[++s].i = Ip;
 142     sector[s].j = Iq;
 143     sector[s].ml = 2;
 144   }
 145   else if (FcM==Fc) { /* grumpf we found a Multiloop */
 146     int fm, u;
 147     /* backtrack in fM2 */
 148     fm = fM2[Mi+1];
 149     for (u=Mi+TURN+1; u<length-TURN; u++)
 150       if (fm == fML[indx[u]+Mi+1] + fML[indx[length]+u+1]) {
 151         sector[++s].i=Mi+1;
 152         sector[s].j=u;
 153         sector[s].ml = 1;
 154         sector[++s].i=u+1;
 155         sector[s].j=length;
 156         sector[s].ml = 1;
 157         break;
 158       }
 159     sector[++s].i = 1;
 160     sector[s].j = Mi;
 161     sector[s].ml = 1;
 162   }
 163   backtrack(strings, s);
 164
 165 #ifdef PAREN
 166   parenthesis_structure(structure, base_pair2, length);
 167 #else
 168   letter_structure(structure, base_pair2, length);
 169 #endif
 170
 171   /*
 172   *  Backward compatibility:
 173   *  This block may be removed if deprecated functions
 174   *  relying on the global variable "base_pair" vanishs from within the package!
 175   */
 176   base_pair = base_pair2;
 177   /*
 178   {
 179     if(base_pair) free(base_pair);
 180     base_pair = (bondT *)space(sizeof(bondT) * (1+length/2));
 181     memcpy(base_pair, base_pair2, sizeof(bondT) * (1+length/2));
 182   }
 183   */
 184   free(fM2);
 185   free(type);
 186   free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);
 187
 188   return (float) Fc/(n_seq*100.);
 189 }