7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
10 #include "dp_lib_header.h"
13 int cfasta_cdna_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
16 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
17 /*TG_MODE=0---> gop and gep*/
18 /*TG_MODE=1---> --- gep*/
19 /*TG_MODE=2---> --- ---*/
22 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
28 static char **group_list;
32 int max_n_chosen_diag;
35 /********Prepare Penalties******/
37 maximise=CL->maximise;
40 /********************************/
44 group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
46 B=dna_aln2_3frame_cdna_aln(A, ns, l_s);
47 l0=strlen(B->seq_al[0]);
48 l1=strlen(B->seq_al[3]);
49 tot_diag=evaluate_diagonals_cdna ( B, ns, l_s, CL, maximise,n_groups,group_list, ktup);
51 max_n_chosen_diag=100;
52 n_chosen_diag=step=10 ;
53 n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);
56 diag=extract_N_diag (l0,l1, tot_diag, n_chosen_diag,2);
57 score =make_fasta_cdna_pair_wise ( A,B, ns, l_s, CL, diag);
63 while (new_score!=score && n_chosen_diag< max_n_chosen_diag )
67 ungap_sub_aln ( A, ns[0], l_s[0]);
68 ungap_sub_aln ( A, ns[1], l_s[1]);
72 n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);
74 diag =extract_N_diag (l0,l1, tot_diag, n_chosen_diag,3);
75 new_score=make_fasta_cdna_pair_wise ( A, B,ns, l_s, CL, diag);
80 free_int (tot_diag, -1);
86 int fasta_cdna_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
88 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
89 /*TG_MODE=0---> gop and gep*/
90 /*TG_MODE=1---> --- gep*/
91 /*TG_MODE=2---> --- ---*/
96 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
101 static char **group_list;
105 /********Prepare Penalties******/
108 maximise=CL->maximise;
110 /********************************/
117 group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
120 B=dna_aln2_3frame_cdna_aln(A, ns, l_s);
123 l0=strlen ( B->seq_al[0]);
124 l1=strlen ( B->seq_al[3]);
127 tot_diag=evaluate_diagonals_cdna( B, ns, l_s, CL, maximise,n_groups,group_list, ktup);
130 diag=extract_N_diag (l0, l1, tot_diag,20,1);
131 score=make_fasta_cdna_pair_wise ( A,B, ns, l_s, CL, diag);
134 free_int (tot_diag, -1);
139 Dp_Model* initialize_dna_dp_model (Constraint_list *CL)
144 int deltaf1, deltaf0,deltatype;
145 int type, type1, type0;
147 M=vcalloc ( 1, sizeof (Dp_Model));
149 for (M->nstate=0,f0=0; f0<3; f0++)
150 for ( f1=0; f1<3; f1++)M->nstate+=3;
153 M->START=M->nstate++;
156 M->TG_MODE=CL->TG_MODE;
158 M->gop=CL->gop*SCORE_K;
159 M->gep=CL->gep*SCORE_K;
163 M->f_gop=CL->f_gop*SCORE_K;
164 M->f_gep=CL->f_gep*SCORE_K;
167 M->bounded_model=declare_int (M->nstate+1, M->nstate+1);
168 M->model=declare_int (M->nstate+1, M->nstate+1);
169 for ( a=0; a<=M->nstate; a++)
170 for ( b=0; b<= M->nstate; b++)
171 M->model[a][b]=UNDEFINED;
172 M->model_properties=declare_int ( M->nstate, 10);
175 M->TYPE=a++;M->F0=a++;M->F1=a++; M->LEN_I=a++; M->LEN_J=a++; M->DELTA_I=a++;M->DELTA_J=a++;M->EMISSION=a++;M->TERM_EMISSION=a++;
177 M->NON_CODING=a++; M->INSERTION=a++; M->DELETION=a++; M->CODING0=a++; M->CODING1=a++;M->CODING2=a++;
180 for ( a=0,f0=0; f0<3; f0++)
181 for ( f1=0; f1<3; f1++, a+=3)
183 M->model_properties[a+0][M->TYPE]=M->CODING0;
184 M->model_properties[a+0][M->F0]=f0;
185 M->model_properties[a+0][M->F1]=f1;
186 M->model_properties[a+0][M->LEN_I]=1;
187 M->model_properties[a+0][M->LEN_J]=1;
188 M->model_properties[a+0][M->DELTA_I]=-1;
189 M->model_properties[a+0][M->DELTA_J]= 0;
190 M->model_properties[a+0][M->EMISSION]=0;
191 M->model_properties[a+0][M->TERM_EMISSION]=0;
193 M->model_properties[a+1][M->TYPE]=M->DELETION;
194 M->model_properties[a+1][M->F0]=f0;
195 M->model_properties[a+1][M->F1]=f1;
196 M->model_properties[a+1][M->LEN_I]=1;
197 M->model_properties[a+1][M->LEN_J]=0;
198 M->model_properties[a+1][M->DELTA_I]=-1;
199 M->model_properties[a+1][M->DELTA_J]=+1;
200 M->model_properties[a+1][M->EMISSION]=M->gep;
201 M->model_properties[a+1][M->TERM_EMISSION]=(M->TG_MODE==2)?0:M->gep;
203 M->model_properties[a+2][M->TYPE]=M->INSERTION;
204 M->model_properties[a+2][M->F0]=f0;
205 M->model_properties[a+2][M->F1]=f1;
206 M->model_properties[a+2][M->LEN_I]=0;
207 M->model_properties[a+2][M->LEN_J]=1;
208 M->model_properties[a+2][M->DELTA_I]= 0;
209 M->model_properties[a+2][M->DELTA_J]=-1;
210 M->model_properties[a+2][M->EMISSION]=M->gep;
211 M->model_properties[a+2][M->TERM_EMISSION]=(M->TG_MODE==2)?0:M->gep;
214 /*UM" Unmatched State*/
215 M->model_properties[a][M->TYPE]=M->NON_CODING;
216 M->model_properties[a][M->F0]=0;
217 M->model_properties[a][M->F1]=0;
218 M->model_properties[a][M->LEN_I]=1;
219 M->model_properties[a][M->LEN_J]=1;
220 M->model_properties[a][M->DELTA_I]=-1;
221 M->model_properties[a][M->DELTA_J]=0;
222 M->model_properties[a][M->EMISSION]=M->f_gep;
223 M->model_properties[a][M->TERM_EMISSION]=(M->F_TG_MODE==2)?0:M->f_gep;
226 M->model_properties[M->START][M->TYPE]=M->NON_CODING;
227 M->model_properties[a+1][M->F0]=0;
228 M->model_properties[a+1][M->F1]=0;
229 M->model_properties[a+1][M->LEN_I]=0;
230 M->model_properties[a+1][M->LEN_J]=0;
231 M->model_properties[a+1][M->DELTA_I]=0 ;
232 M->model_properties[a+1][M->DELTA_J]=0;
233 M->model_properties[a+1][M->EMISSION]=0;
234 M->model_properties[a+1][M->TERM_EMISSION]=0;
236 M->model_properties[M->END][M->TYPE]=M->NON_CODING;
237 M->model_properties[a+2][M->F0]=0;
238 M->model_properties[a+2][M->F1]=0;
239 M->model_properties[a+2][M->LEN_I]=0;
240 M->model_properties[a+2][M->LEN_J]=0;
241 M->model_properties[a+2][M->DELTA_I]=0 ;
242 M->model_properties[a+2][M->DELTA_J]=0;
243 M->model_properties[a+2][M->EMISSION]=0;
244 M->model_properties[a+2][M->TERM_EMISSION]=0;
246 /*1: SET THE INDEL PENALTIES*/
249 for ( a=0; a< M->START; a++)
251 deltaf0=M->model_properties[M->START][M->F0]-M->model_properties[a][M->F0];
252 deltaf1=M->model_properties[M->START][M->F1]-M->model_properties[a][M->F1];
253 if ( deltaf0==0 && deltaf1==0)deltatype=0;
254 else if ( deltaf0<=0 && deltaf1<=0)deltatype=1;
256 type=M->model_properties[a][M->TYPE];
258 if ( type==M->NON_CODING) M->model[a][M->END]=M->model[M->START][a]=(M->F_TG_MODE==0)?M->f_gop:0;
259 else if ( type==M->CODING0 && deltatype==0)M->model[a][M->END]=M->model[M->START][a]=ALLOWED;
260 else if ( type==M->CODING0 && deltatype==1)M->model[a][M->END]=M->model[M->START][a]=(M->F_TG_MODE==0)?M->f_gop:0;
261 else if ( type==M->INSERTION && deltatype==0)M->model[a][M->END]=M->model[M->START][a]=(M->TG_MODE==0)?M->gop:0;
262 else if ( type==M->INSERTION && deltatype==1)M->model[a][M->END]=M->model[M->START][a]=(M->TG_MODE==0)?M->gop:0+(M->F_TG_MODE==0)?M->f_gop:0;
263 else if ( type==M->DELETION && deltatype==0)M->model[a][M->END]=M->model[M->START][a]=(M->TG_MODE==0)?M->gop:0;
264 else if ( type==M->DELETION && deltatype==1)M->model[a][M->END]=M->model[M->START][a]=(M->TG_MODE==0)?M->gop:0+(M->F_TG_MODE==0)?M->f_gop:0;
265 else M->model[a][M->END]=M->model[M->START][a]=UNDEFINED;
268 if (type==M->NON_CODING ||M->model_properties[a][M->F0] ||M->model_properties[a][M->F1]) M->model[M->START][a]=M->model[a][M->END]=(M->F_TG_MODE==0)?M->f_gop:0;
269 else if (type==M->INSERTION || type==M->DELETION)M->model[M->START][a]=M->model[a][M->END]=( M->TG_MODE==0)?M->gop:0;
270 else M->model[M->START][a]=M->model[a][M->END]=ALLOWED;
273 for ( b=0; b< M->START; b++)
276 deltaf0=M->model_properties[a][M->F0]-M->model_properties[b][M->F0];
277 deltaf1=M->model_properties[a][M->F1]-M->model_properties[b][M->F1];
278 type0=M->model_properties[a][M->TYPE];
279 type1=M->model_properties[b][M->TYPE];
281 if ( deltaf0==0 && deltaf1==0)deltatype=0;
282 else if ( deltaf0<=0 && deltaf1<=0)deltatype=1;
287 if ( type0==M->NON_CODING && type1==M->NON_CODING )M->model[a][b]=UNDEFINED;
288 else if ( type0==M->NON_CODING && type1==M->CODING0 )M->model[a][b]=ALLOWED ;
289 else if ( type0==M->NON_CODING && type1==M->INSERTION )M->model[a][b]=M->gop;
290 else if ( type0==M->NON_CODING && type1==M->DELETION )M->model[a][b]=M->gop;
292 else if ( type0==M->CODING0 && type1==M->NON_CODING )M->model[a][b]=M->f_gop;
293 else if ( type0==M->CODING0 && type1==M->CODING0 && deltatype==0 )M->model[a][b]=ALLOWED;
294 else if ( type0==M->CODING0 && type1==M->CODING0 && deltatype==1 )M->model[a][b]=M->f_gop;
295 else if ( type0==M->CODING0 && type1==M->INSERTION && deltatype==0 )M->model[a][b]=M->gop;
296 else if ( type0==M->CODING0 && type1==M->INSERTION && deltatype==1 )M->model[a][b]=M->gop+M->f_gop;
297 else if ( type0==M->CODING0 && type1==M->DELETION && deltatype==0 )M->model[a][b]=M->gop;
298 else if ( type0==M->CODING0 && type1==M->DELETION && deltatype==1 )M->model[a][b]=M->gop+M->f_gop;
300 else if ( type0==M->INSERTION && type1==M->NON_CODING )M->model[a][b]=M->f_gop;
301 else if ( type0==M->INSERTION && type1==M->CODING0 && deltatype==0 )M->model[a][b]=ALLOWED;
302 else if ( type0==M->INSERTION && type1==M->CODING0 && deltatype==1 )M->model[a][b]=M->f_gop;
303 else if ( type0==M->INSERTION && type1==M->INSERTION && deltatype==0 )M->model[a][b]=ALLOWED;
304 else if ( type0==M->INSERTION && type1==M->INSERTION && deltatype==1 )M->model[a][b]=M->f_gop;
305 else if ( type0==M->INSERTION && type1==M->DELETION && deltatype==0 )M->model[a][b]=M->gop;
306 else if ( type0==M->INSERTION && type1==M->DELETION && deltatype==1 )M->model[a][b]=M->gop+M->f_gop;
308 else if ( type0==M->DELETION && type1==M->NON_CODING )M->model[a][b]=M->f_gop;
309 else if ( type0==M->DELETION && type1==M->CODING0 && deltatype==0 )M->model[a][b]=ALLOWED;
310 else if ( type0==M->DELETION && type1==M->CODING0 && deltatype==1 )M->model[a][b]=M->f_gop;
311 else if ( type0==M->DELETION && type1==M->INSERTION && deltatype==0 )M->model[a][b]=M->gop;
312 else if ( type0==M->DELETION && type1==M->INSERTION && deltatype==1 )M->model[a][b]=M->gop+M->f_gop;
313 else if ( type0==M->DELETION && type1==M->DELETION && deltatype==0 )M->model[a][b]=ALLOWED;
314 else if ( type0==M->DELETION && type1==M->DELETION && deltatype==1 )M->model[a][b]=M->f_gop;
316 else {M->model[a][b]=UNDEFINED;}
322 /*2 SET THE FRAMESHIFT PENALTIES
324 for ( a=0; a< M->START; a++)
326 type=M->model_properties[a][M->TYPE];
328 for ( b=0; b< M->START; b++)
330 deltaf0=M->model_properties[a][M->F0]-M->model_properties[b][M->F0];
331 deltaf1=M->model_properties[a][M->F1]-M->model_properties[b][M->F1];
336 if (b==M->UM) M->model[a][b]+=M->f_gop;
337 else if (a==M->UM) M->model[a][b]+=ALLOWED;
338 else if (deltaf1==0 && deltaf0==0)M->model[a][b]+=ALLOWED;
339 else if (deltaf1<=0 && deltaf0<=0)M->model[a][b]+=M->f_gop;
340 else M->model[a][b]=UNDEFINED;
344 M->model[M->UM][M->UM]=UNDEFINED;
348 for (c=0,a=0, d=0; a< M->START; a++)
349 for ( b=0; b<M->START; b++, d++)
351 if (M->model[a][b]!=UNDEFINED)
353 M->bounded_model[b][1+M->bounded_model[b][0]++]=a;
359 int make_fasta_cdna_pair_wise (Alignment *B,Alignment *A,int*in_ns, int **l_s,Constraint_list *CL, int *diag)
367 int deltaf0, deltaf1, delta;
369 int ala, alb, aa0, aa1;
376 int debug_cdna_fasta=0;
379 int state,prev_state;
384 l0=strlen ( B->seq_al[l_s[0][0]]);
385 l1=strlen ( B->seq_al[l_s[1][0]]);
387 al=declare_char (2, l0+l1+1);
388 B=realloc_aln2 (B,B->nseq,l0+l1+1);
391 free_int (B->cdna_cache, -1);
392 B->cdna_cache=declare_int(1, l0+l1+1);
394 if ( !M)M=initialize_dna_dp_model (CL);
399 tl_s=declare_int (2, 2);tns=vcalloc(2, sizeof(int));tl_s[0][0]=0;tl_s[1][0]=3;tns[0]=tns[1]=1;
400 DPR=make_fast_dp_pair_wise (A,tns, tl_s,CL,M);
401 vfree(tns);free_int(tl_s, -1);
408 while ( (k=DPR->traceback[a++])!=M->START);
409 while ( (k=DPR->traceback[a++])!=M->END)
412 f0=M->model_properties[k][M->F0];
413 f1=M->model_properties[k][M->F1];
415 len_i=M->model_properties[k][M->LEN_I];
416 len_j=M->model_properties[k][M->LEN_J];
418 type=M->model_properties[k][M->TYPE];
422 if (type==M->CODING0)
424 deltaf0=(aa0*3+f0)-ala;
425 deltaf1=(aa1*3+f1)-alb;
427 delta=MAX(deltaf0, deltaf1);
429 for (nr1=0, nr2=0,c=0; c<delta; c++, nr1++, nr2++,p++)
431 if (nr1<deltaf0 && ala<l0)al[0][p]=B->seq_al[l_s[0][0]][ala++];
434 if (nr2<deltaf1 && alb<l1)al[1][p]=B->seq_al[l_s[1][0]][alb++];
437 B->cdna_cache[0][p]=M->NON_CODING;
438 if ( is_gap(al[1][p]) && is_gap(al[0][p]))p--;
439 else if ( debug_cdna_fasta)fprintf (stderr, "\nUM: %c %c", al[0][p], al[1][p]);
441 for ( c=0; c< 3; c++, p++)
443 if ( c==0)B->cdna_cache[0][p]=M->CODING0;
444 else if ( c==1)B->cdna_cache[0][p]=M->CODING1;
445 else if ( c==2)B->cdna_cache[0][p]=M->CODING2;
446 if (ala<l0)al[0][p]=B->seq_al[l_s[0][0]][ala++];
449 if (alb<l1)al[1][p]=B->seq_al[l_s[1][0]][alb++];
452 if ( is_gap(al[1][p]) && is_gap(al[0][p]))p--;
453 else if ( debug_cdna_fasta)fprintf (stderr, "\n%d: %c %c",k, al[0][p], al[1][p]);
461 deltaf0=(aa0*3+f0)-ala;
462 deltaf1=(aa1*3+f1)-alb;
463 delta=MAX(deltaf0, deltaf1);
464 for (nr1=0, nr2=0,c=0; c<delta; c++, nr1++, nr2++,p++)
466 if (nr1<deltaf0 && ala<l0)al[0][p]=B->seq_al[l_s[0][0]][ala++];
469 if (nr2<deltaf1 && alb<l1)al[1][p]=B->seq_al[l_s[1][0]][alb++];
472 B->cdna_cache[0][p]=M->NON_CODING;
473 if ( is_gap(al[1][p]) && is_gap(al[0][p]))p--;
474 else if ( debug_cdna_fasta)fprintf (stderr, "\nUM: %c %c", al[0][p], al[1][p]);
478 /*End New traceback*/
487 sprintf( B->seq_al[l_s[0][0]], "%s", al[0]);
488 sprintf( B->seq_al[l_s[1][0]], "%s", al[1]);
489 B->len_aln=strlen (al[0]);
495 if ( debug_cdna_fasta)
497 fprintf ( stderr, "\nA-A=%d, %d", CL->M['a'-'A']['a'-'A'], CL->M['a'-'A']['a'-'A'] *SCORE_K);
498 for ( a=1; a<diag[0]; a++)
500 fprintf ( stderr, "\nchosen diag: %d", diag[a]);
503 fprintf ( stderr, "\n GOP=%d GEP=%d TG_MODE=%d", M->gop, M->gep, M->TG_MODE);
504 fprintf ( stderr, "\nF_GOP=%d F_GEP=%d F_TG_MODE=%d", M->gop, M->gep, M->F_TG_MODE);
506 DA=copy_aln (B, NULL);
507 DA=realloc_aln2 (DA,6,(DA->len_aln+1));
510 for ( a=0; a<B->len_aln; a++)
513 fprintf ( stderr, "\n%d", DA->cdna_cache[0][a]);
514 if (DA->cdna_cache[0][a]>=M->CODING0)DA->seq_al[DA->nseq][a]=DA->cdna_cache[0][a]-M->nstate+'0';
515 else DA->seq_al[DA->nseq][a]=DA->cdna_cache[0][a]-M->nstate+'0';
517 if (DA->cdna_cache[0][a]==M->CODING0)
519 DA->seq_al[DA->nseq+1][a]=translate_dna_codon (DA->seq_al[0]+a,'*');
520 DA->seq_al[DA->nseq+2][a]=translate_dna_codon (DA->seq_al[1]+a,'*');
524 DA->seq_al[DA->nseq+1][a]='-';
525 DA->seq_al[DA->nseq+2][a]='-';
536 for (prev_state=M->START,a=0; a< DA->len_aln;)
538 state=DA->cdna_cache[0][a];
539 t=M->model[prev_state][state];
540 if ( DA->cdna_cache[0][a]==M->CODING0)
542 a1=translate_dna_codon (A->seq_al[0]+a,'x');
543 a2=translate_dna_codon (A->seq_al[1]+a,'x');
545 if ( a1!='x' && a2!='x')
547 e=CL->M[a1-'A'][a2-'A']*SCORE_K;
550 else if ( DA->cdna_cache[0][a]>M->CODING0);
553 e=M->model_properties[B->cdna_cache[0][a]][M->EMISSION];
555 if ( e==UNDEFINED || t==UNDEFINED) fprintf ( stderr, "\nPROBLEM %d\n", a);
557 fprintf ( stderr, "\n[%c..%c: %d(e)+%d(t)=%d]", A->seq_al[0][a], A->seq_al[1][a], e,t,e+t);
561 if (B->cdna_cache[0][a]==M->NON_CODING)a++;
568 for ( a=0; a<B->len_aln; a++)
571 if ( B->cdna_cache[0][a]<M->CODING0)B->cdna_cache[0][a]=0;
572 else B->cdna_cache[0][a]=1;
582 Dp_Result * make_fast_dp_pair_wise (Alignment *A,int*ns, int **l_s, Constraint_list *CL,Dp_Model *M)
588 int l0, l1, len_al,len_diag;
589 static int max_len_al, max_len_diag;
598 static int *Mat, *LMat, *trace;
600 int state, cur_state, prev_state;
602 int last_i=0, last_j=0;
604 int len_i, len_j, len;
620 l0=strlen (A->seq_al[l_s[0][0]]);
621 l1=strlen (A->seq_al[l_s[1][0]]);
625 if ( (len_al>max_len_al || len_diag>max_len_diag))
631 max_len_diag=max_len_al=0;
637 max_len_diag=len_diag;
638 mI=max_len_al*max_len_diag;
642 Mat =vcalloc ( M->nstate*max_len_al*max_len_diag, sizeof (int));
643 LMat =vcalloc ( M->nstate*max_len_al*max_len_diag, sizeof (int));
644 trace=vcalloc ( M->nstate*max_len_al*max_len_diag, sizeof (int));
648 prev=vcalloc ( M->nstate, sizeof (int));
649 DPR=vcalloc ( 1, sizeof ( Dp_Result));
650 DPR->traceback=vcalloc (max_len_al, sizeof (int));
652 /*PREPARE THE EVALUATION*/
655 fprintf ( stderr, "\nERROR: function make_fasta_cdna_pair_wise can only handle two sequences at a time [FATAL:%s]",PROGRAM);
660 /*INITIALIZATION OF THE DP MATRICES*/
664 for (j=0; j<=ndiag;j++)
666 for ( state=0; state<M->nstate; state++)
668 Mat [state*mI+i*mJ+j]=UNDEFINED;
669 LMat [state*mI+i*mJ+j]=UNDEFINED;
670 trace [state*mI+i*mJ+j]=M->START;
677 for (i=0; i<=l0; i++)
678 for ( j=0; j<=ndiag; j++)
680 pos_j=M->diag[j]-l0+i;
682 if (!(pos_j==0 || pos_i==0))continue;
683 if ( pos_j<0 || pos_i<0)continue;
684 if ( pos_i==0 && pos_j==0)
686 for ( a=0; a< M->nstate; a++)
689 LMat [a*mI+i*mJ+j]=0;
690 trace[a*mI+i*mJ+j]=M->START;
696 for ( state=0; state<M->START; state++)
698 if (pos_j==0 && M->model_properties[state][M->LEN_J])continue;
699 if (pos_i==0 && M->model_properties[state][M->LEN_I])continue;
702 t=M->model[M->START][state];
703 e=M->model_properties[state][M->TERM_EMISSION];
704 Mat [state*mI+i*mJ+j]=t+e*l;
705 LMat [state*mI+i*mJ+j]=l;
706 trace [state*mI+i*mJ+j]=M->START;
711 /*DYNAMIC PROGRAMMING: Forward Pass*/
717 for (j=1; j<=ndiag;j++)
719 pos_j=M->diag[j]-l0+i;
722 if (pos_j<=0 || pos_j>l1 )continue;
726 for (cur_state=0; cur_state<M->START; cur_state++)
728 if (M->model_properties[cur_state][M->DELTA_J])
730 prev_j=j+M->model_properties[cur_state][M->DELTA_J];
731 prev_i=i+M->model_properties[cur_state][M->DELTA_I]*FABS((M->diag[j]-M->diag[prev_j]));
736 prev_i=i+M->model_properties[cur_state][M->DELTA_I];
738 len_i=FABS((i-prev_i));
739 len_j=FABS((M->diag[prev_j]-M->diag[j]));
740 len=MAX(len_i, len_j);
741 a1=A->seq_al[M->model_properties[cur_state][M->F0] ][pos_i-1];
742 a2=A->seq_al[M->model_properties[cur_state][M->F1]+3][pos_j-1];
744 if (M->model_properties[cur_state][M->TYPE]==M->CODING0)
746 if ( a1=='o' || a2=='o')em=-(mat['w'-'A']['w'-'A'])*SCORE_K;
747 else if (a1=='x' || a2=='x')em=UNDEFINED;
748 else if ( a1==0 || a2==0)exit (0);
751 em=(mat[a1-'A'][a2-'A'])*SCORE_K;
756 em=M->model_properties[cur_state][M->EMISSION];
761 for (pc=best_pc=UNDEFINED, model_index=1; model_index<=M->bounded_model[cur_state][0]; model_index++)
763 prev_state=M->bounded_model[cur_state][model_index];
765 if(prev_i<0 || prev_j<0 ||prev_i>l0 || prev_j>ndiag || len==UNDEFINED)prev_score=UNDEFINED;
766 else prev_score=Mat[prev_state*mI+prev_i*mJ+prev_j];
767 t=M->model[prev_state][cur_state];
770 if (prev_score==UNDEFINED || len==UNDEFINED)e=UNDEFINED;
771 else if (len==0|| e==UNDEFINED)e=UNDEFINED;
774 if (is_defined_int(3,prev_score,e, t))
780 /*Identify the best previous score*/
781 if (best_pc==UNDEFINED || (pc>best_pc && pc!=UNDEFINED))
783 prev[cur_state]=prev_state;
789 Mat[cur_state*mI+i*mJ+j]=best_pc;
793 if ( Mat[cur_state*mI+i*mJ+j]==UNDEFINED)
795 LMat[cur_state*mI+i*mJ+j]=UNDEFINED;
796 trace[cur_state*mI+i*mJ+j]=UNDEFINED;
800 else if ( prev[cur_state]==cur_state)
802 LMat [cur_state*mI+i*mJ+j]= LMat [cur_state*mI+prev_i*mJ+prev_j]+len;
803 trace[cur_state*mI+i*mJ+j]= trace[cur_state*mI+prev_i*mJ+prev_j];
807 LMat[cur_state*mI+i*mJ+j]=len;
808 trace[cur_state*mI+i*mJ+j]=prev[cur_state];
817 for (pc=best_pc=UNDEFINED, state=0; state<M->START; state++)
819 t=M->model[state][M->END];
820 e=M->model_properties[state][M->TERM_EMISSION];
821 l=LMat[state*mI+i*mJ+j];
824 if (!is_defined_int(4,t,e,Mat[state*mI+i*mJ+j],l))Mat[state*mI+i*mJ+j]=UNDEFINED;
825 else Mat[state*mI+i*mJ+j]+=t+e*(l);
826 pc=Mat[state*mI+i*mJ+j];
829 if (best_pc==UNDEFINED || (pc>best_pc && pc!=UNDEFINED))
846 next_k=trace[k*mI+i*mJ+j];
853 DPR->traceback[len++]=k;
855 new_i+=M->model_properties[k][M->DELTA_I]*l;
858 if ( M->model_properties[k][M->DELTA_J])
860 while ( next_k!=M->START && FABS((M->diag[j]-M->diag[new_j]))!=l)new_j+=M->model_properties[k][M->DELTA_J];
868 DPR->traceback[DPR->len++]=M->START;
869 invert_list_int (DPR->traceback,DPR->len);
870 DPR->traceback[DPR->len]=M->END;
881 int ** evaluate_diagonals_cdna ( Alignment *B, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup)
896 if ( ns[0]!=1 || ns[1]!=1)
898 fprintf ( stderr, "\nERROR 2 SEQUENCES ONLY [FATAL:%s", PROGRAM);
906 l0=strlen ( B->seq_al[0]);
907 l1=strlen ( B->seq_al[3]);
908 n_tot_diag=(l0+l1-1);
910 tot_diag=declare_int ( n_tot_diag+1, 2);
911 for ( c=0; c<= n_tot_diag; c++)tot_diag[c][0]=c;
913 for (f1=0; f1< 3; f1++)
915 for ( f2=0; f2< 3; f2++)
925 diag=evaluate_diagonals_for_two_sequences( s1, s2, maximise,NULL,ktup);
926 for (c=1; c<=(p1+p2-1); c++)
928 tot_diag[diag[c][0]][1]+=diag[c][1]*diag[c][1];
937 sort_int (tot_diag+1, 2, 1,0, n_tot_diag-1);
946 /*********************************COPYRIGHT NOTICE**********************************/
947 /*© Centro de Regulacio Genomica */
949 /*Cedric Notredame */
950 /*Tue Oct 27 10:12:26 WEST 2009. */
951 /*All rights reserved.*/
952 /*This file is part of T-COFFEE.*/
954 /* T-COFFEE is free software; you can redistribute it and/or modify*/
955 /* it under the terms of the GNU General Public License as published by*/
956 /* the Free Software Foundation; either version 2 of the License, or*/
957 /* (at your option) any later version.*/
959 /* T-COFFEE is distributed in the hope that it will be useful,*/
960 /* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
961 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
962 /* GNU General Public License for more details.*/
964 /* You should have received a copy of the GNU General Public License*/
965 /* along with Foobar; if not, write to the Free Software*/
966 /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
967 /*............................................... |*/
968 /* If you need some more information*/
969 /* cedric.notredame@europe.com*/
970 /*............................................... |*/
974 /*********************************COPYRIGHT NOTICE**********************************/