binaries/src/ViennaRNA/lib/ali_plex.c

   1 /* Last changed Time-stamp: <2007-10-30 14:06:22 htafer> */
   2 /*
   3            compute the duplex structure of two RNA strands,
   4                 allowing only inter-strand base pairs.
   5          see cofold() for computing hybrid structures without
   6                              restriction.
   7                              Ivo Hofacker
   8                           Vienna RNA package
   9
  10 */
  11
  12
  13 /*
  14   library containing the function used in rnaplex
  15   the program rnaplex uses the following function
  16   Lduplexfold: finds high scoring segments
  17   it stores the end-position of these segments in an array
  18   and call then for each of these positions the duplexfold function
  19   which allows to make backtracking for each of the high scoring position
  20   It allows to find suboptimal partially overlapping (depends on a a parameter)
  21   duplexes between a long RNA and a shorter one.
  22   Contrarly to RNAduplex, the energy model is not in E~log(N),
  23   where N is the length of an interial loop but used an affine model,
  24   where the extension and begin parameter are fitted to the energy
  25   parameter used by RNAduplex. This allows to check for duplex between a short RNA(20nt)
  26   and a long one at the speed of 1Mnt/s. At this speed the whole genome (3Gnt) can be analyzed for one siRNA
  27   in about 50 minutes.
  28   The algorithm is based on an idea by Durbin and Eddy:when the alginment reach a value larger than a
  29   given threshold this value is stored in an array. When the alignment score goes
  30   then under this threshold, the alignemnent begin from this value, in that way the backtracking allow us
  31   to find all non-overlapping high-scoring segments.
  32   For more information check "durbin, biological sequence analysis"
  33 */
  34
  35 #include <config.h>
  36 #include <stdio.h>
  37 #include <stdlib.h>
  38 #include <math.h>
  39 #include <ctype.h>
  40 #include <string.h>
  41 #include "utils.h"
  42 #include "energy_par.h"
  43 #include "fold_vars.h"
  44 #include "fold.h"
  45 #include "pair_mat.h"
  46 #include "params.h"
  47 #include "loop_energies.h"
  48 #include "plex.h"
  49 #include "ali_plex.h"
  50
  51 /*@unused@*/
  52 static char rcsid[] UNUSED = "$Id: plex.c,v 1.14 2007/06/12 12:50:16 htafer Exp $";
  53
  54
  55 #define PUBLIC
  56 #define PRIVATE static
  57
  58 #define STACK_BULGE1  1   /* stacking energies for bulges of size 1 */
  59 #define NEW_NINIO     1   /* new asymetry penalty */
  60 #define ARRAY 32          /*array size*/
  61 #define UNIT 100
  62 #define MINPSCORE -2 * UNIT
  63 /**
  64 *** Due to the great similarity between functions,
  65 *** more annotation can be found in plex.c
  66 **/
  67
  68 PRIVATE short *encode_seq(const char *seq);
  69 PRIVATE void  update_dfold_params(void);
  70 /**
  71 *** aliduplexfold(_XS)/alibacktrack(_XS) computes duplex interaction with standard energy and considers extension_cost
  72 *** alifind_max(_XS)/aliplot_max(_XS) find suboptimals and MFE
  73 **/
  74 PRIVATE duplexT aliduplexfold(const char *s1[], const char *s2[], const int extension_cost);
  75 PRIVATE char *    alibacktrack(int i, int j, const short *s1[], const short *s2[], const int extension_cost);
  76 PRIVATE void      alifind_max(const int *position, const int *position_j,const int delta, const int threshold,
  77                            const int alignment_length, const char *s1[], const char *s2[], const int extension_cost, const int fast);
  78 PRIVATE void      aliplot_max(const int max, const int max_pos, const int max_pos_j,
  79                            const int alignment_length, const char *s1[], const char *s2[], const int extension_cost, const int fast);
  80 PRIVATE duplexT aliduplexfold_XS(const char *s1[], const char *s2[],const int **access_s1,
  81                                  const int **access_s2, const int i_pos, const int j_pos, const int threshold,const int i_flag, const int j_flag);
  82 PRIVATE char *    alibacktrack_XS(int i, int j, const short *s1[], const short *s2[], const int** access_s1, const int** access_s2,const int i_flag, const int j_flag);
  83 PRIVATE void  alifind_max_XS(const int *position, const int *position_j,
  84                                  const int delta,  const int threshold, const int alignment_length,
  85                                  const char* s1[], const char* s2[],
  86                                  const int **access_s1, const int **access_s2, const int fast);
  87 PRIVATE void aliplot_max_XS(const int max, const int max_pos, const int max_pos_j,
  88                             const int alignment_length, const char *s1[], const char* s2[],
  89                             const int **access_s1, const int **access_s2, const int fast);
  90
  91 /**
  92 *** computes covariance score
  93 **/
  94
  95 PRIVATE int covscore(const int *types, int n_seq);
  96
  97 extern double cv_fact; /* from alifold.c, default 1 */
  98 extern double nc_fact;
  99
 100
 101 /*@unused@*/
 102
 103 #define MAXSECTORS      500     /* dimension for a backtrack array */
 104 #define LOCALITY        0.      /* locality parameter for base-pairs */
 105
 106 PRIVATE paramT *P = NULL;
 107 PRIVATE int   **c = NULL;
 108 PRIVATE int  **lc = NULL, **lin = NULL, **lbx = NULL, **lby = NULL,**linx = NULL, **liny = NULL;
 109
 110
 111
 112
 113 PRIVATE int   n1,n2;
 114 PRIVATE int n3, n4;
 115 PRIVATE int delay_free=0;
 116
 117
 118 /*-----------------------------------------------------------------------duplexfold_XS---------------------------------------------------------------------------*/
 119
 120
 121
 122 /*----------------------------------------------ALIDUPLEXFOLD-----------------------------------------------------------------------------------------------------------*/
 123 PRIVATE duplexT aliduplexfold(const char *s1[], const char *s2[], const int extension_cost) {
 124   int i, j, s, n_seq, l1, Emin=INF, i_min=0, j_min=0;
 125   char *struc;
 126   duplexT mfe;
 127   short **S1, **S2;
 128   int *type;
 129   n3 = (int) strlen(s1[0]);
 130   n4 = (int) strlen(s2[0]);
 131   for (s=0; s1[s]!=NULL; s++);
 132   n_seq = s;
 133   for (s=0; s2[s]!=NULL; s++);
 134   if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
 135
 136   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
 137     update_fold_params();  if(P) free(P); P = scale_parameters();
 138     make_pair_matrix();
 139   }
 140
 141   c = (int **) space(sizeof(int *) * (n3+1));
 142   for (i=1; i<=n3; i++) c[i] = (int *) space(sizeof(int) * (n4+1));
 143
 144   S1 = (short **) space((n_seq+1)*sizeof(short *));
 145   S2 = (short **) space((n_seq+1)*sizeof(short *));
 146   for (s=0; s<n_seq; s++) {
 147     if (strlen(s1[s]) != n3) nrerror("uneqal seqence lengths");
 148     if (strlen(s2[s]) != n4) nrerror("uneqal seqence lengths");
 149     S1[s] = encode_seq(s1[s]);
 150     S2[s] = encode_seq(s2[s]);
 151   }
 152   type = (int *) space(n_seq*sizeof(int));
 153
 154   for (i=1; i<=n3; i++) {
 155     for (j=n4; j>0; j--) {
 156       int k,l,E,psc;
 157       for (s=0; s<n_seq; s++) {
 158         type[s] = pair[S1[s][i]][S2[s][j]];
 159       }
 160       psc = covscore(type, n_seq);
 161       for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
 162       c[i][j] = (psc>=MINPSCORE) ? (n_seq*(P->DuplexInit + 2*extension_cost)) : INF;
 163       if (psc<MINPSCORE) continue;
 164       for (s=0; s<n_seq; s++) {
 165         c[i][j] += E_ExtLoop(type[s], (i>1) ? S1[s][i-1] : -1, (j<n4) ? S2[s][j+1] : -1, P) + 2*extension_cost;
 166       }
 167       for (k=i-1; k>0 && k>i-MAXLOOP-2; k--) {
 168         for (l=j+1; l<=n4; l++) {
 169           int type2;
 170           if (i-k+l-j-2>MAXLOOP) break;
 171           if (c[k][l]>INF/2) continue;
 172           for (E=s=0; s<n_seq; s++) {
 173             type2 = pair[S1[s][k]][S2[s][l]];
 174             if (type2==0) type2=7;
 175             E += E_IntLoop(i-k-1, l-j-1, type2, rtype[type[s]],
 176                            S1[s][k+1], S2[s][l-1], S1[s][i-1], S2[s][j+1],P) + (i-k+l-j)*extension_cost;
 177           }
 178           c[i][j] = MIN2(c[i][j], c[k][l]+E);
 179         }
 180       }
 181       c[i][j] -= psc;
 182       E = c[i][j];
 183       for (s=0; s<n_seq; s++) {
 184         E += E_ExtLoop(rtype[type[s]], (j>1) ? S2[s][j-1] : -1, (i<n3) ? S1[s][i+1] : -1, P) +2*extension_cost;
 185       }
 186       if (E<Emin) {
 187         Emin=E; i_min=i; j_min=j;
 188       }
 189     }
 190   }
 191   struc = alibacktrack(i_min, j_min, (const short int**) S1, (const short int**) S2 , extension_cost);
 192   if (i_min<n3) i_min++;
 193   if (j_min>1 ) j_min--;
 194   l1 = strchr(struc, '&')-struc;
 195   int size;
 196   size=strlen(struc)-1;
 197   Emin-=size * n_seq * extension_cost;
 198   mfe.i = i_min;
 199   mfe.j = j_min;
 200   mfe.energy = (float) (Emin/(100.*n_seq));
 201   mfe.structure = struc;
 202   if (!delay_free) {
 203     for (i=1; i<=n3; i++) free(c[i]);
 204     free(c);
 205   }
 206   for (s=0; s<n_seq; s++) {
 207     free(S1[s]); free(S2[s]);
 208   }
 209   free(S1); free(S2); free(type);
 210   return mfe;
 211 }
 212
 213
 214 PRIVATE char *alibacktrack(int i, int j, const short *S1[], const short *S2[], const int extension_cost) {
 215   /* backtrack structure going backwards from i, and forwards from j
 216      return structure in bracket notation with & as separator */
 217   int k, l, *type, type2, E, traced, i0, j0, s, n_seq;
 218   char *st1, *st2, *struc;
 219
 220   n3 = (int) S1[0][0];
 221   n4 = (int) S2[0][0];
 222
 223   for (s=0; S1[s]!=NULL; s++);
 224   n_seq = s;
 225   for (s=0; S2[s]!=NULL; s++);
 226   if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");
 227
 228   st1 = (char *) space(sizeof(char)*(n3+1));
 229   st2 = (char *) space(sizeof(char)*(n4+1));
 230   type = (int *) space(n_seq*sizeof(int));
 231
 232   i0=MIN2(i+1,n3); j0=MAX2(j-1,1);
 233
 234   while (i>0 && j<=n4) {
 235     int psc;
 236     E = c[i][j]; traced=0;
 237     st1[i-1] = '(';
 238     st2[j-1] = ')';
 239     for (s=0; s<n_seq; s++) {
 240       type[s] = pair[S1[s][i]][S2[s][j]];
 241     }
 242     psc = covscore(type, n_seq);
 243     for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
 244     E += psc;
 245     for (k=i-1; k>0 && k>i-MAXLOOP-2; k--) {
 246       for (l=j+1; l<=n4; l++) {
 247         int LE;
 248         if (i-k+l-j-2>MAXLOOP) break;
 249         if (c[k][l]>INF/2) continue;
 250         for (s=LE=0; s<n_seq; s++) {
 251           type2 = pair[S1[s][k]][S2[s][l]];
 252           if (type2==0) type2=7;
 253           LE += E_IntLoop(i-k-1, l-j-1, type2, rtype[type[s]],
 254                           S1[s][k+1], S2[s][l-1], S1[s][i-1], S2[s][j+1],P)+(i-k+l-j)*extension_cost;
 255         }
 256         if (E == c[k][l]+LE) {
 257           traced=1;
 258           i=k; j=l;
 259           break;
 260         }
 261       }
 262       if (traced) break;
 263     }
 264     if (!traced) {
 265       for (s=0; s<n_seq; s++) {
 266         E -= E_ExtLoop(type[s], (i>1) ? S1[s][i-1] : -1, (j<n4) ? S2[s][j+1] : -1, P) + 2*extension_cost;
 267       }
 268       if (E != n_seq*P->DuplexInit + n_seq*2*extension_cost) {
 269         nrerror("backtrack failed in aliduplex");
 270       } else break;
 271     }
 272   }
 273   if (i>1)  i--;
 274   if (j<n4) j++;
 275
 276   struc = (char *) space(i0-i+1+j-j0+1+2);
 277   for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
 278   for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = '.';
 279   strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&");
 280   strcat(struc, st2+j0-1);
 281
 282   /* printf("%s %3d,%-3d : %3d,%-3d\n", struc, i,i0,j0,j);  */
 283   free(st1); free(st2); free(type);
 284
 285   return struc;
 286 }
 287
 288 duplexT** aliLduplexfold(const char *s1[], const char *s2[], const int threshold, const int extension_cost, const int alignment_length, const int delta, const int fast,const int il_a, const int il_b, const int b_a, const int b_b)
 289 {
 290   short **S1, **S2;
 291   int *type, type2;
 292   int i, j,s,n_seq;
 293   s=0;
 294   int bopen=b_b;
 295   int bext=b_a+extension_cost;
 296   int iopen=il_b;
 297   int iext_s=2*(il_a+extension_cost);/* iext_s 2 nt nucleotide extension of interior loop, on i and j side */
 298   int iext_ass=50+il_a+extension_cost;/* iext_ass assymetric extension of interior loop, either on i or on j side. */
 299   int min_colonne=INF; /* enthaelt das maximum einer kolonne */
 300   int i_length;
 301   int max_pos;/* get position of the best hit */
 302   int max_pos_j;
 303   int temp;
 304   int min_j_colonne;
 305   int max=INF;
 306   /* FOLLOWING NEXT 4 LINE DEFINES AN ARRAY CONTAINING POSITION OF THE SUBOPT IN S1 */
 307   int *position; /* contains the position of the hits with energy > E */
 308   int *position_j;
 309
 310
 311   n1 = (int) strlen(s1[0]);
 312   n2 = (int) strlen(s2[0]);
 313   for (s=0; s1[s]; s++);
 314   n_seq = s;
 315   for (s=0; s2[s]; s++);
 316   if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
 317
 318   position = (int *) space ((delta+(n1)+4+delta) * sizeof(int));
 319   position_j= (int *) space ((delta+(n1)+4+delta) * sizeof(int));
 320
 321   if ((!P) || (fabs(P->temperature - temperature)>1e-6)){
 322     update_dfold_params();
 323   }
 324
 325   lc   = (int**) space (sizeof(int *) * 5);
 326   lin  = (int**) space (sizeof(int *) * 5);
 327   lbx  = (int**) space (sizeof(int *) * 5);
 328   lby  = (int**) space (sizeof(int *) * 5);
 329   linx = (int**) space (sizeof(int *) * 5);
 330   liny = (int**) space (sizeof(int *) * 5);
 331
 332   for (i=0; i<=4; i++){
 333     lc[i]  = (int *) space(sizeof(int) * (n2+5));
 334     lin[i] = (int *) space(sizeof(int) * (n2+5));
 335     lbx[i] = (int *) space(sizeof(int) * (n2+5));
 336     lby[i] = (int *) space(sizeof(int) * (n2+5));
 337     linx[i]= (int *) space(sizeof(int) * (n2+5));
 338     liny[i]= (int *) space(sizeof(int) * (n2+5));
 339   }
 340
 341
 342   S1 = (short **) space((n_seq+1)*sizeof(short *));
 343   S2 = (short **) space((n_seq+1)*sizeof(short *));
 344   for (s=0; s<n_seq; s++) {
 345     if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
 346     if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
 347     S1[s] = encode_seq(s1[s]);
 348     S2[s] = encode_seq(s2[s]);
 349   }
 350   type = (int *) space(n_seq*sizeof(int));
 351   /**
 352   *** array initialization
 353   **/
 354   for(j=n2;j>=0;j--) {
 355     lbx[0][j]=lbx[1][j]=lbx[2][j]=lbx[3][j]    = lbx[4][j] =INF;
 356     lin[0][j]=lin[1][j]=lin[2][j]=lin[3][j]    = lin[4][j] =INF;
 357     lc[0][j] =lc[1][j] =lc[2][j] = lc[3][j]    =  lc[4][j] =INF;
 358     lby[0][j]=lby[1][j]=lby[2][j]=lby[3][j]    = lby[4][j] =INF;
 359     liny[0][j]=liny[1][j]=liny[2][j]=liny[3][j]=liny[4][j]=INF;
 360     linx[0][j]=linx[1][j]=linx[2][j]=linx[3][j]=linx[4][j]=INF;
 361   }
 362   i=10;
 363   i_length= n1 - 9 ;
 364   while(i < i_length) {
 365     int idx=i%5;
 366     int idx_1=(i-1)%5;
 367     int idx_2=(i-2)%5;
 368     int idx_3=(i-3)%5;
 369     int idx_4=(i-4)%5;
 370     j=n2-9;
 371     while (9 < --j) {
 372       int psc;
 373       for (s=0; s<n_seq; s++) {
 374         type[s] = pair[S1[s][i]][S2[s][j]];
 375       }
 376       psc = covscore(type, n_seq);
 377       for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
 378       lc[idx][j] = (psc>=MINPSCORE) ? (n_seq*P->DuplexInit + 2*n_seq*extension_cost) : INF;
 379       /**
 380       *** Update matrix. It is the average over all sequence of a given structure element
 381       *** c_stack -> stacking of c
 382       *** c_10, c01 -> stack from bulge
 383       *** c_nm -> arrives in stack from nxm loop
 384       *** c_in -> arrives in stack from interior loop
 385       *** c_bx -> arrives in stack from large bulge on target
 386       *** c_by -> arrives in stack from large bulge on query
 387       ***
 388       **/
 389       int c_stack, c_10, c_01, c_11, c_22, c_21, c_12, c_23, c_32, c_in, c_in2x, c_in2y, c_bx, c_by, c_inx, c_iny;  /* matrix c */
 390       int in, in_x, in_y, in_xy; /*  in begin, in_x assymetric, in_y assymetric, in_xy symetric; */
 391       int inx, inx_x;
 392       int iny, iny_y;
 393       int bx, bx_x;
 394       int by, by_y;
 395       in=lc[idx_1][j+1]; in_x=lin[idx_1][j]; in_y=lin[idx][j+1]; in_xy=lin[idx_1][j+1];
 396       inx=lc[idx_1][j+1]; inx_x=linx[idx_1][j];
 397       iny=lc[idx_1][j+1]; iny_y=liny[idx][j+1];
 398       bx=lc[idx_1][j]; bx_x=lbx[idx_1][j];
 399       by=lc[idx][j+1]; by_y=lby[idx][j+1];
 400       c_stack=lc[idx_1][j+1]; c_01=lc[idx_1][j+2];c_10=lc[idx_2][j+1];
 401       c_12=lc[idx_2][j+3];c_21=lc[idx_3][j+2];c_11=lc[idx_2][j+2];
 402       c_22=lc[idx_3][j+3];c_32=lc[idx_4][j+3];c_23=lc[idx_3][j+4];
 403       c_in=lin[idx_3][j+3];c_in2x=lin[idx_4][j+2];c_in2y=lin[idx_2][j+4];
 404       c_inx=linx[idx_3][j+1]; c_iny=liny[idx_1][j+3];
 405       c_bx=lbx[idx_2][j+1];c_by=lby[idx_1][j+2];
 406       for (s=0; s<n_seq; s++) {
 407         type2 = pair[S2[s][j+1]][S1[s][i-1]];
 408         in   +=P->mismatchI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s;
 409         in_x +=iext_ass;
 410         in_y +=iext_ass;
 411         in_xy+=iext_s;
 412         inx  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s;
 413         inx_x+=iext_ass;
 414         iny  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s;
 415         iny_y+=iext_ass;
 416         type2=pair[S2[s][j]][S1[s][i-1]];
 417         bx   +=bopen+bext+(type2>2?P->TerminalAU:0);
 418         bx_x +=bext;
 419         type2=pair[S2[s][j+1]][S1[s][i]];
 420         by   +=bopen+bext+(type2>2?P->TerminalAU:0);
 421         by_y +=bext;
 422       }
 423       lin [idx][j]=MIN2(in, MIN2(in_x, MIN2(in_y, in_xy)));
 424       linx[idx][j]=MIN2(inx_x, inx);
 425       liny[idx][j]=MIN2(iny_y, iny);
 426       lby[idx][j] =MIN2(by, by_y);
 427       lbx[idx][j] =MIN2(bx, bx_x);
 428
 429       if (psc<MINPSCORE) continue;
 430       for (s=0; s<n_seq; s++) {
 431         lc[idx][j]+=E_ExtLoop(type[s], S1[s][i-1],S2[s][j+1], P) + 2*extension_cost;
 432       }
 433       for (s=0; s<n_seq; s++) {
 434         type2=pair[S1[s][i-1]][S2[s][j+1]];if (type2==0) type2=7;
 435         c_stack+=E_IntLoop(0,0,type2, rtype[type[s]],S1[s][i], S2[s][j], S1[s][i-1], S2[s][j+1], P)+2*extension_cost;
 436         type2=pair[S1[s][i-1]][S2[s][j+2]];if (type2==0) type2=7;
 437         c_01   +=E_IntLoop(0,1,type2, rtype[type[s]],S1[s][i], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+3*extension_cost;
 438         type2=pair[S1[s][i-2]][S2[s][j+1]]; if (type2==0) type2=7;
 439         c_10   +=E_IntLoop(1,0,type2, rtype[type[s]],S1[s][i-1], S2[s][j], S1[s][i-1], S2[s][j+1], P)+3*extension_cost;
 440         type2=pair[S1[s][i-2]][S2[s][j+2]]; if (type2==0) type2=7;
 441         c_11   +=E_IntLoop(1,1,type2, rtype[type[s]],S1[s][i-1], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+4*extension_cost;
 442         type2 = pair[S1[s][i-3]][S2[s][j+3]];if (type2==0) type2=7;
 443         c_22   +=E_IntLoop(2,2,type2, rtype[type[s]],S1[s][i-2], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+6*extension_cost;
 444         type2 = pair[S1[s][i-3]][S2[s][j+2]];if (type2==0) type2=7;
 445         c_21   +=E_IntLoop(2,1,type2, rtype[type[s]],S1[s][i-2], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+5*extension_cost;
 446         type2 = pair[S1[s][i-2]][S2[s][j+3]];if (type2==0) type2=7;
 447         c_12   +=E_IntLoop(1,2,type2, rtype[type[s]],S1[s][i-1], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+5*extension_cost;
 448         type2 = pair[S1[s][i-4]][S2[s][j+3]];if (type2==0) type2=7;
 449         c_32   +=E_IntLoop(3,2,type2, rtype[type[s]],S1[s][i-3], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+7*extension_cost;
 450         type2 = pair[S1[s][i-3]][S2[s][j+4]];if (type2==0) type2=7;
 451         c_23   +=E_IntLoop(2,3,type2, rtype[type[s]],S1[s][i-2], S2[s][j+3], S1[s][i-1], S2[s][j+1], P)+7*extension_cost;
 452         c_in   +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+2*extension_cost+2*iext_s;
 453         c_in2x +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+2*extension_cost;
 454         c_in2y +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+2*extension_cost;
 455         c_inx  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+2*extension_cost;
 456         c_iny  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+2*extension_cost;
 457         int bAU;
 458         bAU=(type[s]>2?P->TerminalAU:0);
 459         c_bx   +=2*extension_cost+bext+bAU;
 460         c_by   +=2*extension_cost+bext+bAU;
 461       }
 462       lc[idx][j] =MIN2(lc[idx][j],
 463                        MIN2(c_stack,
 464                             MIN2(c_10,
 465                                  MIN2(c_01,
 466                                       MIN2(c_11,
 467                                            MIN2(c_21,
 468                                                 MIN2(c_12,
 469                                                      MIN2(c_22,
 470                                                           MIN2(c_23,
 471                                                                MIN2(c_32,
 472                                                                     MIN2(c_bx,
 473                                                                          MIN2(c_by,
 474                                                                               MIN2(c_in,
 475                                                                                    MIN2(c_in2x,
 476                                                                                         MIN2(c_in2y,
 477                                                                                              MIN2(c_inx,c_iny)
 478                                                                                              )
 479                                                                                         )
 480                                                                                    )
 481                                                                               )
 482                                                                          )
 483                                                                     )
 484                                                                )
 485                                                           )
 486                                                      )
 487                                                 )
 488                                            )
 489                                       )
 490                                  )
 491                             )
 492                        );
 493       lc[idx][j]-=psc;
 494       temp=lc[idx][j];
 495       for (s=0; s<n_seq; s++) {
 496         temp+=E_ExtLoop(rtype[type[s]], S2[s][j-1],S1[s][i+1],P)+2*extension_cost;
 497       }
 498       if(min_colonne > temp){
 499         min_colonne=temp;
 500         min_j_colonne=j;
 501       }
 502     }
 503     if(max>=min_colonne){
 504             max=min_colonne;
 505             max_pos=i;
 506         max_pos_j=min_j_colonne;
 507     }
 508     position[i+delta]=min_colonne;min_colonne=INF;
 509     position_j[i+delta]=min_j_colonne;
 510     i++;
 511   }
 512   /* printf("MAX:%d ",max); */
 513   for (s=0; s<n_seq; s++) {free(S1[s]);free(S2[s]);}
 514   free(S1); free(S2);
 515   if(max<threshold){
 516     alifind_max(position, position_j, delta, threshold, alignment_length, s1, s2, extension_cost, fast);
 517   }
 518   aliplot_max(max, max_pos, max_pos_j,alignment_length, s1, s2, extension_cost,fast);
 519   for (i=0; i<=4; i++) {free(lc[i]);free(lin[i]);free(lbx[i]);free(lby[i]);free(linx[i]);free(liny[i]);}
 520   /* free(lc[0]);free(lin[0]);free(lbx[0]);free(lby[0]);free(linx[0]);free(liny[0]); */
 521   free(lc);free(lin);free(lbx);free(lby);free(linx);free(liny);
 522   free(position);
 523   free(position_j);
 524   free(type);
 525   return NULL;
 526 }
 527
 528
 529 PRIVATE void alifind_max(const int *position, const int *position_j,
 530                          const int delta, const int threshold, const int alignment_length,
 531                          const char *s1[], const char *s2[],
 532                          const int extension_cost, const int fast){
 533   int n_seq=0;
 534   for (n_seq=0; s1[n_seq]!=NULL; n_seq++);
 535   int pos=n1-9;
 536   if(fast==1){
 537     while(10<pos--){
 538       int temp_min=0;
 539       if(position[pos+delta]<(threshold)){
 540         int search_range;
 541         search_range=delta+1;
 542         while(--search_range){
 543           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
 544             temp_min=search_range;
 545           }
 546         }
 547         pos-=temp_min;
 548         int max_pos_j;
 549         max_pos_j=position_j[pos+delta];
 550         int max;
 551         max=position[pos+delta];
 552         printf("target upper bound %d: query lower bound %d  (%5.2f) \n", pos-10, max_pos_j-10, ((double)max)/(n_seq*100));
 553         pos=MAX2(10,pos-delta);
 554       }
 555     }
 556   }
 557   else{
 558     pos=n1-9;
 559     while(pos-- > 10) {
 560       /* printf("delta %d position:%d value:%d\n", delta, pos, position[pos]); */
 561       int temp_min=0;
 562       if(position[pos+delta]<(threshold)){
 563         int search_range;
 564         search_range=delta+1;
 565         while(--search_range){
 566           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
 567             temp_min=search_range;
 568           }
 569         }
 570         pos-=temp_min;
 571         int max_pos_j;
 572         max_pos_j=position_j[pos+delta];
 573         /* printf("%d %d %d\n", pos, max_pos_j,position[pos+delta]); */
 574         int begin_t=MAX2(11, pos-alignment_length+1);
 575         int end_t  =MIN2(n1-10, pos+1);
 576         int begin_q=MAX2(11, max_pos_j-1);
 577         int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
 578         char **s3, **s4;
 579         s3 = (char**) space(sizeof(char*)*(n_seq+1));
 580         s4 = (char**) space(sizeof(char*)*(n_seq+1));
 581         int i;
 582         for(i=0; i<n_seq; i++){
 583           s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
 584           s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
 585           strncpy(s3[i], (s1[i]+begin_t-1), end_t - begin_t +1);
 586           strncpy(s4[i], (s2[i]+begin_q-1), end_q - begin_q +1);
 587           s3[i][end_t - begin_t +1]='\0';
 588           s4[i][end_q - begin_q +1]='\0';
 589         }
 590         duplexT test;
 591         test = aliduplexfold((const char**)s3, (const char**)s4, extension_cost);
 592         /* printf("test %d threshold %d",test.energy*100,(threshold/n_seq)); */
 593         if(test.energy * 100 < (int) (threshold/n_seq)){
 594           int l1=strchr(test.structure, '&')-test.structure;
 595                    printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", test.structure,
 596                  begin_t -10 +test.i-l1,
 597                  begin_t -10 +test.i-1,
 598                  begin_q -10 + test.j - 1,
 599                  begin_q-11 + test.j +strlen(test.structure) -l1 -2 , test.energy);
 600           pos=MAX2(10,pos-delta);
 601
 602         }
 603         for(i=0;i<n_seq;i++){
 604           free(s3[i]);free(s4[i]);
 605         }
 606         free(s3);free(s4);
 607         free(test.structure);
 608       }
 609     }
 610   }
 611 }
 612 PRIVATE void aliplot_max(const int max, const int max_pos, const int max_pos_j, const int alignment_length, const char *s1[], const char *s2[], const int extension_cost, const int fast)
 613 {
 614   int n_seq;
 615   for (n_seq=0; !(s1[n_seq]==NULL); n_seq++);
 616   n1 = strlen(s1[0]); /* get length of alignment */
 617   n2 = strlen(s2[0]); /* get length of alignment */
 618   if(fast==1){
 619     printf("target upper bound %d: query lower bound %d (%5.2f)\n",
 620            max_pos-10, max_pos_j-10, (double) ((double)max)/(100*n_seq));
 621   }
 622   else{
 623     int begin_t=MAX2(11, max_pos-alignment_length+1);
 624     int end_t  =MIN2(n1-10, max_pos+1);
 625     int begin_q=MAX2(11, max_pos_j-1);
 626     int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
 627     char **s3, **s4;
 628     s3 = (char**) space(sizeof(char*)*(n_seq+1));
 629     s4 = (char**) space(sizeof(char*)*(n_seq+1));
 630     int i;
 631     for(i=0; i<n_seq; i++){
 632       s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
 633       s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
 634       strncpy(s3[i], (s1[i]+begin_t-1), end_t - begin_t +1);
 635       strncpy(s4[i], (s2[i]+begin_q-1), end_q - begin_q +1);
 636       s3[i][end_t - begin_t +1]='\0';
 637       s4[i][end_q - begin_q +1]='\0';
 638     }
 639     duplexT test;
 640     s3[n_seq]=s4[n_seq]=NULL;
 641     test = aliduplexfold((const char**) s3,(const char**) s4, extension_cost);
 642     int l1=strchr(test.structure, '&')-test.structure;
 643     printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n",
 644            test.structure,
 645            begin_t -10 +test.i-l1,
 646            begin_t -10 +test.i-1,
 647            begin_q-10 + test.j - 1,
 648            begin_q -11 + test.j +strlen(test.structure) - l1 - 2,
 649            test.energy);
 650     for(i=0; i<n_seq ; i++){
 651       free(s3[i]);free(s4[i]);
 652     }
 653     free(s3);free(s4);
 654     free(test.structure);
 655   }
 656 }
 657
 658 PRIVATE duplexT aliduplexfold_XS(const char *s1[], const char *s2[],
 659                                  const int **access_s1, const int **access_s2,
 660                                  const int i_pos, const int j_pos, const int threshold,
 661                                  const int i_flag, const int j_flag){
 662   int i,j,s,p,q, Emin=INF, l_min=0, k_min=0;
 663   char *struc;
 664   short **S1,**S2;
 665   int *type,*type2;
 666   struc=NULL;
 667   duplexT mfe;
 668   int n_seq;
 669   n3 = (int) strlen(s1[0]);
 670   n4 = (int) strlen(s2[0]);
 671   for (s=0; s1[s]!=NULL; s++);
 672   n_seq = s;
 673   for (s=0; s2[s]!=NULL; s++);
 674   /* printf("%d \n",i_pos); */
 675   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
 676     update_fold_params();  if(P) free(P); P = scale_parameters();
 677     make_pair_matrix();
 678   }
 679   c = (int **) space(sizeof(int *) * (n3+1));
 680   for (i=0; i<=n3; i++) c[i] = (int *) space(sizeof(int) * (n4+1));
 681   for (i=0; i<=n3; i++){
 682     for(j=0;j<=n4;j++){
 683       c[i][j]=INF;
 684     }
 685   }
 686   S1 = (short **) space((n_seq+1)*sizeof(short *));
 687   S2 = (short **) space((n_seq+1)*sizeof(short *));
 688   for (s=0; s<n_seq; s++) {
 689     if (strlen(s1[s]) != n3) nrerror("uneqal seqence lengths");
 690     if (strlen(s2[s]) != n4) nrerror("uneqal seqence lengths");
 691     S1[s] = encode_seq(s1[s]);
 692     S2[s] = encode_seq(s2[s]);
 693   }
 694   type =  (int *) space(n_seq*sizeof(int));
 695   type2 = (int *) space(n_seq*sizeof(int));
 696   int type3, E, k,l;
 697   i=n3-i_flag; j=1+j_flag;
 698   for (s=0; s<n_seq; s++) {
 699     type[s] = pair[S1[s][i]][S2[s][j]];
 700   }
 701   c[i][j] = n_seq*P->DuplexInit - covscore(type,n_seq);
 702   for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
 703   for (s=0; s<n_seq; s++) {
 704     c[i][j]+=E_ExtLoop(rtype[type[s]], (j_flag ? S2[s][j-1] : -1) , (i_flag ? S1[s][i+1] : -1),  P);
 705   }
 706   k_min=i; l_min=j; Emin=c[i][j];
 707   for (k=i; k>1 ; k--) {
 708     if(k<i) c[k+1][0]=INF;
 709     for (l=j; l<=n4-1; l++) {
 710       if(!(k==i && l==j)){
 711         c[k][l]=INF;
 712       }
 713       int psc2;
 714       for(s=0;s<n_seq;s++){
 715         type2[s] = pair[S1[s][k]][S2[s][l]];
 716       }
 717       psc2=covscore(type2, n_seq);
 718       if (psc2<MINPSCORE) continue;
 719       for (s=0; s<n_seq; s++) if (type2[s]==0) type2[s]=7;
 720       for (p=k+1; p<= n3 -i_flag && p<k+MAXLOOP-1; p++) {
 721         for (q = l-1; q >= 1+j_flag; q--) {
 722           if (p-k+l-q-2>MAXLOOP) break;
 723           for(E=s=0;s<n_seq;s++){
 724             type3=pair[S1[s][p]][S2[s][q]];
 725             if(type3==0) type3=7;
 726             E += E_IntLoop(p-k-1, l-q-1, type2[s], rtype[type3],
 727                            S1[s][k+1], S2[s][l-1], S1[s][p-1], S2[s][q+1],P);
 728           }
 729           c[k][l] = MIN2(c[k][l], c[p][q]+E);
 730         }
 731       }
 732       c[k][l]-=psc2;
 733       E = c[k][l];
 734       E+=n_seq*(access_s1[i-k+1][i_pos]+access_s2[l-1][j_pos+(l-1)-1]);
 735       for (s=0; s<n_seq; s++) {
 736         E+=E_ExtLoop(type2[s], (k>1) ? S1[s][k-1] : -1, (l<n4) ? S2[s][l+1] : -1, P);
 737       }
 738       if (E<Emin) {
 739         Emin=E; k_min=k; l_min=l;
 740       }
 741     }
 742   }
 743   if(Emin > threshold-1){
 744     mfe.structure=NULL;
 745     mfe.energy=INF;
 746     for (i=0; i<=n3; i++) free(c[i]);
 747     free(c);
 748     for(i=0; i<=n_seq;i++){
 749       free(S1[i]);
 750       free(S2[i]);
 751     }
 752     free(S1); free(S2); /* free(SS1); free(SS2); */
 753     free(type);free(type2);
 754     return mfe;
 755     } else{
 756     struc = alibacktrack_XS(k_min, l_min,(const short int**)S1,(const short int**)S2,access_s1, access_s2,i_flag,j_flag);
 757     }
 758   int dx_5, dx_3, dy_5, dy_3,dGx,dGy,bonus_x, bonus_y,temp_dangle;
 759   dx_5=0; dx_3=0; dy_5=0; dy_3=0;dGx=0;dGy=0;bonus_x=0; bonus_y=0;temp_dangle=0;
 760   dGx =n_seq*(access_s1[i-k_min+1][i_pos]);dx_3=0; dx_5=0;bonus_x=0;
 761   dGy =n_seq*access_s2[l_min-j+1][j_pos + (l_min-1)-1];
 762   mfe.tb=i_pos -9 - i + k_min -1 -dx_5;
 763   mfe.te=i_pos -9 -1 + dx_3;
 764   mfe.qb=j_pos -9 -1 - dy_5;
 765   mfe.qe=j_pos + l_min -3 -9 + dy_3;
 766   mfe.ddG=(double) (Emin *0.01);
 767   mfe.dG1=(double) (dGx*(0.01));
 768   mfe.dG2=(double) (dGy*(0.01));
 769   mfe.energy= mfe.ddG - mfe.dG1 - mfe.dG2;
 770   mfe.structure = struc;
 771   for (i=0; i<=n3; i++) free(c[i]);
 772   free(c);
 773   for(i=0; i<=n_seq;i++){
 774     free(S1[i]);
 775     free(S2[i]);
 776   }
 777   free(S1); free(S2); free(type);free(type2);
 778   return mfe;
 779 }
 780
 781 PRIVATE char *alibacktrack_XS(int i, int j, const short *S1[], const short *S2[],
 782                               const int** access_s1, const int ** access_s2, const int i_flag, const int j_flag) {
 783   int n3,n4,k, l, *type, type2, E, traced, i0, j0,s,n_seq,psc;
 784   char *st1=NULL, *st2=NULL, *struc=NULL;
 785
 786   n3 = (int) S1[0][0];
 787   n4 = (int) S2[0][0];
 788   for (s=0; S1[s]!=NULL; s++);
 789   n_seq = s;
 790   for (s=0; S2[s]!=NULL; s++);
 791   if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");
 792
 793   st1 = (char *) space(sizeof(char)*(n3+1));
 794   st2 = (char *) space(sizeof(char)*(n4+1));
 795   type = (int *) space(n_seq*sizeof(int));
 796
 797   i0=i;/*MAX2(i-1,1);*/j0=j;/*MIN2(j+1,n4);*/
 798   while (i<=n3-i_flag && j>=1+j_flag) {
 799     E = c[i][j]; traced=0;
 800     st1[i-1] = '(';
 801     st2[j-1] = ')';
 802     for (s=0; s<n_seq; s++) {
 803       type[s] = pair[S1[s][i]][S2[s][j]];
 804     }
 805     psc = covscore(type,n_seq);
 806     for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
 807     E += psc;
 808     for (k=i+1; k<=n3 && k>i-MAXLOOP-2; k++) {
 809       for (l=j-1; l>=1; l--) {
 810         int LE;
 811         if (i-k+l-j-2>MAXLOOP) break;
 812         for (s=LE=0; s<n_seq; s++) {
 813           type2 = pair[S1[s][k]][S2[s][l]];
 814           if (type2==0) type2=7;
 815           LE += E_IntLoop(k-i-1, j-l-1, type[s], rtype[type2],
 816                           S1[s][i+1], S2[s][j-1], S1[s][k-1], S2[s][l+1],P);
 817         }
 818         if (E == c[k][l]+LE) {
 819           traced=1;
 820           i=k; j=l;
 821           break;
 822         }
 823       }
 824       if (traced) break;
 825     }
 826     if (!traced) {
 827       for (s=0; s<n_seq; s++) {
 828         if (type[s]>2) E -= P->TerminalAU;
 829       }
 830       break;
 831       if (E != n_seq*P->DuplexInit) {
 832         nrerror("backtrack failed in fold duplex bal");
 833       } else break;
 834     }
 835   }
 836   struc = (char *) space(i-i0+1+j0-j+1+2);
 837   for (k=MAX2(i0,1); k<=i; k++) if (!st1[k-1]) st1[k-1] = '.';
 838   for (k=j; k<=j0; k++) if (!st2[k-1]) st2[k-1] = '.';
 839   strcpy(struc, st1+MAX2(i0-1,0)); strcat(struc, "&");
 840   strcat(struc, st2+j-1);
 841   free(st1);
 842   free(st2);
 843   free(type);
 844   return struc;
 845 }
 846
 847 duplexT** aliLduplexfold_XS(const char*s1[], const char* s2[], const int **access_s1, const int **access_s2, const int threshold, const int alignment_length, const int delta,  const int fast,const int il_a, const int il_b, const int b_a, const int b_b)
 848 {
 849   short **S1, **S2;
 850   int *type,type2;
 851   int i, j,s,n_seq;
 852   s=0;
 853   int bopen=b_b;
 854   int bext=b_a;
 855   int iopen=il_b;
 856   int iext_s=2*(il_a);/* iext_s 2 nt nucleotide extension of interior loop, on i and j side */
 857   int iext_ass=50+il_a;/* iext_ass assymetric extension of interior loop, either on i or on j side. */
 858   int min_colonne=INF; /* enthaelt das maximum einer kolonne */
 859   int i_length;
 860   int max_pos;/* get position of the best hit */
 861   int max_pos_j;
 862   int temp;
 863   int min_j_colonne;
 864   int max=INF;
 865   /* FOLLOWING NEXT 4 LINE DEFINES AN ARRAY CONTAINING POSITION OF THE SUBOPT IN S1 */
 866   int *position; /* contains the position of the hits with energy > E */
 867   int *position_j;
 868   int maxPenalty[4];
 869
 870   n1 = (int) strlen(s1[0]);
 871   n2 = (int) strlen(s2[0]);
 872   for (s=0; s1[s]; s++);
 873   n_seq = s;
 874   for (s=0; s2[s]; s++);
 875   if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
 876
 877   position = (int *) space ((delta+(n1)+4+delta) * sizeof(int));
 878   position_j= (int *) space ((delta+(n1)+4+delta) * sizeof(int));
 879
 880   /**
 881   *** extension penalty, computed only once, further reduce the computation time
 882   **/
 883   if ((!P) || (fabs(P->temperature - temperature)>1e-6)){
 884     update_dfold_params();
 885   }
 886   maxPenalty[0]=(int) -1*P->stack[2][2]/2;
 887   maxPenalty[1]=(int) -1*P->stack[2][2];
 888   maxPenalty[2]=(int) -3*P->stack[2][2]/2;
 889   maxPenalty[3]=(int) -2*P->stack[2][2];
 890
 891
 892   lc   = (int**) space (sizeof(int *) * 5);
 893   lin  = (int**) space (sizeof(int *) * 5);
 894   lbx  = (int**) space (sizeof(int *) * 5);
 895   lby  = (int**) space (sizeof(int *) * 5);
 896   linx = (int**) space (sizeof(int *) * 5);
 897   liny = (int**) space (sizeof(int *) * 5);
 898
 899   for (i=0; i<=4; i++){
 900     lc[i]  = (int *) space(sizeof(int) * (n2+5));
 901     lin[i] = (int *) space(sizeof(int) * (n2+5));
 902     lbx[i] = (int *) space(sizeof(int) * (n2+5));
 903     lby[i] = (int *) space(sizeof(int) * (n2+5));
 904     linx[i]= (int *) space(sizeof(int) * (n2+5));
 905     liny[i]= (int *) space(sizeof(int) * (n2+5));
 906   }
 907
 908
 909   S1 = (short **) space((n_seq+1)*sizeof(short *));
 910   S2 = (short **) space((n_seq+1)*sizeof(short *));
 911   for (s=0; s<n_seq; s++) {
 912     if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
 913     if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
 914     S1[s] = encode_seq(s1[s]);
 915     S2[s] = encode_seq(s2[s]);
 916   }
 917   type = (int *) space(n_seq*sizeof(int));
 918   /**
 919   *** array initialization
 920   **/
 921
 922   for(j=n2+4;j>=0;j--) {
 923     lbx[0][j]=lbx[1][j]=lbx[2][j]=lbx[3][j]    = lbx[4][j] =INF;
 924     lin[0][j]=lin[1][j]=lin[2][j]=lin[3][j]    = lin[4][j] =INF;
 925     lc[0][j] =lc[1][j] =lc[2][j] = lc[3][j]    =  lc[4][j] =INF;
 926     lby[0][j]=lby[1][j]=lby[2][j]=lby[3][j]    = lby[4][j] =INF;
 927     liny[0][j]=liny[1][j]=liny[2][j]=liny[3][j]=liny[4][j]=INF;
 928     linx[0][j]=linx[1][j]=linx[2][j]=linx[3][j]=linx[4][j]=INF;
 929   }
 930   i=10;
 931   i_length= n1 - 9 ;
 932   int di1,di2,di3,di4; /* contains accessibility penalty */
 933   while(i < i_length) {
 934     int idx=i%5;
 935     int idx_1=(i-1)%5;
 936     int idx_2=(i-2)%5;
 937     int idx_3=(i-3)%5;
 938     int idx_4=(i-4)%5;
 939
 940     di1 = access_s1[5][i]   - access_s1[4][i-1];
 941     di2 = access_s1[5][i-1] - access_s1[4][i-2] + di1;
 942     di3 = access_s1[5][i-2] - access_s1[4][i-3] + di2;
 943     di4 = access_s1[5][i-3] - access_s1[4][i-4] + di3;
 944     di1=MIN2(di1,maxPenalty[0]);
 945     di2=MIN2(di2,maxPenalty[1]);
 946     di3=MIN2(di3,maxPenalty[2]);
 947     di4=MIN2(di3,maxPenalty[3]);
 948     j=n2-9;
 949     while (9 < --j) {
 950       int dj1,dj2,dj3,dj4;
 951       dj1 = access_s2[5][j+4] - access_s2[4][j+4];
 952       dj2 = access_s2[5][j+5] - access_s2[4][j+5] + dj1;
 953       dj3 = access_s2[5][j+6] - access_s2[4][j+6] + dj2;
 954       dj4 = access_s2[5][j+7] - access_s2[4][j+7] + dj3;
 955       dj1=MIN2(dj1,maxPenalty[0]);
 956       dj2=MIN2(dj2,maxPenalty[1]);
 957       dj3=MIN2(dj3,maxPenalty[2]);
 958       dj4=MIN2(dj4,maxPenalty[3]);
 959       int psc;
 960       for (s=0; s<n_seq; s++) { /* initialize type1 */
 961         type[s] = pair[S1[s][i]][S2[s][j]];
 962       }
 963       psc = covscore(type, n_seq);
 964       for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
 965       lc[idx][j] = ((psc >= MINPSCORE) ? n_seq*(P->DuplexInit + access_s1[1][i] + access_s2[1][j]) : INF);
 966       int c_stack, c_10, c_01, c_11, c_22, c_21, c_12, c_23, c_32, c_in, c_in2x, c_in2y, c_bx, c_by, c_inx, c_iny;  /* matrix c */
 967       int in,  in_x, in_y, in_xy; /*  in begin, in_x assymetric, in_y assymetric, in_xy symetric; */
 968       int inx, inx_x;
 969       int iny, iny_y;
 970       int bx,  bx_x;
 971       int by,  by_y;
 972       in=lc[idx_1][j+1]; in_x=lin[idx_1][j]; in_y=lin[idx][j+1]; in_xy=lin[idx_1][j+1];
 973       inx=lc[idx_1][j+1]; inx_x=linx[idx_1][j];
 974       iny=lc[idx_1][j+1]; iny_y=liny[idx][j+1];
 975       bx=lc[idx_1][j]; bx_x=lbx[idx_1][j];
 976       by=lc[idx][j+1]; by_y=lby[idx][j+1];
 977       c_stack=lc[idx_1][j+1]; c_01=lc[idx_1][j+2];c_10=lc[idx_2][j+1];
 978       c_12=lc[idx_2][j+3];c_21=lc[idx_3][j+2];c_11=lc[idx_2][j+2];
 979       c_22=lc[idx_3][j+3];c_32=lc[idx_4][j+3];c_23=lc[idx_3][j+4];
 980       c_in=lin[idx_3][j+3];c_in2x=lin[idx_4][j+2];c_in2y=lin[idx_2][j+4];
 981       c_inx=linx[idx_3][j+1]; c_iny=liny[idx_1][j+3];
 982       c_bx=lbx[idx_2][j+1];c_by=lby[idx_1][j+2];
 983       for (s=0; s<n_seq; s++) {
 984         type2 = pair[S2[s][j+1]][S1[s][i-1]];
 985         in   +=P->mismatchI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s+di1+dj1;
 986         in_x +=iext_ass+di1;
 987         in_y +=iext_ass+dj1;
 988         in_xy+=iext_s+di1+dj1;
 989         inx  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s+di1+dj1;
 990         inx_x+=iext_ass+di1;
 991         iny  +=P->mismatch1nI[type2][S2[s][j]][S1[s][i]]+iopen+iext_s+di1+dj1;
 992         iny_y+=iext_ass+dj1;
 993         type2=pair[S2[s][j]][S1[s][i-1]];
 994         bx   +=bopen+bext+(type2>2?P->TerminalAU:0)+di1;
 995         bx_x +=bext+di1;
 996         type2=pair[S2[s][j+1]][S1[s][i]];
 997         by   +=bopen+bext+(type2>2?P->TerminalAU:0)+dj1;
 998         by_y +=bext+dj1;
 999       }
1000       lin[idx][j]=MIN2(in, MIN2(in_x, MIN2(in_y, in_xy)));
1001       linx[idx][j]=MIN2(inx_x, inx);
1002       liny[idx][j]=MIN2(iny_y, iny);
1003       lby[idx][j]=MIN2(by, by_y);
1004       lbx[idx][j]=MIN2(bx, bx_x);
1005       if (psc<MINPSCORE) continue;
1006       for (s=0; s<n_seq; s++) {
1007         lc[idx][j]+=E_ExtLoop(type[s], S1[s][i-1],S2[s][j+1],P);
1008       }
1009       for (s=0; s<n_seq; s++) {
1010         type2=pair[S1[s][i-1]][S2[s][j+1]];if (type2==0) type2=7;
1011         c_stack+=E_IntLoop(0,0,type2, rtype[type[s]],S1[s][i], S2[s][j], S1[s][i-1], S2[s][j+1], P)+di1+dj1;
1012         type2=pair[S1[s][i-1]][S2[s][j+2]];if (type2==0) type2=7;
1013         c_01   +=E_IntLoop(0,1,type2, rtype[type[s]],S1[s][i], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+di1+dj2;
1014         type2=pair[S1[s][i-2]][S2[s][j+1]];if (type2==0) type2=7;
1015         c_10   +=E_IntLoop(1,0,type2, rtype[type[s]],S1[s][i-1], S2[s][j], S1[s][i-1], S2[s][j+1], P)+di2+dj1;
1016         type2=pair[S1[s][i-2]][S2[s][j+2]];if (type2==0) type2=7;
1017         c_11   +=E_IntLoop(1,1,type2, rtype[type[s]],S1[s][i-1], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+di2+dj2;
1018         type2=pair[S1[s][i-3]][S2[s][j+3]];if (type2==0) type2=7;
1019         c_22   +=E_IntLoop(2,2,type2, rtype[type[s]],S1[s][i-2], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+di3+dj3;
1020         type2= pair[S1[s][i-3]][S2[s][j+2]];if (type2==0) type2=7;
1021         c_21   +=E_IntLoop(2,1,type2, rtype[type[s]],S1[s][i-2], S2[s][j+1], S1[s][i-1], S2[s][j+1], P)+di3+dj2;
1022         type2= pair[S1[s][i-2]][S2[s][j+3]];if (type2==0) type2=7;
1023         c_12   +=E_IntLoop(1,2,type2, rtype[type[s]],S1[s][i-1], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+di2+dj3;
1024         type2 = pair[S1[s][i-4]][S2[s][j+3]];if (type2==0) type2=7;
1025         c_32   +=E_IntLoop(3,2,type2, rtype[type[s]],S1[s][i-3], S2[s][j+2], S1[s][i-1], S2[s][j+1], P)+di4+dj3;
1026         type2 = pair[S1[s][i-3]][S2[s][j+4]];if (type2==0) type2=7;
1027         c_23   +=E_IntLoop(2,3,type2, rtype[type[s]],S1[s][i-2], S2[s][j+3], S1[s][i-1], S2[s][j+1], P)+dj4+di3;
1028         c_in   +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+di3+dj3+2*iext_s;
1029         c_in2x +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+di4+dj2;
1030         c_in2y +=P->mismatchI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_s+2*iext_ass+di2+dj4;
1031         c_inx  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+di3+dj1;
1032         c_iny  +=P->mismatch1nI[rtype[type[s]]][S1[s][i-1]][S2[s][j+1]]+iext_ass+iext_ass+di1+dj3;
1033         int bAU;
1034         bAU=(type[s]>2?P->TerminalAU:0);
1035         c_bx   +=di2+dj1+bext+bAU;
1036         c_by   +=di1+dj2+bext+bAU;
1037       }
1038       lc[idx][j] =MIN2(lc[idx][j],
1039                        MIN2(c_stack,
1040                             MIN2(c_10,
1041                                  MIN2(c_01,
1042                                       MIN2(c_11,
1043                                            MIN2(c_21,
1044                                                 MIN2(c_12,
1045                                                      MIN2(c_22,
1046                                                           MIN2(c_23,
1047                                                                MIN2(c_32,
1048                                                                     MIN2(c_bx,
1049                                                                          MIN2(c_by,
1050                                                                               MIN2(c_in,
1051                                                                                    MIN2(c_in2x,
1052                                                                                         MIN2(c_in2y,
1053                                                                                              MIN2(c_inx,c_iny)
1054                                                                                              )
1055                                                                                         )
1056                                                                                    )
1057                                                                               )
1058                                                                          )
1059                                                                     )
1060                                                                )
1061                                                           )
1062                                                      )
1063                                                 )
1064                                            )
1065                                       )
1066                                  )
1067                             )
1068                        );
1069       lc[idx][j]-=psc;
1070       temp=lc[idx][j];
1071
1072       for (s=0; s<n_seq; s++) {
1073         temp+=E_ExtLoop(rtype[type[s]], S2[s][j-1],S1[s][i+1],P);
1074       }
1075       if(min_colonne > temp){
1076         min_colonne=temp;
1077         min_j_colonne=j;
1078       }
1079     }
1080     if(max>=min_colonne){
1081             max=min_colonne;
1082             max_pos=i;
1083         max_pos_j=min_j_colonne;
1084     }
1085     position[i+delta]=min_colonne;min_colonne=INF;
1086     position_j[i+delta]=min_j_colonne;
1087     i++;
1088   }
1089   /* printf("MAX: %d",max); */
1090   for (s=0; s<n_seq; s++) {free(S1[s]);free(S2[s]);}
1091   free(S1); free(S2);
1092   if(max<threshold){
1093     alifind_max_XS(position, position_j, delta, threshold, alignment_length, s1, s2,access_s1,access_s2, fast);
1094   }
1095   aliplot_max_XS(max, max_pos, max_pos_j,alignment_length, s1, s2, access_s1, access_s2, fast);
1096   for (i=0; i<=4; i++) {free(lc[i]);free(lin[i]);free(lbx[i]);free(lby[i]);free(linx[i]);free(liny[i]);}
1097   /* free(lc[0]);free(lin[0]);free(lbx[0]);free(lby[0]);free(linx[0]);free(liny[0]); */
1098   free(lc);free(lin);free(lbx);free(lby);free(linx);free(liny);
1099   free(position);
1100   free(position_j);
1101   free(type);
1102   return NULL;
1103 }
1104
1105
1106
1107
1108 PRIVATE void alifind_max_XS(const int *position, const int *position_j,
1109                                 const int delta, const int threshold, const int alignment_length,
1110                                 const char *s1[], const char *s2[],
1111                                 const int **access_s1, const int **access_s2, const int fast){
1112
1113   int n_seq=0;
1114   for (n_seq=0; s1[n_seq]!=NULL; n_seq++);
1115   int pos=n1-9;
1116   if(fast==1){
1117     while(10 < pos--){
1118       int temp_min=0;
1119       if(position[pos+delta]<(threshold)){
1120         int search_range;
1121         search_range=delta+1;
1122         while(--search_range){
1123           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
1124             temp_min=search_range;
1125           }
1126         }
1127         pos-=temp_min;
1128         int max_pos_j;
1129         max_pos_j=position_j[pos+delta];
1130         int max;
1131         max=position[pos+delta];
1132         /*         printf("target upper bound %d: query lower bound %d  (%5.2f) \n", pos-10, max_pos_j-10, ((double)max)/100); */
1133         pos=MAX2(10,pos-delta);
1134       }
1135     }
1136   }
1137   else{
1138     pos=n1-9;
1139     while( pos-- > 10 ) {
1140       /* printf("delta %d position:%d value:%d\n", delta, pos, position[pos]); */
1141       int temp_min=0;
1142       if(position[pos+delta]<(threshold)){
1143         int search_range;
1144         search_range=delta+1;
1145         while(--search_range){
1146            if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
1147             temp_min=search_range;
1148           }
1149         }
1150         pos-=temp_min;
1151         int max_pos_j;
1152         max_pos_j=position_j[pos+delta]; /* position on j */
1153         int begin_t=MAX2(11, pos-alignment_length);
1154         int end_t  =MIN2(n1-10, pos+1);
1155         int begin_q=MAX2(11, max_pos_j-1);
1156         int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
1157         int i_flag;
1158         int j_flag;
1159         i_flag = (end_t   ==  pos+1?1:0);
1160         j_flag = (begin_q == max_pos_j-1?1:0);
1161         char **s3, **s4;
1162         s3 = (char**) space(sizeof(char*)*(n_seq+1));
1163         s4 = (char**) space(sizeof(char*)*(n_seq+1));
1164         int i;
1165         for(i=0; i<n_seq; i++){
1166           s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
1167           s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
1168           strncpy(s3[i], (s1[i]+begin_t), end_t - begin_t +1);
1169           strncpy(s4[i], (s2[i]+begin_q), end_q - begin_q +1);
1170           s3[i][end_t - begin_t +1]='\0';
1171           s4[i][end_q - begin_q +1]='\0';
1172         }
1173         duplexT test;
1174         test = aliduplexfold_XS((const char**) s3, (const char**) s4, access_s1, access_s2,pos, max_pos_j,threshold,i_flag,j_flag);
1175         /*         printf("position %d approximation %d test %d threshold %d\n", pos, position[pos+delta], (int)test.energy,(int)(threshold/n_seq)); */
1176         if(test.energy*100  < (int) (threshold/n_seq)){
1177         printf("%s %3d,%-3d: %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f)\n", test.structure,
1178                test.tb,test.te,test.qb,test.qe, test.ddG/n_seq, test.energy/n_seq, test.dG1/n_seq, test.dG2/n_seq);
1179         free(test.structure);
1180         pos=MAX2(10,pos-delta);
1181         }
1182         for(i=0;i<n_seq;i++){
1183           free(s3[i]);free(s4[i]);
1184         }
1185         free(s3);free(s4);
1186         /* free(test.structure); */
1187       }
1188     }
1189   }
1190 }
1191
1192 PRIVATE void aliplot_max_XS(const int max, const int max_pos, const int max_pos_j,const int alignment_length,
1193                             const char *s1[], const char* s2[],const int **access_s1, const int **access_s2,
1194                             const int fast){
1195
1196   int n_seq;
1197   for (n_seq=0; !(s1[n_seq]==NULL); n_seq++);
1198   n1 = strlen(s1[0]); /* get length of alignment */
1199   n2 = strlen(s2[0]); /* get length of alignme */
1200
1201   if(fast){
1202     printf("target upper bound %d: query lower bound %d (%5.2f)\n",
1203            max_pos-10, max_pos_j-10, (double) ((double)max)/(100*n_seq));
1204   }
1205   else{
1206     int begin_t=MAX2(11, max_pos-alignment_length); /* only get the position that binds.. */
1207     int end_t  =MIN2(n1-10, max_pos+1);              /* ..no dangles */
1208     int begin_q=MAX2(11, max_pos_j-1);
1209     int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
1210     int i_flag;
1211     int j_flag;
1212     i_flag = (end_t   == max_pos+1?1:0);
1213     j_flag = (begin_q == max_pos_j-1?1:0);
1214     char **s3, **s4;
1215     s3 = (char**) space(sizeof(char*)*(n_seq+1));
1216     s4 = (char**) space(sizeof(char*)*(n_seq+1));
1217     int i;
1218     for(i=0; i<n_seq; i++){
1219       s3[i] = (char*) space(sizeof(char)*(end_t-begin_t+2));
1220       s4[i] = (char*) space(sizeof(char)*(end_q-begin_q+2));
1221       strncpy(s3[i], (s1[i]+begin_t), end_t - begin_t +1);
1222       strncpy(s4[i], (s2[i]+begin_q), end_q - begin_q +1);
1223       s3[i][end_t - begin_t +1]='\0';
1224       s4[i][end_q - begin_q +1]='\0';
1225     }
1226     duplexT test;
1227     test = aliduplexfold_XS((const char**) s3, (const char**) s4, access_s1, access_s2,max_pos, max_pos_j,INF,i_flag,j_flag);
1228     printf("%s %3d,%-3d: %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f)\n", test.structure,
1229            test.tb,test.te,test.qb,test.qe, test.ddG/n_seq, test.energy/n_seq, test.dG1/n_seq, test.dG2/n_seq);
1230     free(test.structure);
1231     for(i=0;i<n_seq;i++){
1232       free(s3[i]);free(s4[i]);
1233     }
1234     free(s3);free(s4);
1235   }
1236 }
1237
1238 PRIVATE int covscore(const int *types, int n_seq) {
1239   /* calculate co-variance bonus for a pair depending on  */
1240   /* compensatory/consistent mutations and incompatible seqs */
1241   /* should be 0 for conserved pairs, >0 for good pairs      */
1242 #define NONE -10000 /* score for forbidden pairs */
1243   int k,l,s,score, pscore;
1244   int dm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
1245                 {0,0,2,2,1,2,2} /* CG */,
1246                 {0,2,0,1,2,2,2} /* GC */,
1247                 {0,2,1,0,2,1,2} /* GU */,
1248                 {0,1,2,2,0,2,1} /* UG */,
1249                 {0,2,2,1,2,0,2} /* AU */,
1250                 {0,2,2,2,1,2,0} /* UA */};
1251
1252   int pfreq[8]={0,0,0,0,0,0,0,0};
1253   for (s=0; s<n_seq; s++)
1254     pfreq[types[s]]++;
1255
1256   if (pfreq[0]*2>n_seq) return NONE;
1257   for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
1258     for (l=k+1; l<=6; l++)
1259       /* scores for replacements between pairtypes    */
1260       /* consistent or compensatory mutations score 1 or 2  */
1261       score += pfreq[k]*pfreq[l]*dm[k][l];
1262
1263   /* counter examples score -1, gap-gap scores -0.25   */
1264   pscore = cv_fact *
1265     ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
1266   return pscore;
1267 }
1268
1269
1270 PRIVATE void update_dfold_params(void)
1271 {
1272   if(P) free(P);
1273   P = scale_parameters();
1274   make_pair_matrix();
1275 }
1276
1277
1278
1279 PRIVATE short * encode_seq(const char *sequence) {
1280   unsigned int i,l;
1281   short *S;
1282   l = strlen(sequence);
1283   S = (short *) space(sizeof(short)*(l+2));
1284   S[0] = (short) l;
1285
1286   /* make numerical encoding of sequence */
1287   for (i=1; i<=l; i++)
1288     S[i]= (short) encode_char(toupper(sequence[i-1]));
1289
1290   /* for circular folding add first base at position n+1 */
1291   S[l+1] = S[1];
1292
1293   return S;
1294 }