7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
10 #include "dp_lib_header.h"
13 /*******************************************************************************/
14 /* myers and Miller */
16 /* makes DP between the the ns[0] sequences and the ns[1] */
18 /* for MODE, see the function get_dp_cost */
19 /*******************************************************************************/
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 */
27 last = sapp[-1] -= (k); \
29 last = *sapp++ = -(k); \
31 /* Append "Insert k" op */
34 { sapp[-1] = (k); *sapp++ = last; } \
36 last = *sapp++ = (k); \
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 )
43 int a,b, i, j, l,l1, l2, len;
49 /********Prepare Penalties******/
51 /********************************/
54 pos=aln2pos_simple ( A,-1, ns, l_s);
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));
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);
68 i=0;j=0;sapp=S; len=0;
69 while (!(i==l1 && j==l2))
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;}
79 A=realloc_aln2 ( A,A->max_n_seq,len+1);
80 char_buf=declare_char (A->max_n_seq,len+1);
82 i=0;j=0;sapp=S; len=0;
83 while (!(i==l1 && j==l2))
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];
97 for ( a=0; a<l; a++, j++, len++)
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];
108 for ( a=0; a<l; a++, i++, len++)
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]='-';
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]]);}
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);
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)
143 /* Forward cost-only vectors */
146 /* Reverse cost-only vectors */
147 int midi, midj, type; /* Midpoint, type, and cost */
150 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
151 /*TG_MODE=0---> gop and gep*/
152 /*TG_MODE=1---> --- gep*/
153 /*TG_MODE=2---> --- ---*/
166 CC=vcalloc (L, sizeof (int));
167 DD=vcalloc (L, sizeof (int));
168 RR=vcalloc (L, sizeof (int));
169 SS=vcalloc (L, sizeof (int));
193 if (N <= 0){if (M > 0) DEL(M);return gap(M);}
203 if (tb > te) tb = te;
204 midc = (tb+h) + gap(N);
207 for (j = 1; j <= N; j++)
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);
220 {if (midj > 1) INS(midj-1);
222 if (midj < N) INS(N-midj);
227 /* Divide: Find optimum midpoint (midi,midj) of cost midc */
230 midi = M/2; /* Forward phase: */
231 CC[0] = 0; /* Compute C(M/2,k) & D(M/2,k) for all k */
233 for (j = 1; j <= N; j++)
238 for (i = 1; i <= midi; i++)
244 for (j = 1; j <= N; j++)
248 if ((c = c + m) > (e = e + h)) e = c;
249 if ((c = CC[j] + m) > (d = DD[j] + h)) d = c;
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);
264 RR[N] = 0; /* Reverse phase: */
268 for (j = N-1; j >= 0; j--)
273 for (i = M-1; i >= midi; i--)
277 for (j = N-1; j >= 0; j--)
279 if ((c = c + m) > (e = e + h)) e = c;
280 if ((c = RR[j] + m) > (d = SS[j] + h)) d = c;
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);
295 midc = CC[0]+RR[0]; /* Find optimal midpoint */
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]))
305 for (j = N; j >= 0; j--)
306 if ((c = DD[j] + SS[j] - g) > midc)
312 /* Conquer: recursively around midpoint */
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);
322 diff (A,ns, l_s, s1,midi-1, s2, midj, tb,0, CL, pos);
324 diff (A,ns, l_s, s1+midi+1, M-midi-1,s2+midj, N-midj,0,te, CL, pos);
343 /******************************COPYRIGHT NOTICE*******************************/
344 /*© Centro de Regulacio Genomica */
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.*/
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.*/
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.*/
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 /*............................................... |*/
371 /******************************COPYRIGHT NOTICE*******************************/