binaries/src/ViennaRNA/lib/alifold.c

   1 /* Last changed Time-stamp: <2009-02-24 15:17:17 ivo> */
   2 /*
   3                   minimum free energy folding
   4                   for a set of aligned sequences
   5
   6                   c Ivo Hofacker
   7
   8                   Vienna RNA package
   9 */
  10
  11 /**
  12 *** \file alifold.c
  13 **/
  14
  15
  16 #include <config.h>
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <math.h>
  20 #include <ctype.h>
  21 #include <string.h>
  22 #include <limits.h>
  23
  24 #include "fold.h"
  25 #include "utils.h"
  26 #include "energy_par.h"
  27 #include "fold_vars.h"
  28 #include "pair_mat.h"
  29 #include "params.h"
  30 #include "ribo.h"
  31 #include "aln_util.h"
  32 #include "loop_energies.h"
  33 #include "gquad.h"
  34 #include "alifold.h"
  35
  36 #ifdef _OPENMP
  37 #include <omp.h>
  38 #endif
  39
  40 #define PAREN
  41 #define STACK_BULGE1  1     /* stacking energies for bulges of size 1 */
  42 #define NEW_NINIO     1     /* new asymetry penalty */
  43 #define MAXSECTORS    500   /* dimension for a backtrack array */
  44 #define LOCALITY        0.  /* locality parameter for base-pairs */
  45 #define UNIT 100
  46 #define MINPSCORE -2 * UNIT
  47
  48 /*
  49 #################################
  50 # GLOBAL VARIABLES              #
  51 #################################
  52 */
  53 PUBLIC double cv_fact=1.; /* should be made static to not interfere with other threads */
  54 PUBLIC double nc_fact=1.; /* should be made static to not interfere with other threads */
  55
  56 /*
  57 #################################
  58 # PRIVATE VARIABLES             #
  59 #################################
  60 */
  61 PRIVATE short           **S     = NULL;
  62 PRIVATE short           **S5    = NULL;     /*S5[s][i] holds next base 5' of i in sequence s*/
  63 PRIVATE short           **S3    = NULL;     /*Sl[s][i] holds next base 3' of i in sequence s*/
  64 PRIVATE char            **Ss    = NULL;
  65 PRIVATE unsigned short  **a2s   = NULL;
  66 PRIVATE paramT          *P      = NULL;
  67 PRIVATE int             *indx   = NULL;     /* index for moving in the triangle matrices c[] and fMl[]*/
  68 PRIVATE int             *c      = NULL;     /* energy array, given that i-j pair */
  69 PRIVATE int             *cc     = NULL;     /* linear array for calculating canonical structures */
  70 PRIVATE int             *cc1    = NULL;     /*   "     "        */
  71 PRIVATE int             *f5     = NULL;     /* energy of 5' end */
  72 PRIVATE int             *fML    = NULL;     /* multi-loop auxiliary energy array */
  73 PRIVATE int             *Fmi    = NULL;     /* holds row i of fML (avoids jumps in memory) */
  74 PRIVATE int             *DMLi   = NULL;     /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j])  */
  75 PRIVATE int             *DMLi1  = NULL;     /*             MIN(fML[i+1,k]+fML[k+1,j])  */
  76 PRIVATE int             *DMLi2  = NULL;     /*             MIN(fML[i+2,k]+fML[k+1,j])  */
  77 PRIVATE int             *pscore = NULL;     /* precomputed array of pair types */
  78 PRIVATE int             init_length = -1;
  79 PRIVATE sect            sector[MAXSECTORS]; /* stack of partial structures for backtracking */
  80 PRIVATE bondT           *base_pair2 = NULL;
  81 PRIVATE int             circular    = 0;
  82 PRIVATE int             with_gquad  = 0;
  83
  84 PRIVATE int             *ggg        = NULL; /* minimum free energies of the gquadruplexes */
  85 PRIVATE char            *cons_seq   = NULL;
  86 PRIVATE short           *S_cons     = NULL;
  87
  88 #ifdef _OPENMP
  89
  90 #pragma omp threadprivate(S, S5, S3, Ss, a2s, P, indx, c, cc, cc1, f5, fML, Fmi, DMLi, DMLi1, DMLi2,\
  91                           pscore, init_length, sector, base_pair2,\
  92                           ggg, with_gquad, cons_seq, S_cons)
  93
  94 #endif
  95
  96 /*
  97 #################################
  98 # PRIVATE FUNCTION DECLARATIONS #
  99 #################################
 100 */
 101
 102 PRIVATE void  init_alifold(int length);
 103 PRIVATE void  get_arrays(unsigned int size);
 104 PRIVATE void  make_pscores(const short *const *S, const char **AS, int n_seq, const char *structure);
 105 PRIVATE int   fill_arrays(const char **strings);
 106 PRIVATE void  backtrack(const char **strings, int s);
 107
 108 PRIVATE void    energy_of_alistruct_pt(const char **sequences,short * ptable, int n_seq, int *energy);
 109 PRIVATE void    stack_energy_pt(int i, const char **sequences, short *ptable,  int n_seq, int *energy);
 110 PRIVATE int     ML_Energy_pt(int i, int n_seq, short *pt);
 111 PRIVATE int     EL_Energy_pt(int i, int n_seq, short *pt);
 112
 113 PRIVATE void    en_corr_of_loop_gquad(int i,
 114                                       int j,
 115                                       const char **sequences,
 116                                       const char *structure,
 117                                       short *pt,
 118                                       int *loop_idx,
 119                                       int n_seq,
 120                                       int en[2]);
 121
 122 /*
 123 #################################
 124 # BEGIN OF FUNCTION DEFINITIONS #
 125 #################################
 126 */
 127
 128 /* unsafe function that will be replaced by a threadsafe companion in the future */
 129 PRIVATE void init_alifold(int length){
 130
 131 #ifdef _OPENMP
 132 /* Explicitly turn off dynamic threads */
 133   omp_set_dynamic(0);
 134 #endif
 135
 136   if (length < 1) nrerror("initialize_fold: argument must be greater 0");
 137   free_alifold_arrays();
 138   get_arrays((unsigned) length);
 139   init_length = length;
 140
 141   indx = get_indx((unsigned)length);
 142
 143   update_alifold_params();
 144 }
 145
 146 PRIVATE void get_arrays(unsigned int size){
 147   if(size >= (unsigned int)sqrt((double)INT_MAX))
 148     nrerror("get_arrays@alifold.c: sequence length exceeds addressable range");
 149
 150   c       = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 151   fML     = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 152   pscore  = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 153   f5      = (int *) space(sizeof(int)*(size+2));
 154   cc      = (int *) space(sizeof(int)*(size+2));
 155   cc1     = (int *) space(sizeof(int)*(size+2));
 156   Fmi     = (int *) space(sizeof(int)*(size+1));
 157   DMLi    = (int *) space(sizeof(int)*(size+1));
 158   DMLi1   = (int *) space(sizeof(int)*(size+1));
 159   DMLi2   = (int *) space(sizeof(int)*(size+1));
 160   if(base_pair2) free(base_pair2);
 161
 162   base_pair2 = (bondT *) space(sizeof(bondT)*(1+size/2+200)); /* add a guess of how many G's may be involved in a G quadruplex */
 163 }
 164
 165 PUBLIC  void  free_alifold_arrays(void){
 166   if(indx)        free(indx);
 167   if(c)           free(c);
 168   if(fML)         free(fML);
 169   if(f5)          free(f5);
 170   if(cc)          free(cc);
 171   if(cc1)         free(cc1);
 172   if(pscore)      free(pscore);
 173   if(base_pair2)  free(base_pair2);
 174   if(Fmi)         free(Fmi);
 175   if(DMLi)        free(DMLi);
 176   if(DMLi1)       free(DMLi1);
 177   if(DMLi2)       free(DMLi2);
 178   if(P)           free(P);
 179   if(ggg)         free(ggg);
 180   if(cons_seq)    free(cons_seq);
 181   if(S_cons)      free(S_cons);
 182
 183   indx = c = fML = f5 = cc = cc1 = Fmi = DMLi = DMLi1 = DMLi2 = ggg = NULL;
 184   pscore      = NULL;
 185   base_pair   = NULL;
 186   base_pair2  = NULL;
 187   P           = NULL;
 188   init_length = 0;
 189   cons_seq    = NULL;
 190   S_cons      = NULL;
 191 }
 192
 193
 194 PUBLIC void alloc_sequence_arrays(const char **sequences, short ***S, short ***S5, short ***S3, unsigned short ***a2s, char ***Ss, int circ){
 195   unsigned int s, n_seq, length;
 196   if(sequences[0] != NULL){
 197     length = strlen(sequences[0]);
 198     for (s=0; sequences[s] != NULL; s++);
 199     n_seq = s;
 200     *S    = (short **)          space((n_seq+1) * sizeof(short *));
 201     *S5   = (short **)          space((n_seq+1) * sizeof(short *));
 202     *S3   = (short **)          space((n_seq+1) * sizeof(short *));
 203     *a2s  = (unsigned short **) space((n_seq+1) * sizeof(unsigned short *));
 204     *Ss   = (char **)           space((n_seq+1) * sizeof(char *));
 205     for (s=0; s<n_seq; s++) {
 206       if(strlen(sequences[s]) != length) nrerror("uneqal seqence lengths");
 207       (*S5)[s]  = (short *)         space((length + 2) * sizeof(short));
 208       (*S3)[s]  = (short *)         space((length + 2) * sizeof(short));
 209       (*a2s)[s] = (unsigned short *)space((length + 2) * sizeof(unsigned short));
 210       (*Ss)[s]  = (char *)          space((length + 2) * sizeof(char));
 211       (*S)[s]   = (short *)         space((length + 2) * sizeof(short));
 212       encode_ali_sequence(sequences[s], (*S)[s], (*S5)[s], (*S3)[s], (*Ss)[s], (*a2s)[s], circ);
 213     }
 214     (*S5)[n_seq]  = NULL;
 215     (*S3)[n_seq]  = NULL;
 216     (*a2s)[n_seq] = NULL;
 217     (*Ss)[n_seq]  = NULL;
 218     (*S)[n_seq]   = NULL;
 219   }
 220   else nrerror("alloc_sequence_arrays: no sequences in the alignment!");
 221 }
 222
 223 PUBLIC void free_sequence_arrays(unsigned int n_seq, short ***S, short ***S5, short ***S3, unsigned short ***a2s, char ***Ss){
 224   unsigned int s;
 225   for (s=0; s<n_seq; s++) {
 226     free((*S)[s]);
 227     free((*S5)[s]);
 228     free((*S3)[s]);
 229     free((*a2s)[s]);
 230     free((*Ss)[s]);
 231   }
 232   free(*S);   *S    = NULL;
 233   free(*S5);  *S5   = NULL;
 234   free(*S3);  *S3   = NULL;
 235   free(*a2s); *a2s  = NULL;
 236   free(*Ss);  *Ss   = NULL;
 237 }
 238
 239 PUBLIC void update_alifold_params(void){
 240   if(P) free(P);
 241   P = scale_parameters();
 242   make_pair_matrix();
 243   if (init_length < 0) init_length=0;
 244 }
 245
 246 PUBLIC float alifold(const char **strings, char *structure){
 247   int  length, energy, s, n_seq;
 248
 249   circular = 0;
 250   length = (int) strlen(strings[0]);
 251
 252 #ifdef _OPENMP
 253   /* always init everything since all global static variables are uninitialized when entering a thread */
 254   init_alifold(length);
 255 #else
 256   if (length>init_length) init_alifold(length);
 257 #endif
 258   if (fabs(P->temperature - temperature)>1e-6)  update_alifold_params();
 259
 260   for (s=0; strings[s]!=NULL; s++);
 261   n_seq       = s;
 262   with_gquad  = P->model_details.gquad;
 263
 264   alloc_sequence_arrays(strings, &S, &S5, &S3, &a2s, &Ss, circular);
 265   make_pscores((const short **) S, strings, n_seq, structure);
 266
 267   energy = fill_arrays((const char **)strings);
 268
 269   backtrack((const char **)strings, 0);
 270
 271 #ifdef PAREN
 272   parenthesis_structure(structure, base_pair2, length);
 273 #else
 274   letter_structure(structure, base_pair2, length);
 275 #endif
 276
 277   /*
 278   *  Backward compatibility:
 279   *  This block may be removed if deprecated functions
 280   *  relying on the global variable "base_pair" vanishs from within the package!
 281   */
 282   base_pair = base_pair2;
 283   /*
 284   {
 285     if(base_pair) free(base_pair);
 286     base_pair = (bondT *)space(sizeof(bondT) * (1+length/2));
 287     memcpy(base_pair, base_pair2, sizeof(bondT) * (1+length/2));
 288   }
 289   */
 290   free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);
 291
 292   if (backtrack_type=='C')
 293     return (float) c[indx[length]+1]/(n_seq*100.);
 294   else if (backtrack_type=='M')
 295     return (float) fML[indx[length]+1]/(n_seq*100.);
 296   else
 297     return (float) f5[length]/(n_seq*100.);
 298 }
 299
 300
 301 /**
 302 *** the actual forward recursion to fill the energy arrays
 303 **/
 304 PRIVATE int fill_arrays(const char **strings) {
 305   int   i, j, k, p, q, length, energy, new_c;
 306   int   decomp, MLenergy, new_fML;
 307   int   s, n_seq, *type, type_2, tt;
 308
 309   /* count number of sequences */
 310   for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
 311
 312   type = (int *) space(n_seq*sizeof(int));
 313   length = strlen(strings[0]);
 314
 315   /* init energies */
 316
 317   if(with_gquad){
 318     cons_seq = consensus(strings);
 319     /* make g-island annotation of the consensus */
 320     S_cons = encode_sequence(cons_seq, 0);
 321     ggg = get_gquad_ali_matrix(S_cons, S, n_seq, P);
 322   }
 323
 324   for (j=1; j<=length; j++){
 325     Fmi[j]=DMLi[j]=DMLi1[j]=DMLi2[j]=INF;
 326     for (i=(j>TURN?(j-TURN):1); i<j; i++) {
 327       c[indx[j]+i] = fML[indx[j]+i] = INF;
 328     }
 329   }
 330
 331   /* begin recursions */
 332   for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
 333     for (j = i+TURN+1; j <= length; j++) {
 334       int ij, psc, l1, maxq, minq, up, c0;
 335       ij = indx[j]+i;
 336
 337       for (s=0; s<n_seq; s++) {
 338         type[s] = pair[S[s][i]][S[s][j]];
 339         if (type[s]==0) type[s]=7;
 340       }
 341
 342       psc = pscore[indx[j]+i];
 343       if (psc>=MINPSCORE) {   /* a pair to consider */
 344         int stackEnergy = INF;
 345         /* hairpin ----------------------------------------------*/
 346
 347
 348         for (new_c=s=0; s<n_seq; s++) {
 349           if ((a2s[s][j-1]-a2s[s][i])<3) new_c+=600;
 350           else  new_c += E_Hairpin(a2s[s][j-1]-a2s[s][i],type[s],S3[s][i],S5[s][j],Ss[s]+(a2s[s][i-1]), P);
 351        }
 352         /*--------------------------------------------------------
 353           check for elementary structures involving more than one
 354           closing pair.
 355           --------------------------------------------------------*/
 356
 357         for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
 358           minq = j-i+p-MAXLOOP-2;
 359           if (minq<p+1+TURN) minq = p+1+TURN;
 360           for (q = minq; q < j; q++) {
 361             if (pscore[indx[q]+p]<MINPSCORE) continue;
 362             for (energy = s=0; s<n_seq; s++) {
 363               type_2 = pair[S[s][q]][S[s][p]]; /* q,p not p,q! */
 364               if (type_2 == 0) type_2 = 7;
 365               energy += E_IntLoop(a2s[s][p-1]-a2s[s][i], a2s[s][j-1]-a2s[s][q], type[s], type_2,
 366                                    S3[s][i], S5[s][j],
 367                                    S5[s][p], S3[s][q], P);
 368             }
 369             new_c = MIN2(new_c, energy + c[indx[q]+p]);
 370             if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
 371           } /* end q-loop */
 372         } /* end p-loop */
 373
 374         /* multi-loop decomposition ------------------------*/
 375         decomp = DMLi1[j-1];
 376         if(dangles){
 377           for(s=0; s<n_seq; s++){
 378             tt = rtype[type[s]];
 379             decomp += E_MLstem(tt, S5[s][j], S3[s][i], P);
 380           }
 381         }
 382         else{
 383           for(s=0; s<n_seq; s++){
 384             tt = rtype[type[s]];
 385             decomp += E_MLstem(tt, -1, -1, P);
 386           }
 387         }
 388         MLenergy = decomp + n_seq*P->MLclosing;
 389         new_c = MIN2(new_c, MLenergy);
 390
 391         if(with_gquad){
 392           decomp = 0;
 393           for(s=0;s<n_seq;s++){
 394             tt = type[s];
 395             if(dangles == 2)
 396               decomp += P->mismatchI[tt][S3[s][i]][S5[s][j]];
 397             if(tt > 2)
 398               decomp += P->TerminalAU;
 399           }
 400           for(p = i + 2; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
 401             l1    = p - i - 1;
 402             if(l1>MAXLOOP) break;
 403             if(S_cons[p] != 3) continue;
 404             minq  = j - i + p - MAXLOOP - 2;
 405             c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 406             minq  = MAX2(c0, minq);
 407             c0    = j - 1;
 408             maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 409             maxq  = MIN2(c0, maxq);
 410             for(q = minq; q < maxq; q++){
 411               if(S_cons[q] != 3) continue;
 412               c0    = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[l1 + j - q - 1];
 413               new_c = MIN2(new_c, c0);
 414             }
 415           }
 416
 417           p = i + 1;
 418           if(S_cons[p] == 3){
 419             if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
 420               minq  = j - i + p - MAXLOOP - 2;
 421               c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 422               minq  = MAX2(c0, minq);
 423               c0    = j - 3;
 424               maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 425               maxq  = MIN2(c0, maxq);
 426               for(q = minq; q < maxq; q++){
 427                 if(S_cons[q] != 3) continue;
 428                 c0  = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[j - q - 1];
 429                 new_c   = MIN2(new_c, c0);
 430               }
 431             }
 432           }
 433           q = j - 1;
 434           if(S_cons[q] == 3)
 435             for(p = i + 4; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
 436               l1    = p - i - 1;
 437               if(l1>MAXLOOP) break;
 438               if(S_cons[p] != 3) continue;
 439               c0  = decomp + ggg[indx[q] + p] + n_seq * P->internal_loop[l1];
 440               new_c   = MIN2(new_c, c0);
 441             }
 442         }
 443
 444         new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
 445
 446         cc[j] = new_c - psc; /* add covariance bonnus/penalty */
 447         if (noLonelyPairs)
 448           c[ij] = cc1[j-1]+stackEnergy-psc;
 449         else
 450           c[ij] = cc[j];
 451
 452       } /* end >> if (pair) << */
 453
 454       else c[ij] = INF;
 455       /* done with c[i,j], now compute fML[i,j] */
 456       /* free ends ? -----------------------------------------*/
 457
 458       new_fML = fML[ij+1]+n_seq*P->MLbase;
 459       new_fML = MIN2(fML[indx[j-1]+i]+n_seq*P->MLbase, new_fML);
 460       energy = c[ij];
 461       if(dangles){
 462         for (s=0; s<n_seq; s++) {
 463           energy += E_MLstem(type[s], S5[s][i], S3[s][j], P);
 464         }
 465       }
 466       else{
 467         for (s=0; s<n_seq; s++) {
 468           energy += E_MLstem(type[s], -1, -1, P);
 469         }
 470       }
 471       new_fML = MIN2(energy, new_fML);
 472
 473       if(with_gquad){
 474         decomp = ggg[indx[j] + i] + n_seq * E_MLstem(0, -1, -1, P);
 475         new_fML = MIN2(new_fML, decomp);
 476       }
 477
 478       /* modular decomposition -------------------------------*/
 479       for (decomp = INF, k = i+1+TURN; k <= j-2-TURN; k++)
 480         decomp = MIN2(decomp, Fmi[k]+fML[indx[j]+k+1]);
 481
 482       DMLi[j] = decomp;               /* store for use in ML decompositon */
 483       new_fML = MIN2(new_fML,decomp);
 484
 485       /* coaxial stacking deleted */
 486
 487       fML[ij] = Fmi[j] = new_fML;     /* substring energy */
 488     } /* END for j */
 489
 490     {
 491       int *FF; /* rotate the auxilliary arrays */
 492       FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
 493       FF = cc1; cc1=cc; cc=FF;
 494       for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
 495     }
 496   } /* END for i */
 497   /* calculate energies of 5' and 3' fragments */
 498
 499   f5[TURN + 1] = 0;
 500   switch(dangles){
 501     case 0:   for(j = TURN + 2; j <= length; j++){
 502                 f5[j] = f5[j-1];
 503                 if (c[indx[j]+1]<INF){
 504                   energy = c[indx[j]+1];
 505                   for(s = 0; s < n_seq; s++){
 506                     tt = pair[S[s][1]][S[s][j]];
 507                     if(tt==0) tt=7;
 508                     energy += E_ExtLoop(tt, -1, -1, P);
 509                   }
 510                   f5[j] = MIN2(f5[j], energy);
 511                 }
 512
 513                 if(with_gquad){
 514                   if(ggg[indx[j]+1] < INF)
 515                     f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
 516                 }
 517
 518                 for(i = j - TURN - 1; i > 1; i--){
 519                   if(c[indx[j]+i]<INF){
 520                     energy = f5[i-1] + c[indx[j]+i];
 521                     for(s = 0; s < n_seq; s++){
 522                       tt = pair[S[s][i]][S[s][j]];
 523                       if(tt==0) tt=7;
 524                       energy += E_ExtLoop(tt, -1, -1, P);
 525                     }
 526                     f5[j] = MIN2(f5[j], energy);
 527                   }
 528
 529                   if(with_gquad){
 530                     if(ggg[indx[j]+i] < INF)
 531                       f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
 532                   }
 533
 534                 }
 535               }
 536               break;
 537     default:  for(j = TURN + 2; j <= length; j++){
 538                 f5[j] = f5[j-1];
 539                 if (c[indx[j]+1]<INF) {
 540                   energy = c[indx[j]+1];
 541                   for(s = 0; s < n_seq; s++){
 542                     tt = pair[S[s][1]][S[s][j]];
 543                     if(tt==0) tt=7;
 544                     energy += E_ExtLoop(tt, -1, (j<length) ? S3[s][j] : -1, P);
 545                   }
 546                   f5[j] = MIN2(f5[j], energy);
 547                 }
 548
 549                 if(with_gquad){
 550                   if(ggg[indx[j]+1] < INF)
 551                     f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
 552                 }
 553
 554                 for(i = j - TURN - 1; i > 1; i--){
 555                   if (c[indx[j]+i]<INF) {
 556                     energy = f5[i-1] + c[indx[j]+i];
 557                     for(s = 0; s < n_seq; s++){
 558                       tt = pair[S[s][i]][S[s][j]];
 559                       if(tt==0) tt=7;
 560                       energy += E_ExtLoop(tt, S5[s][i], (j < length) ? S3[s][j] : -1, P);
 561                     }
 562                     f5[j] = MIN2(f5[j], energy);
 563                   }
 564
 565                   if(with_gquad){
 566                     if(ggg[indx[j]+i] < INF)
 567                       f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
 568                   }
 569
 570                 }
 571               }
 572               break;
 573   }
 574   free(type);
 575   return(f5[length]);
 576 }
 577
 578 #include "alicircfold.inc"
 579
 580 /**
 581 *** backtrack in the energy matrices to obtain a structure with MFE
 582 **/
 583 PRIVATE void backtrack(const char **strings, int s) {
 584   /*------------------------------------------------------------------
 585     trace back through the "c", "f5" and "fML" arrays to get the
 586     base pairing list. No search for equivalent structures is done.
 587     This inverts the folding procedure, hence it's very fast.
 588     ------------------------------------------------------------------*/
 589    /* normally s=0.
 590      If s>0 then s items have been already pushed onto the sector stack */
 591   int   i, j, k, p, q, length, energy, up, c0, l1, minq, maxq;
 592   int   type_2, tt, mm;
 593   int   b=0, cov_en = 0;
 594   int   n_seq;
 595   int *type;
 596   length = strlen(strings[0]);
 597   for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
 598   type = (int *) space(n_seq*sizeof(int));
 599
 600   if (s==0) {
 601     sector[++s].i = 1;
 602     sector[s].j = length;
 603     sector[s].ml = (backtrack_type=='M') ? 1 : ((backtrack_type=='C')?2:0);
 604   }
 605   while (s>0) {
 606     int ss, ml, fij, fi, cij, traced, i1, j1, d3, d5, jj=0, gq=0;
 607     int canonical = 1;     /* (i,j) closes a canonical structure */
 608     i  = sector[s].i;
 609     j  = sector[s].j;
 610     ml = sector[s--].ml;   /* ml is a flag indicating if backtracking is to
 611                               occur in the fML- (1) or in the f-array (0) */
 612     if (ml==2) {
 613       base_pair2[++b].i = i;
 614       base_pair2[b].j   = j;
 615       cov_en += pscore[indx[j]+i];
 616       goto repeat1;
 617     }
 618
 619     if (j < i+TURN+1) continue; /* no more pairs in this interval */
 620
 621     fij = (ml)? fML[indx[j]+i] : f5[j];
 622     fi  = (ml)?(fML[indx[j-1]+i]+n_seq*P->MLbase):f5[j-1];
 623
 624     if (fij == fi) {  /* 3' end is unpaired */
 625       sector[++s].i = i;
 626       sector[s].j   = j-1;
 627       sector[s].ml  = ml;
 628       continue;
 629     }
 630
 631     if (ml == 0) { /* backtrack in f5 */
 632       switch(dangles){
 633         case 0:   /* j or j-1 is paired. Find pairing partner */
 634                   for (i=j-TURN-1,traced=0; i>=1; i--) {
 635                     int cc, en;
 636                     jj = i-1;
 637                     if (c[indx[j]+i]<INF) {
 638                       en = c[indx[j]+i] + f5[i-1];
 639                       for(ss = 0; ss < n_seq; ss++){
 640                         type[ss] = pair[S[ss][i]][S[ss][j]];
 641                         if (type[ss]==0) type[ss] = 7;
 642                         en += E_ExtLoop(type[ss], -1, -1, P);
 643                       }
 644                       if (fij == en) traced=j;
 645                     }
 646
 647                     if(with_gquad){
 648                       if(fij == f5[i-1] + ggg[indx[j]+i]){
 649                         /* found the decomposition */
 650                         traced = j; jj = i - 1; gq = 1;
 651                         break;
 652                       }
 653                     }
 654
 655                     if (traced) break;
 656                   }
 657                   break;
 658         default:  /* j or j-1 is paired. Find pairing partner */
 659                   for (i=j-TURN-1,traced=0; i>=1; i--) {
 660                     int cc, en;
 661                     jj = i-1;
 662                     if (c[indx[j]+i]<INF) {
 663                       en = c[indx[j]+i] + f5[i-1];
 664                       for(ss = 0; ss < n_seq; ss++){
 665                         type[ss] = pair[S[ss][i]][S[ss][j]];
 666                         if (type[ss]==0) type[ss] = 7;
 667                         en += E_ExtLoop(type[ss], (i>1) ? S5[ss][i]: -1, (j < length) ? S3[ss][j] : -1, P);
 668                       }
 669                       if (fij == en) traced=j;
 670                     }
 671
 672                     if(with_gquad){
 673                       if(fij == f5[i-1] + ggg[indx[j]+i]){
 674                         /* found the decomposition */
 675                         traced = j; jj = i - 1; gq = 1;
 676                         break;
 677                       }
 678                     }
 679
 680                     if (traced) break;
 681                   }
 682                   break;
 683       }
 684
 685       if (!traced) nrerror("backtrack failed in f5");
 686       /* push back the remaining f5 portion */
 687       sector[++s].i = 1;
 688       sector[s].j   = jj;
 689       sector[s].ml  = ml;
 690
 691       /* trace back the base pair found */
 692       j=traced;
 693
 694       if(with_gquad && gq){
 695         /* goto backtrace of gquadruplex */
 696         goto repeat_gquad;
 697       }
 698
 699       base_pair2[++b].i = i;
 700       base_pair2[b].j   = j;
 701       cov_en += pscore[indx[j]+i];
 702       goto repeat1;
 703     }
 704     else { /* trace back in fML array */
 705       if (fML[indx[j]+i+1]+n_seq*P->MLbase == fij) { /* 5' end is unpaired */
 706         sector[++s].i = i+1;
 707         sector[s].j   = j;
 708         sector[s].ml  = ml;
 709         continue;
 710       }
 711
 712       if(with_gquad){
 713         if(fij == ggg[indx[j]+i] + n_seq * E_MLstem(0, -1, -1, P)){
 714           /* go to backtracing of quadruplex */
 715           goto repeat_gquad;
 716         }
 717       }
 718
 719       cij = c[indx[j]+i];
 720       if(dangles){
 721         for(ss = 0; ss < n_seq; ss++){
 722           tt = pair[S[ss][i]][S[ss][j]];
 723           if(tt==0) tt=7;
 724           cij += E_MLstem(tt, S5[ss][i], S3[ss][j], P);
 725         }
 726       }
 727       else{
 728         for(ss = 0; ss < n_seq; ss++){
 729           tt = pair[S[ss][i]][S[ss][j]];
 730           if(tt==0) tt=7;
 731           cij += E_MLstem(tt, -1, -1, P);
 732         }
 733       }
 734
 735       if (fij==cij){
 736         /* found a pair */
 737         base_pair2[++b].i = i;
 738         base_pair2[b].j   = j;
 739         cov_en += pscore[indx[j]+i];
 740         goto repeat1;
 741       }
 742
 743       for (k = i+1+TURN; k <= j-2-TURN; k++)
 744         if (fij == (fML[indx[k]+i]+fML[indx[j]+k+1]))
 745           break;
 746
 747       sector[++s].i = i;
 748       sector[s].j   = k;
 749       sector[s].ml  = ml;
 750       sector[++s].i = k+1;
 751       sector[s].j   = j;
 752       sector[s].ml  = ml;
 753
 754       if (k>j-2-TURN) nrerror("backtrack failed in fML");
 755       continue;
 756     }
 757
 758   repeat1:
 759
 760     /*----- begin of "repeat:" -----*/
 761     if (canonical)  cij = c[indx[j]+i];
 762
 763     for (ss=0; ss<n_seq; ss++) {
 764       type[ss] = pair[S[ss][i]][S[ss][j]];
 765       if (type[ss]==0) type[ss] = 7;
 766     }
 767
 768     if (noLonelyPairs)
 769       if (cij == c[indx[j]+i]) {
 770         /* (i.j) closes canonical structures, thus
 771            (i+1.j-1) must be a pair                */
 772         for (ss=0; ss<n_seq; ss++) {
 773           type_2 = pair[S[ss][j-1]][S[ss][i+1]];  /* j,i not i,j */
 774           if (type_2==0) type_2 = 7;
 775           cij -= P->stack[type[ss]][type_2];
 776         }
 777         cij += pscore[indx[j]+i];
 778         base_pair2[++b].i = i+1;
 779         base_pair2[b].j   = j-1;
 780         cov_en += pscore[indx[j-1]+i+1];
 781         i++; j--;
 782         canonical=0;
 783         goto repeat1;
 784       }
 785     canonical = 1;
 786     cij += pscore[indx[j]+i];
 787
 788     {int cc=0;
 789     for (ss=0; ss<n_seq; ss++) {
 790         if ((a2s[ss][j-1]-a2s[ss][i])<3) cc+=600;
 791         else cc += E_Hairpin(a2s[ss][j-1]-a2s[ss][i], type[ss], S3[ss][i], S5[ss][j], Ss[ss]+a2s[ss][i-1], P);
 792       }
 793     if (cij == cc) /* found hairpin */
 794       continue;
 795     }
 796     for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
 797       minq = j-i+p-MAXLOOP-2;
 798       if (minq<p+1+TURN) minq = p+1+TURN;
 799       for (q = j-1; q >= minq; q--) {
 800
 801         if (c[indx[q]+p]>=INF) continue;
 802
 803         for (ss=energy=0; ss<n_seq; ss++) {
 804           type_2 = pair[S[ss][q]][S[ss][p]];  /* q,p not p,q */
 805           if (type_2==0) type_2 = 7;
 806           energy += E_IntLoop(a2s[ss][p-1]-a2s[ss][i],a2s[ss][j-1]-a2s[ss][q],
 807                                type[ss], type_2,
 808                                S3[ss][i], S5[ss][j],
 809                                S5[ss][p], S3[ss][q], P);
 810
 811         }
 812         traced = (cij == energy+c[indx[q]+p]);
 813         if (traced) {
 814           base_pair2[++b].i = p;
 815           base_pair2[b].j   = q;
 816           cov_en += pscore[indx[q]+p];
 817           i = p, j = q;
 818           goto repeat1;
 819         }
 820       }
 821     }
 822
 823     /* end of repeat: --------------------------------------------------*/
 824
 825     /* (i.j) must close a multi-loop */
 826
 827     i1 = i+1;
 828     j1 = j-1;
 829
 830     if(with_gquad){
 831       /*
 832         The case that is handled here actually resembles something like
 833         an interior loop where the enclosing base pair is of regular
 834         kind and the enclosed pair is not a canonical one but a g-quadruplex
 835         that should then be decomposed further...
 836       */
 837       mm = 0;
 838       for(s=0;s<n_seq;s++){
 839         tt = type[s];
 840         if(tt == 0) tt = 7;
 841         if(dangles == 2)
 842           mm += P->mismatchI[tt][S3[s][i]][S5[s][j]];
 843         if(tt > 2)
 844           mm += P->TerminalAU;
 845       }
 846
 847       for(p = i + 2;
 848         p < j - VRNA_GQUAD_MIN_BOX_SIZE;
 849         p++){
 850         if(S_cons[p] != 3) continue;
 851         l1    = p - i - 1;
 852         if(l1>MAXLOOP) break;
 853         minq  = j - i + p - MAXLOOP - 2;
 854         c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 855         minq  = MAX2(c0, minq);
 856         c0    = j - 1;
 857         maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 858         maxq  = MIN2(c0, maxq);
 859         for(q = minq; q < maxq; q++){
 860           if(S_cons[q] != 3) continue;
 861           c0  = mm + ggg[indx[q] + p] + n_seq * P->internal_loop[l1 + j - q - 1];
 862           if(cij == c0){
 863             i=p;j=q;
 864             goto repeat_gquad;
 865           }
 866         }
 867       }
 868       p = i1;
 869       if(S_cons[p] == 3){
 870         if(p < j - VRNA_GQUAD_MIN_BOX_SIZE){
 871           minq  = j - i + p - MAXLOOP - 2;
 872           c0    = p + VRNA_GQUAD_MIN_BOX_SIZE - 1;
 873           minq  = MAX2(c0, minq);
 874           c0    = j - 3;
 875           maxq  = p + VRNA_GQUAD_MAX_BOX_SIZE + 1;
 876           maxq  = MIN2(c0, maxq);
 877           for(q = minq; q < maxq; q++){
 878             if(S_cons[q] != 3) continue;
 879             if(cij == mm + ggg[indx[q] + p] + n_seq * P->internal_loop[j - q - 1]){
 880               i = p; j=q;
 881               goto repeat_gquad;
 882             }
 883           }
 884         }
 885       }
 886       q = j1;
 887       if(S_cons[q] == 3)
 888         for(p = i1 + 3; p < j - VRNA_GQUAD_MIN_BOX_SIZE; p++){
 889           l1    = p - i - 1;
 890           if(l1>MAXLOOP) break;
 891           if(S_cons[p] != 3) continue;
 892           if(cij == mm + ggg[indx[q] + p] + n_seq * P->internal_loop[l1]){
 893             i = p; j = q;
 894             goto repeat_gquad;
 895           }
 896         }
 897     }
 898
 899     mm = n_seq*P->MLclosing;
 900     if(dangles){
 901       for(ss = 0; ss < n_seq; ss++){
 902         tt = rtype[type[ss]];
 903         mm += E_MLstem(tt, S5[ss][j], S3[ss][i], P);
 904       }
 905     }
 906     else{
 907       for(ss = 0; ss < n_seq; ss++){
 908         tt = rtype[type[ss]];
 909         mm += E_MLstem(tt, -1, -1, P);
 910       }
 911     }
 912     sector[s+1].ml  = sector[s+2].ml = 1;
 913
 914     for (k = i1+TURN+1; k < j1-TURN-1; k++){
 915       if(cij == fML[indx[k]+i1] + fML[indx[j1]+k+1] + mm) break;
 916     }
 917
 918     if (k<=j-3-TURN) { /* found the decomposition */
 919       sector[++s].i = i1;
 920       sector[s].j   = k;
 921       sector[++s].i = k+1;
 922       sector[s].j   = j1;
 923     } else {
 924         nrerror("backtracking failed in repeat");
 925     }
 926
 927     continue; /* this is a workarround to not accidentally proceed in the following block */
 928
 929   repeat_gquad:
 930     /*
 931       now we do some fancy stuff to backtrace the stacksize and linker lengths
 932       of the g-quadruplex that should reside within position i,j
 933     */
 934     {
 935       int cnt1, cnt2, cnt3, cnt4, l[3], L, size;
 936       size = j-i+1;
 937
 938       for(L=0; L < VRNA_GQUAD_MIN_STACK_SIZE;L++){
 939         if(S_cons[i+L] != 3) break;
 940         if(S_cons[j-L] != 3) break;
 941       }
 942
 943       if(L == VRNA_GQUAD_MIN_STACK_SIZE){
 944         /* continue only if minimum stack size starting from i is possible */
 945         for(; L<=VRNA_GQUAD_MAX_STACK_SIZE;L++){
 946           if(S_cons[i+L-1] != 3) break; /* break if no more consecutive G's 5' */
 947           if(S_cons[j-L+1] != 3) break; /* break if no more consecutive G'1 3' */
 948           for(    l[0] = VRNA_GQUAD_MIN_LINKER_LENGTH;
 949                   (l[0] <= VRNA_GQUAD_MAX_LINKER_LENGTH)
 950               &&  (size - 4*L - 2*VRNA_GQUAD_MIN_LINKER_LENGTH - l[0] >= 0);
 951               l[0]++){
 952             /* check whether we find the second stretch of consecutive G's */
 953             for(cnt1 = 0; (cnt1 < L) && (S_cons[i+L+l[0]+cnt1] == 3); cnt1++);
 954             if(cnt1 < L) continue;
 955             for(    l[1] = VRNA_GQUAD_MIN_LINKER_LENGTH;
 956                     (l[1] <= VRNA_GQUAD_MAX_LINKER_LENGTH)
 957                 &&  (size - 4*L - VRNA_GQUAD_MIN_LINKER_LENGTH - l[0] - l[1] >= 0);
 958                 l[1]++){
 959               /* check whether we find the third stretch of consectutive G's */
 960               for(cnt1 = 0; (cnt1 < L) && (S_cons[i+2*L+l[0]+l[1]+cnt1] == 3); cnt1++);
 961               if(cnt1 < L) continue;
 962
 963               /*
 964                 the length of the third linker now depends on position j as well
 965                 as the other linker lengths... so we do not have to loop too much
 966               */
 967               l[2] = size - 4*L - l[0] - l[1];
 968               if(l[2] < VRNA_GQUAD_MIN_LINKER_LENGTH) break;
 969               if(l[2] > VRNA_GQUAD_MAX_LINKER_LENGTH) continue;
 970               /* check for contribution */
 971               if(ggg[indx[j]+i] == E_gquad_ali(i, L, l, (const short **)S, n_seq, P)){
 972                 int a;
 973                 /* fill the G's of the quadruplex into base_pair2 */
 974                 for(a=0;a<L;a++){
 975                   base_pair2[++b].i = i+a;
 976                   base_pair2[b].j   = i+a;
 977                   base_pair2[++b].i = i+L+l[0]+a;
 978                   base_pair2[b].j   = i+L+l[0]+a;
 979                   base_pair2[++b].i = i+L+l[0]+L+l[1]+a;
 980                   base_pair2[b].j   = i+L+l[0]+L+l[1]+a;
 981                   base_pair2[++b].i = i+L+l[0]+L+l[1]+L+l[2]+a;
 982                   base_pair2[b].j   = i+L+l[0]+L+l[1]+L+l[2]+a;
 983                 }
 984                 goto repeat_gquad_exit;
 985               }
 986             }
 987           }
 988         }
 989       }
 990       nrerror("backtracking failed in repeat_gquad");
 991     }
 992   repeat_gquad_exit:
 993     asm("nop");
 994
 995   }
 996
 997   /* fprintf(stderr, "covariance energy %6.2f\n", cov_en/100.);  */
 998
 999   base_pair2[0].i = b;    /* save the total number of base pairs */
