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 int myers_miller_pair_wise (Alignment *A,int *ns, int **l_s,Constraint_list *CL )
  41         {
  42         int **pos;
  43         int a,b, i, j, l,l1, l2, len;
  44         int *S;
  45         char ** char_buf;
  46         int score;
  47
  48
  49         /********Prepare Penalties******/
  50
  51         /********************************/
  52
  53
  54         pos=aln2pos_simple ( A,-1, ns, l_s);
  55
  56
  57         l1=strlen (A->seq_al[l_s[0][0]]);
  58         l2=strlen (A->seq_al[l_s[1][0]]);
  59         S=vcalloc (l1+l2+1, sizeof (int));
  60         last=0;
  61         sapp=S;
  62
  63         score=diff (A,ns, l_s, 0, l1, 0, l2, 0, 0, CL, pos);
  64         diff (NULL,ns, l_s, 0, l1, 0, l2, 0, 0, CL, pos);
  65
  66
  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,ma;
 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                  ma=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
 257                  s = CC[j];
 258                  CC[j] = c;
 259                  DD[j] = d;
 260                  }
 261              }
 262          DD[0] = CC[0];
 263
 264          RR[N] = 0;            /* Reverse phase:                          */
 265          t = te;
 266
 267
 268          for (j = N-1; j >= 0; j--)
 269              { RR[j] = t = t+h;
 270              SS[j] = t+g;
 271              }
 272          t = te;
 273          for (i = M-1; i >= midi; i--)
 274              { s = RR[N];
 275              RR[N] = c = t = t+h;
 276              e = t+g;
 277              for (j = N-1; j >= 0; j--)
 278                   {
 279                   if ((c =   c   + m) > (e =   e   + h)) e = c;
 280                   if ((c = RR[j] + m) > (d = SS[j] + h)) d = c;
 281
 282                   ma=c = s + (CL->get_dp_cost) (A, pos, ns[0], l_s[0],i+s1, pos, ns[1], l_s[1],j+s2,CL);
 283
 284                   if (e > c) c = e;
 285                   if (d > c) c = d;
 286
 287
 288                   s = RR[j];
 289                   RR[j] = c;
 290                   SS[j] = d;
 291
 292                   }
 293              }
 294          SS[N] = RR[N];
 295          midc = CC[0]+RR[0];        /* Find optimal midpoint */
 296          midj = 0;
 297          type = 1;
 298          for (j = 0; j <= N; j++)
 299              if ((c = CC[j] + RR[j]) >= midc)
 300                if (c > midc || (CC[j] != DD[j] && RR[j] == SS[j]))
 301                     {
 302                     midc = c;
 303                     midj = j;
 304                     }
 305          for (j = N; j >= 0; j--)
 306              if ((c = DD[j] + SS[j] - g) > midc)
 307                  {midc = c;
 308                  midj = j;
 309                  type = 2;
 310                  }
 311          }
 312 /* Conquer: recursively around midpoint */
 313
 314   if (type == 1)
 315     {
 316
 317     diff (A,ns, l_s, s1,midi, s2, midj, tb, CL->gop*SCORE_K, CL, pos);
 318     diff (A,ns, l_s, s1+midi,M-midi, s2+midj, N-midj, CL->gop*SCORE_K,te, CL, pos);
 319     }
 320   else
 321     {
 322       diff (A,ns, l_s, s1,midi-1, s2, midj, tb,0, CL, pos);
 323       DEL(2);
 324       diff (A,ns, l_s, s1+midi+1, M-midi-1,s2+midj, N-midj,0,te, CL, pos);
 325     }
 326   return midc;
 327   }
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343 /******************************COPYRIGHT NOTICE*******************************/
 344 /*© Centro de Regulacio Genomica */
 345 /*and */
 346 /*Cedric Notredame */
 347 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
 348 /*All rights reserved.*/
 349 /*This file is part of T-COFFEE.*/
 350 /**/
 351 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
 352 /*    it under the terms of the GNU General Public License as published by*/
 353 /*    the Free Software Foundation; either version 2 of the License, or*/
 354 /*    (at your option) any later version.*/
 355 /**/
 356 /*    T-COFFEE is distributed in the hope that it will be useful,*/
 357 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
 358 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
 359 /*    GNU General Public License for more details.*/
 360 /**/
 361 /*    You should have received a copy of the GNU General Public License*/
 362 /*    along with Foobar; if not, write to the Free Software*/
 363 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
 364 /*...............................................                                                                                      |*/
 365 /*  If you need some more information*/
 366 /*  cedric.notredame@europe.com*/
 367 /*...............................................                                                                                                                                     |*/
 368 /**/
 369 /**/
 370 /*      */
 371 /******************************COPYRIGHT NOTICE*******************************/