binaries/src/tcoffee/t_coffee_source/util_dp_mm_nw.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6
   7 #include "io_lib_header.h"
   8 #include "util_lib_header.h"
   9 #include "define_header.h"
  10 #include "dp_lib_header.h"
  11
  12
  13 /*******************************************************************************/
  14 /*                myers and Miller                                             */
  15 /*                                                                             */
  16 /*      makes DP between the the ns[0] sequences and the ns[1]                 */
  17 /*                                                                             */
  18 /*      for MODE, see the function  get_dp_cost                                */
  19 /*******************************************************************************/
  20
  21 #define gap(k)  ((k) <= 0 ? 0 : g+h*(k))    /* k-symbol indel cost */
  22 static int *sapp;                /* Current script append ptr */
  23 static int  last;                /* Last script op appended */
  24                        /* Append "Delete k" op */
  25 #define DEL(k)                \
  26 { if (last < 0)                \
  27     last = sapp[-1] -= (k);        \
  28   else                    \
  29     last = *sapp++ = -(k);        \
  30 }
  31                         /* Append "Insert k" op */
  32 #define INS(k)                \
  33 { if (last < 0)                \
  34     { sapp[-1] = (k); *sapp++ = last; }    \
  35   else                    \
  36     last = *sapp++ = (k);        \
  37 }
  38
  39 #define REP { last = *sapp++ = 0; }        /* Append "Replace" op */
  40
  41 int myers_miller_pair_wise (Alignment *A,int *ns, int **l_s,Constraint_list *CL )
  42         {
  43         int **pos;
  44         int a,b, i, j, l,l1, l2, len;
  45         int *S;
  46         char ** char_buf;
  47         int score;
  48
  49
  50         /********Prepare Penalties******/
  51
  52         /********************************/
  53
  54
  55         pos=aln2pos_simple ( A,-1, ns, l_s);
  56
  57
  58         l1=strlen (A->seq_al[l_s[0][0]]);
  59         l2=strlen (A->seq_al[l_s[1][0]]);
  60         S=vcalloc (l1+l2+1, sizeof (int));
  61         last=0;
  62         sapp=S;
  63
  64
  65         score=diff (A,ns, l_s, 0, l1, 0, l2, 0, 0, CL, pos);
  66         diff (NULL,ns, l_s, 0, l1, 0, l2, 0, 0, CL, pos);
  67
  68         i=0;j=0;sapp=S; len=0;
  69         while (!(i==l1 && j==l2))
  70               {
  71                   if (*sapp==0){i++; j++;len++;}
  72                   else if ( *sapp<0){i-=*sapp;len-=*sapp;}
  73                   else if ( *sapp>0){j+=*sapp;len+=*sapp;}
  74                   sapp++;
  75               }
  76
  77
  78
  79         A=realloc_aln2  ( A,A->max_n_seq,len+1);
  80         char_buf=declare_char (A->max_n_seq,len+1);
  81
  82         i=0;j=0;sapp=S; len=0;
  83         while (!(i==l1 && j==l2))
  84               {
  85
  86                    if (*sapp==0)
  87                       {
  88                           for (b=0; b< ns[0]; b++)
  89                               char_buf[l_s[0][b]][len]=A->seq_al[l_s[0][b]][i];
  90                           for (b=0; b< ns[1]; b++)
  91                               char_buf[l_s[1][b]][len]=A->seq_al[l_s[1][b]][j];
  92                           i++; j++;len++;
  93                       }
  94                   else if ( *sapp>0)
  95                         {
  96                             l=*sapp;
  97                             for ( a=0; a<l; a++, j++, len++)
  98                                 {
  99                                 for (b=0; b< ns[0]; b++)
 100                                     char_buf[l_s[0][b]][len]='-';
 101                                 for (b=0; b< ns[1]; b++)
 102                                     char_buf[l_s[1][b]][len]=A->seq_al[l_s[1][b]][j];
 103                                 }
 104                         }
 105                   else if ( *sapp<0)
 106                         {
 107                             l=-*sapp;
 108                             for ( a=0; a<l; a++, i++, len++)
 109                                 {
 110                                 for (b=0; b< ns[0]; b++)
 111                                     char_buf[l_s[0][b]][len]=A->seq_al[l_s[0][b]][i];;
 112                                 for (b=0; b< ns[1]; b++)
 113                                     char_buf[l_s[1][b]][len]='-';
 114                                 }
 115                         }
 116
 117                   sapp++;
 118               }
 119
 120
 121         A->len_aln=len;
 122         A->nseq=ns[0]+ns[1];
 123
 124         for ( a=0; a< ns[0]; a++){char_buf[l_s[0][a]][len]='\0'; sprintf ( A->seq_al[l_s[0][a]], "%s", char_buf[l_s[0][a]]);}
 125         for ( a=0; a< ns[1]; a++){char_buf[l_s[1][a]][len]='\0'; sprintf ( A->seq_al[l_s[1][a]], "%s", char_buf[l_s[1][a]]);}
 126
 127
 128         vfree (S);
 129         free_char ( char_buf, -1);
 130         l1=strlen (A->seq_al[l_s[0][0]]);
 131         l2=strlen (A->seq_al[l_s[1][0]]);
 132         if ( l1!=l2) exit(1);
 133
 134         free_int (pos, -1);
 135         return score;
 136         }
 137
 138
 139 int diff (Alignment *A, int *ns, int **l_s, int s1, int M,int s2, int N , int tb, int te, Constraint_list *CL, int **pos)
 140         {
 141          static int *CC;
 142          static int *DD;
 143              /* Forward cost-only vectors */
 144          static int *RR;
 145          static int *SS;
 146              /* Reverse cost-only vectors */
 147          int   midi, midj, type;    /* Midpoint, type, and cost */
 148          int  midc;
 149
 150 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 151 /*TG_MODE=0---> gop and gep*/
 152 /*TG_MODE=1---> ---     gep*/
 153 /*TG_MODE=2---> ---     ---*/
 154
 155
 156
 157
 158
 159
 160
 161          if ( !CC)
 162             {
 163               int L;
 164               L=M+N+1;
 165
 166               CC=vcalloc (L, sizeof (int));
 167               DD=vcalloc (L, sizeof (int));
 168               RR=vcalloc (L, sizeof (int));
 169               SS=vcalloc (L, sizeof (int));
 170             }
 171
 172          if ( A==NULL)
 173             {
 174               vfree(CC);
 175               vfree(DD);
 176               vfree(RR);
 177               vfree(SS);
 178               CC=DD=RR=SS=NULL;
 179               return 0;
 180             }
 181
 182          {
 183          int   i, j;
 184          int   c, e, d, s;
 185          int t, g,h,m;
 186
 187
 188
 189          g=CL->gop*SCORE_K;
 190          h=CL->gep*SCORE_K;
 191          m=g+h;
 192
 193          if (N <= 0){if (M > 0) DEL(M);return gap(M);}
 194          if (M <= 1)
 195                {
 196
 197                if (M <= 0)
 198                   {INS(N);
 199                   return gap(N);
 200                   }
 201
 202
 203                if (tb > te) tb = te;
 204                midc = (tb+h) + gap(N);
 205                midj = 0;
 206
 207                for (j = 1; j <= N; j++)
 208                    {
 209
 210                      c = gap(j-1) +(CL->get_dp_cost) (A, pos, ns[0], l_s[0],s1, pos, ns[1], l_s[1],j-1+s2,CL)+ gap(N-j);
 211
 212                      if (c > midc)
 213                        { midc = c;
 214                          midj = j;
 215                        }
 216                    }
 217                if (midj == 0)
 218                   {DEL(1) INS(N)}
 219                else
 220                   {if (midj > 1) INS(midj-1);
 221                   REP;
 222                   if (midj < N) INS(N-midj);
 223                   }
 224
 225                return midc;
 226                }
 227 /* Divide: Find optimum midpoint (midi,midj) of cost midc */
 228
 229
 230          midi = M/2;            /* Forward phase:                          */
 231          CC[0] = 0;            /*   Compute C(M/2,k) & D(M/2,k) for all k */
 232          t = tb;
 233          for (j = 1; j <= N; j++)
 234              { CC[j] = t = t+h;
 235                DD[j] = t+g;
 236              }
 237          t = tb;
 238          for (i = 1; i <= midi; i++)
 239              {
 240              s = CC[0];
 241              CC[0] = c = t = t+h;
 242              e = t+g;
 243
 244              for (j = 1; j <= N; j++)
 245                  {
 246
 247
 248                  if ((c =   c   + m) > (e =   e   + h)) e = c;
 249                  if ((c = CC[j] + m) > (d = DD[j] + h)) d = c;
 250
 251                  c = s + (CL->get_dp_cost) (A, pos, ns[0], l_s[0],i-1+s1, pos, ns[1], l_s[1],j-1+s2,CL);
 252
 253                  if (e > c) c = e;
 254                  if (d > c) c = d;
 255
 256                  s = CC[j];
 257                  CC[j] = c;
 258                  DD[j] = d;
 259                  }
 260              }
 261          DD[0] = CC[0];
 262
 263          RR[N] = 0;            /* Reverse phase:                          */
 264          t = te;
 265
 266
 267          for (j = N-1; j >= 0; j--)
 268              { RR[j] = t = t+h;
 269              SS[j] = t+g;
 270              }
 271          t = te;
 272          for (i = M-1; i >= midi; i--)
 273              { s = RR[N];
 274              RR[N] = c = t = t+h;
 275              e = t+g;
 276              for (j = N-1; j >= 0; j--)
 277                   {
 278                   if ((c =   c   + m) > (e =   e   + h)) e = c;
 279                   if ((c = RR[j] + m) > (d = SS[j] + h)) d = c;
 280
 281                   c = s + (CL->get_dp_cost) (A, pos, ns[0], l_s[0],i+s1, pos, ns[1], l_s[1],j+s2,CL);
 282
 283                   if (e > c) c = e;
 284                   if (d > c) c = d;
 285                   s = RR[j];
 286                   RR[j] = c;
 287                   SS[j] = d;
 288
 289                   }
 290              }
 291          SS[N] = RR[N];
 292          midc = CC[0]+RR[0];        /* Find optimal midpoint */
 293          midj = 0;
 294          type = 1;
 295          for (j = 0; j <= N; j++)
 296              if ((c = CC[j] + RR[j]) >= midc)
 297                if (c > midc || (CC[j] != DD[j] && RR[j] == SS[j]))
 298                     {
 299                     midc = c;
 300                     midj = j;
 301                     }
 302          for (j = N; j >= 0; j--)
 303              if ((c = DD[j] + SS[j] - g) > midc)
 304                  {midc = c;
 305                  midj = j;
 306                  type = 2;
 307                  }
 308          }
 309 /* Conquer: recursively around midpoint */
 310
 311   if (type == 1)
 312     {
 313
 314     diff (A,ns, l_s, s1,midi, s2, midj, tb, CL->gop*SCORE_K, CL, pos);
 315     diff (A,ns, l_s, s1+midi,M-midi, s2+midj, N-midj, CL->gop*SCORE_K,te, CL, pos);
 316     }
 317   else
 318     {
 319       diff (A,ns, l_s, s1,midi-1, s2, midj, tb,0, CL, pos);
 320       DEL(2);
 321       diff (A,ns, l_s, s1+midi+1, M-midi-1,s2+midj, N-midj,0,te, CL, pos);
 322     }
 323   return midc;
 324   }
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340 /*********************************COPYRIGHT NOTICE**********************************/
 341 /*© Centre National de la Recherche Scientifique (CNRS) */
 342 /*and */
 343 /*Cedric Notredame */
 344 /*Wed Sep 21 19:11:38     2005. */
 345 /*All rights reserved.*/
 346 /*This file is part of T-COFFEE.*/
 347 /**/
 348 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
 349 /*    it under the terms of the GNU General Public License as published by*/
 350 /*    the Free Software Foundation; either version 2 of the License, or*/
 351 /*    (at your option) any later version.*/
 352 /**/
 353 /*    T-COFFEE is distributed in the hope that it will be useful,*/
 354 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
 355 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
 356 /*    GNU General Public License for more details.*/
 357 /**/
 358 /*    You should have received a copy of the GNU General Public License*/
 359 /*    along with Foobar; if not, write to the Free Software*/
 360 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
 361 /*...............................................                                                                                      |*/
 362 /*  If you need some more information*/
 363 /*  cedric.notredame@europe.com*/
 364 /*...............................................                                                                                                                                     |*/
 365 /**/
 366 /**/
 367 /*      */
 368 /*********************************COPYRIGHT NOTICE**********************************/
 369 /*********************************COPYRIGHT NOTICE**********************************/
 370 /*© Centro de Regulacio Genomica */
 371 /*and */
 372 /*Cedric Notredame */
 373 /*Tue Oct 27 10:12:26 WEST 2009. */
 374 /*All rights reserved.*/
 375 /*This file is part of T-COFFEE.*/
 376 /**/
 377 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
 378 /*    it under the terms of the GNU General Public License as published by*/
 379 /*    the Free Software Foundation; either version 2 of the License, or*/
 380 /*    (at your option) any later version.*/
 381 /**/
 382 /*    T-COFFEE is distributed in the hope that it will be useful,*/
 383 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
 384 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
 385 /*    GNU General Public License for more details.*/
 386 /**/
 387 /*    You should have received a copy of the GNU General Public License*/
 388 /*    along with Foobar; if not, write to the Free Software*/
 389 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
 390 /*...............................................                                                                                      |*/
 391 /*  If you need some more information*/
 392 /*  cedric.notredame@europe.com*/
 393 /*...............................................                                                                                                                                     |*/
 394 /**/
 395 /**/
 396 /*      */
 397 /*********************************COPYRIGHT NOTICE**********************************/