1000   free(type);
1001 }
1002
1003
1004 PUBLIC void encode_ali_sequence(const char *sequence, short *S, short *s5, short *s3, char *ss, unsigned short *as, int circular){
1005   unsigned int i,l;
1006   unsigned short p;
1007   l     = strlen(sequence);
1008   S[0]  = (short) l;
1009   s5[0] = s5[1] = 0;
1010
1011   /* make numerical encoding of sequence */
1012   for(i=1; i<=l; i++){
1013     short ctemp;
1014     ctemp=(short) encode_char(toupper(sequence[i-1]));
1015     S[i]= ctemp ;
1016   }
1017
1018   if (oldAliEn){
1019     /* use alignment sequences in all energy evaluations */
1020     ss[0]=sequence[0];
1021     for(i=1; i<l; i++){
1022       s5[i] = S[i-1];
1023       s3[i] = S[i+1];
1024       ss[i] = sequence[i];
1025       as[i] = i;
1026     }
1027     ss[l]   = sequence[l];
1028     as[l]   = l;
1029     s5[l]   = S[l-1];
1030     s3[l]   = 0;
1031     S[l+1]  = S[1];
1032     s5[1]   = 0;
1033     if (circular) {
1034       s5[1]   = S[l];
1035       s3[l]   = S[1];
1036       ss[l+1] = S[1];
1037     }
1038   }
1039   else{
1040     if(circular){
1041       for(i=l; i>0; i--){
1042         char c5;
1043         c5 = sequence[i-1];
1044         if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.')) continue;
1045         s5[1] = S[i];
1046         break;
1047       }
1048       for (i=1; i<=l; i++) {
1049         char c3;
1050         c3 = sequence[i-1];
1051         if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.')) continue;
1052         s3[l] = S[i];
1053         break;
1054       }
1055     }
1056     else  s5[1]=s3[l]=0;
1057
1058     for(i=1,p=0; i<=l; i++){
1059       char c5;
1060       c5 = sequence[i-1];
1061       if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.'))
1062         s5[i+1]=s5[i];
1063       else { /* no gap */
1064         ss[p++]=sequence[i-1]; /*start at 0!!*/
1065         s5[i+1]=S[i];
1066       }
1067       as[i]=p;
1068     }
1069     for (i=l; i>=1; i--) {
1070       char c3;
1071       c3 = sequence[i-1];
1072       if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.'))
1073         s3[i-1]=s3[i];
1074       else
1075         s3[i-1]=S[i];
1076     }
1077   }
1078 }
1079
1080 PRIVATE void make_pscores(const short *const* S, const char **AS,
1081                           int n_seq, const char *structure) {
1082   /* calculate co-variance bonus for each pair depending on  */
1083   /* compensatory/consistent mutations and incompatible seqs */
1084   /* should be 0 for conserved pairs, >0 for good pairs      */
1085 #define NONE -10000 /* score for forbidden pairs */
1086   int n,i,j,k,l,s;
1087
1088   int olddm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
1089                   {0,0,2,2,1,2,2} /* CG */,
1090                   {0,2,0,1,2,2,2} /* GC */,
1091                   {0,2,1,0,2,1,2} /* GU */,
1092                   {0,1,2,2,0,2,1} /* UG */,
1093                   {0,2,2,1,2,0,2} /* AU */,
1094                   {0,2,2,2,1,2,0} /* UA */};
1095
1096   float **dm;
1097   n=S[0][0];  /* length of seqs */
1098   if (ribo) {
1099     if (RibosumFile !=NULL) dm=readribosum(RibosumFile);
1100     else dm=get_ribosum(AS,n_seq,n);
1101   }
1102   else { /*use usual matrix*/
1103     dm=(float **)space(7*sizeof(float*));
1104     for (i=0; i<7;i++) {
1105       dm[i]=(float *)space(7*sizeof(float));
1106       for (j=0; j<7; j++)
1107         dm[i][j] = (float) olddm[i][j];
1108     }
1109   }
1110
1111   for (i=1; i<n; i++) {
1112     for (j=i+1; (j<i+TURN+1) && (j<=n); j++)
1113       pscore[indx[j]+i] = NONE;
1114     for (j=i+TURN+1; j<=n; j++) {
1115       int pfreq[8]={0,0,0,0,0,0,0,0};
1116       double score;
1117       for (s=0; s<n_seq; s++) {
1118         int type;
1119         if (S[s][i]==0 && S[s][j]==0) type = 7; /* gap-gap  */
1120         else {
1121           if ((AS[s][i] == '~')||(AS[s][j] == '~')) type = 7;
1122           else type = pair[S[s][i]][S[s][j]];
1123         }
1124         pfreq[type]++;
1125       }
1126       if (pfreq[0]*2+pfreq[7]>n_seq) { pscore[indx[j]+i] = NONE; continue;}
1127       for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
1128         for (l=k; l<=6; l++)
1129           score += pfreq[k]*pfreq[l]*dm[k][l];
1130       /* counter examples score -1, gap-gap scores -0.25   */
1131       pscore[indx[j]+i] = cv_fact *
1132         ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
1133     }
1134   }
1135
1136   if (noLonelyPairs) /* remove unwanted pairs */
1137     for (k=1; k<n-TURN-1; k++)
1138       for (l=1; l<=2; l++) {
1139         int type,ntype=0,otype=0;
1140         i=k; j = i+TURN+l;
1141         type = pscore[indx[j]+i];
1142         while ((i>=1)&&(j<=n)) {
1143           if ((i>1)&&(j<n)) ntype = pscore[indx[j+1]+i-1];
1144           if ((otype<cv_fact*MINPSCORE)&&(ntype<cv_fact*MINPSCORE))  /* too many counterexamples */
1145             pscore[indx[j]+i] = NONE; /* i.j can only form isolated pairs */
1146           otype =  type;
1147           type  = ntype;
1148           i--; j++;
1149         }
1150       }
1151
1152
1153   if (fold_constrained&&(structure!=NULL)) {
1154     int psij, hx, hx2, *stack, *stack2;
1155     stack = (int *) space(sizeof(int)*(n+1));
1156     stack2 = (int *) space(sizeof(int)*(n+1));
1157
1158     for(hx=hx2=0, j=1; j<=n; j++) {
1159       switch (structure[j-1]) {
1160       case 'x': /* can't pair */
1161         for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
1162         for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
1163         break;
1164       case '(':
1165         stack[hx++]=j;
1166         /* fallthrough */
1167       case '[':
1168         stack2[hx2++]=j;
1169         /* fallthrough */
1170       case '<': /* pairs upstream */
1171         for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
1172         break;
1173       case ']':
1174         if (hx2<=0) {
1175           fprintf(stderr, "%s\n", structure);
1176           nrerror("unbalanced brackets in constraints");
1177         }
1178         i = stack2[--hx2];
1179         pscore[indx[j]+i]=NONE;
1180         break;
1181       case ')':
1182         if (hx<=0) {
1183           fprintf(stderr, "%s\n", structure);
1184           nrerror("unbalanced brackets in constraints");
1185         }
1186         i = stack[--hx];
1187         psij = pscore[indx[j]+i]; /* store for later */
1188         for (k=j; k<=n; k++)
1189           for (l=i; l<=j; l++)
1190             pscore[indx[k]+l] = NONE;
1191         for (l=i; l<=j; l++)
1192           for (k=1; k<=i; k++)
1193             pscore[indx[l]+k] = NONE;
1194         for (k=i+1; k<j; k++)
1195           pscore[indx[k]+i] = pscore[indx[j]+k] = NONE;
1196         pscore[indx[j]+i] = (psij>0) ? psij : 0;
1197         /* fallthrough */
1198       case '>': /* pairs downstream */
1199         for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
1200         break;
1201       }
1202     }
1203     if (hx!=0) {
1204       fprintf(stderr, "%s\n", structure);
1205       nrerror("unbalanced brackets in constraint string");
1206     }
1207     free(stack); free(stack2);
1208   }
1209   /*free dm */
1210   for (i=0; i<7;i++) {
1211     free(dm[i]);
1212   }
1213   free(dm);
1214 }
1215
1216 /*--------New scoring part-----------------------------------*/
1217
1218 PUBLIC int get_mpi(char *Alseq[], int n_seq, int length, int *mini) {
1219   int i, j,k;
1220   float ident=0;
1221   int pairnum=0;
1222   int sumident=0;
1223   float minimum=1.;
1224   for(j=0; j<n_seq-1; j++)
1225     for(k=j+1; k<n_seq; k++) {
1226       ident=0;
1227       for (i=1; i<=length; i++){
1228         if (Alseq[k][i]==Alseq[j][i]) ident++;
1229         pairnum++;
1230       }
1231       if ((ident/length)<minimum) minimum=ident/(float)length;
1232       sumident+=ident;
1233     }
1234   mini[0]=(int)(minimum*100.);
1235   if (pairnum>0)   return (int) (sumident*100/pairnum);
1236   else return 0;
1237
1238 }
1239
1240 PUBLIC float **readribosum(char *name){
1241
1242   float **dm;
1243   char *line;
1244   FILE *fp;
1245   int i=0;
1246   int who=0;
1247   float a,b,c,d,e,f;
1248   int translator[7]={0,5,1,2,3,6,4};
1249
1250   fp=fopen(name,"r");
1251   dm=(float **)space(7*sizeof(float*));
1252   for (i=0; i<7;i++) {
1253     dm[i]=(float *)space(7*sizeof(float));
1254   }
1255   while(1) { /*bisma hoit fertisch san*/
1256     line=get_line(fp);
1257     if (*line=='#') continue;
1258     i=0;
1259     i=sscanf(line,"%f %f %f %f %f %f",&a,&b,&c,&d,&e,&f);
1260     if (i==0) break;
1261     dm[translator[++who]][translator[1]]=a;
1262     dm[translator[who]][translator[2]]=b;
1263     dm[translator[who]][translator[3]]=c;
1264     dm[translator[who]][translator[4]]=d;
1265     dm[translator[who]][translator[5]]=e;
1266     dm[translator[who]][translator[6]]=f;
1267     free(line);
1268     if (who==6) break;
1269   }
1270   fclose(fp);
1271   return dm;
1272 }
1273
1274
1275 PRIVATE void en_corr_of_loop_gquad(int i,
1276                                   int j,
1277                                   const char **sequences,
1278                                   const char *structure,
1279                                   short *pt,
1280                                   int *loop_idx,
1281                                   int n_seq,
1282                                   int en[2]){
1283
1284   int pos, energy, en_covar, p, q, r, s, u, type, type2, gq_en[2];
1285   int L, l[3];
1286
1287   energy = en_covar = 0;
1288   q = i;
1289   while((pos = parse_gquad(structure + q-1, &L, l)) > 0){
1290     q += pos-1;
1291     p = q - 4*L - l[0] - l[1] - l[2] + 1;
1292     if(q > j) break;
1293     /* we've found the first g-quadruplex at position [p,q] */
1294     E_gquad_ali_en(p, L, l, (const short **)S, n_seq, gq_en, P);
1295     energy    += gq_en[0];
1296     en_covar  += gq_en[1];
1297     /* check if it's enclosed in a base pair */
1298     if(loop_idx[p] == 0){ q++; continue; /* g-quad in exterior loop */}
1299     else{
1300       energy += E_MLstem(0, -1, -1, P) * n_seq;
1301       /*  find its enclosing pair */
1302       int num_elem, num_g, elem_i, elem_j, up_mis;
1303       num_elem  = 0;
1304       num_g     = 1;
1305       r         = p - 1;
1306       up_mis    = q - p + 1;
1307
1308       /* seek for first pairing base located 5' of the g-quad */
1309       for(r = p - 1; !pt[r] && (r >= i); r--);
1310       if(r < i) nrerror("this should not happen");
1311
1312       if(r < pt[r]){ /* found the enclosing pair */
1313         s = pt[r];
1314       } else {
1315         num_elem++;
1316         elem_i = pt[r];
1317         elem_j = r;
1318         r = pt[r]-1 ;
1319         /* seek for next pairing base 5' of r */
1320         for(; !pt[r] && (r >= i); r--);
1321         if(r < i) nrerror("so nich");
1322         if(r < pt[r]){ /* found the enclosing pair */
1323           s = pt[r];
1324         } else {
1325           /* hop over stems and unpaired nucleotides */
1326           while((r > pt[r]) && (r >= i)){
1327             if(pt[r]){ r = pt[r]; num_elem++;}
1328             r--;
1329           }
1330           if(r < i) nrerror("so nich");
1331           s = pt[r]; /* found the enclosing pair */
1332         }
1333       }
1334       /* now we have the enclosing pair (r,s) */
1335
1336       u = q+1;
1337       /* we know everything about the 5' part of this loop so check the 3' part */
1338       while(u<s){
1339         if(structure[u-1] == '.') u++;
1340         else if (structure[u-1] == '+'){ /* found another gquad */
1341           pos = parse_gquad(structure + u - 1, &L, l);
1342           if(pos > 0){
1343             E_gquad_ali_en(u, L, l, (const short **)S, n_seq, gq_en, P);
1344             energy += gq_en[0] + E_MLstem(0, -1, -1, P) * n_seq;
1345             en_covar += gq_en[1];
1346             up_mis += pos;
1347             u += pos;
1348             num_g++;
1349           }
1350         } else { /* we must have found a stem */
1351           if(!(u < pt[u])) nrerror("wtf!");
1352           num_elem++;
1353           elem_i = u;
1354           elem_j = pt[u];
1355           en_corr_of_loop_gquad(u,
1356                                 pt[u],
1357                                 sequences,
1358                                 structure,
1359                                 pt,
1360                                 loop_idx,
1361                                 n_seq,
1362                                 gq_en);
1363           energy    += gq_en[0];
1364           en_covar  += gq_en[1];
1365           u = pt[u] + 1;
1366         }
1367       }
1368       if(u!=s) nrerror("what the ...");
1369       else{ /* we are done since we've found no other 3' structure element */
1370         switch(num_elem){
1371           /* g-quad was misinterpreted as hairpin closed by (r,s) */
1372           case 0:   /*if(num_g == 1)
1373                       if((p-r-1 == 0) || (s-q-1 == 0))
1374                         nrerror("too few unpaired bases");
1375                     */
1376                     {
1377                       int ee = 0;
1378                       int cnt;
1379                       for(cnt=0;cnt<n_seq;cnt++){
1380                         type = pair[S[cnt][r]][S[cnt][s]];
1381                         if(type == 0) type = 7;
1382                         if ((a2s[cnt][s-1]-a2s[cnt][r])<3) ee+=600;
1383                         else ee += E_Hairpin( a2s[cnt][s-1] - a2s[cnt][r],
1384                                               type,
1385                                               S3[cnt][r],
1386                                               S5[cnt][s],
1387                                               Ss[cnt] + a2s[cnt][r-1],
1388                                               P);
1389                       }
1390                       energy -= ee;
1391                       ee = 0;
1392                       for(cnt=0;cnt < n_seq; cnt++){
1393                         type = pair[S[cnt][r]][S[cnt][s]];
1394                         if(type == 0) type = 7;
1395                         if(dangles == 2)
1396                           ee += P->mismatchI[type][S3[cnt][r]][S5[cnt][s]];
1397                         if(type > 2)
1398                           ee += P->TerminalAU;
1399                       }
1400                       energy += ee;
1401                     }
1402                     energy += n_seq * P->internal_loop[s-r-1-up_mis];
1403                     break;
1404           /* g-quad was misinterpreted as interior loop closed by (r,s) with enclosed pair (elem_i, elem_j) */
1405           case 1:   {
1406                       int ee = 0;
1407                       int cnt;
1408                       for(cnt = 0; cnt<n_seq;cnt++){
1409                         type = pair[S[cnt][r]][S[cnt][s]];
1410                         if(type == 0) type = 7;
1411                         type2 = pair[S[cnt][elem_j]][S[cnt][elem_i]];
1412                         if(type2 == 0) type2 = 7;
1413                         ee += E_IntLoop(a2s[cnt][elem_i-1] - a2s[cnt][r],
1414                                         a2s[cnt][s-1] - a2s[cnt][elem_j],
1415                                         type,
1416                                         type2,
1417                                         S3[cnt][r],
1418                                         S5[cnt][s],
1419                                         S5[cnt][elem_i],
1420                                         S3[cnt][elem_j],
1421                                         P);
1422                       }
1423                       energy -= ee;
1424                       ee = 0;
1425                       for(cnt = 0; cnt < n_seq; cnt++){
1426                         type = pair[S[cnt][s]][S[cnt][r]];
1427                         if(type == 0) type = 7;
1428                         ee += E_MLstem(type, S5[cnt][s], S3[cnt][r], P);
1429                         type = pair[S[cnt][elem_i]][S[cnt][elem_j]];
1430                         if(type == 0) type = 7;
1431                         ee += E_MLstem(type, S5[cnt][elem_i], S3[cnt][elem_j], P);
1432                       }
1433                       energy += ee;
1434                     }
1435                     energy += (P->MLclosing + (elem_i-r-1+s-elem_j-1-up_mis) * P->MLbase) * n_seq;
1436                     break;
1437           /* gquad was misinterpreted as unpaired nucleotides in a multiloop */
1438           default:  energy -= (up_mis) * P->MLbase * n_seq;
1439                     break;
1440         }
1441       }
1442       q = s+1;
1443     }
1444   }
1445   en[0] = energy;
1446   en[1] = en_covar;
1447 }
1448
1449 PUBLIC float energy_of_ali_gquad_structure( const char **sequences,
1450                                             const char *structure,
1451                                             int n_seq,
1452                                             float *energy){
1453
1454   int new=0;
1455   unsigned int s, n;
1456   short *pt;
1457
1458   short           **tempS;
1459   short           **tempS5;     /*S5[s][i] holds next base 5' of i in sequence s*/
1460   short           **tempS3;     /*Sl[s][i] holds next base 3' of i in sequence s*/
1461   char            **tempSs;
1462   unsigned short  **tempa2s;
1463
1464   int *tempindx, *loop_idx;
1465   int *temppscore;
1466
1467   int en_struct[2], gge[2];
1468
1469   if(sequences[0] != NULL){
1470     n = (unsigned int) strlen(sequences[0]);
1471     update_alifold_params();
1472
1473     /*save old memory*/
1474     tempS = S; tempS3 = S3; tempS5 = S5; tempSs = Ss; tempa2s = a2s;
1475     tempindx = indx; temppscore = pscore;
1476
1477     alloc_sequence_arrays(sequences, &S, &S5, &S3, &a2s, &Ss, 0);
1478     pscore  = (int *) space(sizeof(int)*((n+1)*(n+2)/2));
1479     indx    = get_indx(n);
1480     make_pscores((const short *const*)S, sequences, n_seq, NULL);
1481     new     = 1;
1482
1483     pt = make_pair_table(structure);
1484     energy_of_alistruct_pt(sequences,pt, n_seq, &(en_struct[0]));
1485     loop_idx    = make_loop_index_pt(pt);
1486     en_corr_of_loop_gquad(1, n, sequences, structure, pt, loop_idx, n_seq, gge);
1487     en_struct[0] += gge[0];
1488     en_struct[1] += gge[1];
1489
1490     free(loop_idx);
1491     free(pt);
1492     energy[0] = (float)en_struct[0]/(float)(100*n_seq);
1493     energy[1] = (float)en_struct[1]/(float)(100*n_seq);
1494
1495     free(pscore);
1496     free(indx);
1497     free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);
1498
1499     /* restore old memory */
1500     S = tempS; S3 = tempS3; S5 = tempS5; Ss = tempSs; a2s = tempa2s;
1501     indx = tempindx; pscore = temppscore;
1502   }
1503   else nrerror("energy_of_alistruct(): no sequences in alignment!");
1504
1505   return energy[0];
1506
1507 }
1508
1509 PUBLIC  float energy_of_alistruct(const char **sequences, const char *structure, int n_seq, float *energy){
1510   int new=0;
1511   unsigned int s, n;
1512   short *pt;
1513
1514   short           **tempS;
1515   short           **tempS5;     /*S5[s][i] holds next base 5' of i in sequence s*/
1516   short           **tempS3;     /*Sl[s][i] holds next base 3' of i in sequence s*/
1517   char            **tempSs;
1518   unsigned short  **tempa2s;
1519
1520   int *tempindx;
1521   int *temppscore;
1522
1523   int en_struct[2];
1524
1525   if(sequences[0] != NULL){
1526     n = (unsigned int) strlen(sequences[0]);
1527     update_alifold_params();
1528
1529     /*save old memory*/
1530     tempS = S; tempS3 = S3; tempS5 = S5; tempSs = Ss; tempa2s = a2s;
1531     tempindx = indx; temppscore = pscore;
1532
1533     alloc_sequence_arrays(sequences, &S, &S5, &S3, &a2s, &Ss, 0);
1534     pscore  = (int *) space(sizeof(int)*((n+1)*(n+2)/2));
1535     indx    = get_indx(n);
1536     make_pscores((const short *const*)S, sequences, n_seq, NULL);
1537     new     = 1;
1538
1539     pt = make_pair_table(structure);
1540     energy_of_alistruct_pt(sequences,pt, n_seq, &(en_struct[0]));
1541
1542     free(pt);
1543     energy[0] = (float)en_struct[0]/(float)(100*n_seq);
1544     energy[1] = (float)en_struct[1]/(float)(100*n_seq);
1545
1546     free(pscore);
1547     free(indx);
1548     free_sequence_arrays(n_seq, &S, &S5, &S3, &a2s, &Ss);
1549
1550     /* restore old memory */
1551     S = tempS; S3 = tempS3; S5 = tempS5; Ss = tempSs; a2s = tempa2s;
1552     indx = tempindx; pscore = temppscore;
1553   }
1554   else nrerror("energy_of_alistruct(): no sequences in alignment!");
1555
1556   return energy[0];
1557 }
1558
1559 PRIVATE void energy_of_alistruct_pt(const char **sequences, short *pt, int n_seq, int *energy){
1560   int i, length;
1561
1562   length = S[0][0];
1563   energy[0] =  backtrack_type=='M' ? ML_Energy_pt(0, n_seq, pt) : EL_Energy_pt(0, n_seq, pt);
1564   energy[1] = 0;
1565   for (i=1; i<=length; i++) {
1566     if (pt[i]==0) continue;
1567     stack_energy_pt(i, sequences, pt, n_seq, energy);
1568     i=pt[i];
1569   }
1570 }
1571
1572 PRIVATE void stack_energy_pt(int i, const char **sequences, short *pt, int n_seq, int *energy)
1573 {
1574   /* calculate energy of substructure enclosed by (i,j) */
1575   int ee= 0;
1576   int j, p, q, s;
1577   int *type = (int *) space(n_seq*sizeof(int));
1578
1579   j = pt[i];
1580   for (s=0; s<n_seq; s++) {
1581     type[s] = pair[S[s][i]][S[s][j]];
1582     if (type[s]==0) {
1583     type[s]=7;
1584     }
1585   }
1586   p=i; q=j;
1587   while (p<q) { /* process all stacks and interior loops */
1588     int type_2;
1589     while (pt[++p]==0);
1590     while (pt[--q]==0);
1591     if ((pt[q]!=(short)p)||(p>q)) break;
1592     ee=0;
1593     for (s=0; s<n_seq; s++) {
1594       type_2 = pair[S[s][q]][S[s][p]];
1595       if (type_2==0) {
1596         type_2=7;
1597       }
1598       ee += E_IntLoop(a2s[s][p-1]-a2s[s][i], a2s[s][j-1]-a2s[s][q], type[s], type_2, S3[s][i], S5[s][j], S5[s][p], S3[s][q], P);
1599     }
1600     energy[0] += ee;
1601     energy[1] += pscore[indx[j]+i];
1602     i=p; j=q;
1603     for (s=0; s<n_seq; s++) {
1604       type[s] = pair[S[s][i]][S[s][j]];
1605       if (type[s]==0) type[s]=7;
1606     }
1607   }  /* end while */
1608
1609   /* p,q don't pair must have found hairpin or multiloop */
1610
1611   if (p>q) {
1612     ee=0;/* hair pin */
1613     for (s=0; s< n_seq; s++) {
1614       if ((a2s[s][j-1]-a2s[s][i])<3) ee+=600;
1615       else ee += E_Hairpin(a2s[s][j-1]-a2s[s][i], type[s], S3[s][i], S5[s][j], Ss[s]+(a2s[s][i-1]), P);
1616     }
1617     energy[0] += ee;
1618     energy[1] += pscore[indx[j]+i];
1619     free(type);
1620     return;
1621   }
1622   /* (i,j) is exterior pair of multiloop */
1623   energy[1] += pscore[indx[j]+i];
1624   while (p<j) {
1625     /* add up the contributions of the substructures of the ML */
1626     stack_energy_pt(p, sequences, pt, n_seq, energy);
1627     p = pt[p];
1628     /* search for next base pair in multiloop */
1629     while (pt[++p]==0);
1630   }
1631   energy[0] += ML_Energy_pt(i, n_seq, pt);
1632   free(type);
1633 }
1634
1635
1636
1637 PRIVATE int ML_Energy_pt(int i, int n_seq, short *pt){
1638   /* i is the 5'-base of the closing pair */
1639
1640   int   energy, tt, i1, j, p, q, u, s;
1641   short d5, d3;
1642
1643   j = pt[i];
1644   i1  = i;
1645   p   = i+1;
1646   u   = 0;
1647   energy = 0;
1648
1649   do{ /* walk around the multi-loop */
1650     /* hop over unpaired positions */
1651     while (p < j && pt[p]==0) p++;
1652     if(p >= j) break;
1653     /* get position of pairing partner */
1654     q  = pt[p];
1655     /* memorize number of unpaired positions? no, we just approximate here */
1656     u += p-i1-1;
1657
1658     for (s=0; s< n_seq; s++) {
1659       /* get type of base pair P->q */
1660       tt = pair[S[s][p]][S[s][q]];
1661       if (tt==0) tt=7;
1662       d5 = dangles && (a2s[s][p]>1) && (tt!=0) ? S5[s][p] : -1;
1663       d3 = dangles && (a2s[s][q]<a2s[s][S[0][0]]) ? S3[s][q] : -1;
1664       energy += E_MLstem(tt, d5, d3, P);
1665     }
1666     i1  = q;
1667     p   = q + 1;
1668   }while(1);
1669
1670   if(i > 0){
1671     energy  += P->MLclosing * n_seq;
1672     if(dangles){
1673       for (s=0; s<n_seq; s++){
1674         tt = pair[S[s][j]][S[s][i]];
1675         if (tt==0) tt=7;
1676         energy += E_MLstem(tt, S5[s][j], S3[s][i], P);
1677       }
1678     }
1679     else{
1680       for (s=0; s<n_seq; s++){
1681         tt = pair[S[s][j]][S[s][i]];
1682         if (tt==0) tt=7;
1683         energy += E_MLstem(tt, -1, -1, P);
1684       }
1685     }
1686   }
1687   u += j - i1 - 1;
1688   energy += u * P->MLbase * n_seq;
1689   return energy;
1690 }
1691
1692 PRIVATE int EL_Energy_pt(int i, int n_seq, short *pt){
1693   int   energy, tt, i1, j, p, q, s;
1694   short d5, d3;
1695
1696   j = pt[0];
1697
1698   p   = i+1;
1699   energy = 0;
1700
1701   do{ /* walk along the backbone */
1702     /* hop over unpaired positions */
1703     while (p < j && pt[p]==0) p++;
1704     if(p == j) break; /* no more stems */
1705     /* get position of pairing partner */
1706     q  = pt[p];
1707     for (s=0; s< n_seq; s++) {
1708       /* get type of base pair P->q */
1709       tt = pair[S[s][p]][S[s][q]];
1710       if (tt==0) tt=7;
1711       d5 = dangles && (a2s[s][p]>1) && (tt!=0) ? S5[s][p] : -1;
1712       d3 = dangles && (a2s[s][q]<a2s[s][S[0][0]]) ? S3[s][q] : -1;
1713       energy += E_ExtLoop(tt, d5, d3, P);
1714     }
1715     p   = q + 1;
1716   }while(p < j);
1717
1718   return energy;
1719 }
1720