7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
11 #include "dp_lib_header.h"
12 float compute_lambda (int **matrix,char *alphabet);
13 /*********************************************************************************************/
15 /* FUNCTIONS FOR EVALUATING THE CONSISTENCY BETWEEN ALN AND CL */
17 /*********************************************************************************************/
19 /*Fast: score= extended_res/max_extended_residue_for the whole aln
20 slow: score= extended_res/sum all extended score for that residue
21 non_extended score= non_ext /sum all non extended score for that residue
22 heuristic score= extended /sum of extended score of all pairs in the library
23 (i.e. Not ALL the possible pairs)
25 Alignment * main_coffee_evaluate_output2 ( Alignment *IN,Constraint_list *CL, const char *mode );
26 int sub_aln2ecl_raw_score (Alignment *A, Constraint_list *CL, int ns, int *ls)
29 int p1,r1, r2, s1, s2;
33 A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
34 pos=aln2pos_simple ( A,A->nseq);
38 for (p1=0; p1<A->len_aln; p1++)
40 for (s1=0; s1<ns-1; s1++)
42 for (s2=s1+1; s2<ns; s2++)
49 score+= residue_pair_extended_list_pc (CL,ls[s1], r1, ls[s2], r2)*SCORE_K;
56 return (score/(((ns*ns)-ns)/2))/A->len_aln;
58 int aln2ecl_raw_score (Alignment *A, Constraint_list *CL)
61 int p1,r1, r2, s1, s2;
64 A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
65 pos=aln2pos_simple ( A,A->nseq);
69 for (p1=0; p1<A->len_aln; p1++)
71 for (s1=0; s1<A->nseq-1; s1++)
73 for (s2=s1+1; s2<A->nseq; s2++)
78 if (r1>0 && r2>0)score+= residue_pair_extended_list_pc (CL,s1, r1, s2, r2);
84 return (score/(((A->nseq*A->nseq)-A->nseq)/2))/A->len_aln;
86 int node2sub_aln_score (Alignment *A,Constraint_list *CL, char *mode, NT_node T)
88 if ( !T || !T->right ||!T->left)return 0;
93 ns=vcalloc (2, sizeof (int));
94 ls=vcalloc (2, sizeof (int*));
95 ns[0]= (T->left)->nseq;
96 ns[1]=(T->right)->nseq;
97 ls[0]= (T->left)->lseq;
98 ls[1]=(T->right)->lseq;
100 return sub_aln2sub_aln_score (A, CL, mode, ns, ls);
104 int sub_aln2sub_aln_score ( Alignment *A,Constraint_list *CL, const char *mode, int *ns, int **ls)
106 /*Warning: Does Not Take Gaps into account*/
112 /*Make sure evaluation functions update their cache if needed*/
113 A=update_aln_random_tag (A);
114 pos=aln2pos_simple ( A, -1, ns, ls);
115 for (a=0; a< A->len_aln; a++)
116 score+=CL->get_dp_cost (A, pos, ns[0], ls[0], a, pos, ns[1],ls[1], a, CL);
119 score=(int)(((float)score)/(A->len_aln*SCORE_K));
120 score=(int)(CL->residue_index && CL->normalise)?((score*MAXID)/(CL->normalise)):(score);
123 int sub_aln2sub_aln_raw_score ( Alignment *A,Constraint_list *CL, const char *mode, int *ns, int **ls)
125 /*Warning: Does Not Take Gaps into account*/
131 /*Make sure evaluation functions update their cache if needed*/
132 A=update_aln_random_tag (A);
133 pos=aln2pos_simple ( A, -1, ns, ls);
134 for (a=0; a< A->len_aln; a++)
135 score+=CL->get_dp_cost (A, pos, ns[0], ls[0], a, pos, ns[1],ls[1], a, CL);
140 Alignment* main_coffee_evaluate_output_sub_aln ( Alignment *A,Constraint_list *CL, const char *mode, int *n_s, int **l_s)
142 Alignment *SUB1, *SUB2, *SUB3;
143 int a, b, c,*list_seq;
146 if (strm ( CL->evaluate_mode, "no"))return NULL;
149 list_seq=vcalloc (n_s[0]+n_s[1], sizeof (int));
150 for (b=0, a=0; a< 2; a++){for (c=0;c< n_s[a]; c++)list_seq[b++]=l_s[a][c];}
153 SUB1=copy_aln (A, NULL);
154 SUB2=extract_sub_aln (SUB1,n_s[0]+n_s[1],list_seq);
155 SUB3=main_coffee_evaluate_output (SUB2,CL,CL->evaluate_mode);
163 Alignment * overlay_alignment_evaluation ( Alignment *I, Alignment *O)
168 if ( !I || !O) return O;
169 if ( I->len_aln!=O->len_aln)printf_exit (EXIT_FAILURE, stderr, "ERROR: Incompatible alignments in overlay_alignment_evaluation");
171 buf=vcalloc ( MAX(I->len_aln, O->len_aln), sizeof (int));
173 for (a=0; a<O->nseq; a++)
175 if (!strm (I->name[a], O->name[a]))printf_exit (EXIT_FAILURE, stderr, "ERROR: Incompatible alignments in overlay_alignment_evaluation");
176 for (b=0; b<O->len_aln; b++)
179 if ( islower(r))O->seq_al[a][b]=0;
180 else if (r<=9 || (r>='0' && r<='9'))O->seq_al[a][b]=I->seq_al[a][b];
186 Alignment * main_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL, const char *mode )
189 Alignment *TopS=NULL, *LastS=NULL, *CurrentS=NULL;
196 CurrentS= main_coffee_evaluate_output2(IN, CL, mode);
197 if (!TopS)LastS=TopS=CurrentS;
208 Alignment * main_coffee_evaluate_output2 ( Alignment *IN,Constraint_list *CL, const char *mode )
211 /*Make sure evaluation functions update their cache if needed*/
212 IN=update_aln_random_tag (IN);
214 if ( CL->evaluate_residue_pair==evaluate_matrix_score || CL->ne==0 ||strm ( mode , "categories") || strm ( mode , "matrix")|| strm(mode, "sar")|| strstr (mode, "boxshade") )
217 if ( strm ( mode , "categories")) return categories_evaluate_output (IN, CL);
218 else if ( strm ( mode , "matrix"))return matrix_evaluate_output (IN, CL);
219 else if ( strm ( mode, "sar"))return sar_evaluate_output (IN, CL);
220 else if ( strstr ( mode, "boxshade"))return boxshade_evaluate_output (IN, CL, atoi (strstr(mode, "_")+1));
222 else if ( CL->evaluate_residue_pair==evaluate_matrix_score) return matrix_evaluate_output (IN, CL);
223 else if ( CL->ne==0) return matrix_evaluate_output (IN, CL);
225 else if ( strm (mode, "no"))return NULL;
226 else if ( strstr( mode, "triplet"))
228 return triplet_coffee_evaluate_output ( IN,CL);
230 else if ( strstr ( mode, "fast"))
232 return fast_coffee_evaluate_output ( IN,CL);
234 else if ( strstr ( mode, "slow"))
236 return slow_coffee_evaluate_output ( IN,CL);
239 else if ( strstr ( mode, "non_extended"))
241 return non_extended_t_coffee_evaluate_output ( IN,CL);
244 else if ( strm (mode, "sequences"))
246 return coffee_seq_evaluate_output ( IN,CL);
250 fprintf ( stderr, "\nUNKNOWN MODE FOR ALIGNMENT EVALUATION: *%s* [FATAL:%s]",mode, PROGRAM);
259 Alignment * coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
261 fprintf ( stderr, "\n[WARNING:%s]THE FUNCTION coffee_evaluate_output IS NOT ANYMORE SUPPORTED\n", PROGRAM);
262 fprintf ( stderr, "\n[WARNING]fast_coffee_evaluate_output WILL BE USED INSTEAD\n");
264 return fast_coffee_evaluate_output (IN,CL);
266 Alignment * matrix_evaluate_output ( Alignment *IN,Constraint_list *CL)
268 int a,b, c,r, s, r1, r2;
285 Residue x: sum of observed extended X.. /sum of possible X..
289 if ( !CL->M)CL->M=read_matrice ("blosum62mt");
291 OUT=copy_aln (IN, OUT);
294 tot_res=declare_double ( IN->nseq, IN->len_aln);
295 max_res=declare_double ( IN->nseq, IN->len_aln);
297 tot_seq=declare_double ( IN->nseq, 1);
298 max_seq=declare_double ( IN->nseq, 1);
299 tot_col=declare_double ( IN->len_aln,1);
300 max_col=declare_double ( IN->len_aln,1);
304 for (a=0; a< IN->len_aln; a++)
306 for ( b=0; b< IN->nseq; b++)
308 r1=tolower(IN->seq_al[b][a]);
309 if ( is_gap(r1))continue;
310 r= CL->M[r1-'A'][r1-'A'];
312 for ( c=0; c<IN->nseq; c++)
314 r2=tolower(IN->seq_al[c][a]);
315 if (b==c || is_gap (r2))continue;
317 s=CL->M[r2-'A'][r1-'A'];
336 for ( a=0; a< IN->nseq; a++)
338 if ( !max_seq[a][0])continue;
340 OUT->score_seq[a]=(tot_seq[a][0]*100)/max_seq[a][0];
341 for (b=0; b< IN->len_aln; b++)
344 if ( is_gap(r1) || !max_res[a][b])continue;
346 r1=(tot_res[a][b]*10)/max_res[a][b];
349 (OUT)->seq_al[a][b]=r1+'0';
353 for ( a=0; a< IN->len_aln; a++)
355 r1=(max_col[a][0]==0)?0:((tot_col[a][0]*10)/max_col[a][0]);
357 (OUT)->seq_al[OUT->nseq][a]=r1+'0';
359 sprintf ( OUT->name[IN->nseq], "cons");
360 if (max_aln)OUT->score_seq[OUT->nseq]=OUT->score_aln=(100*tot_aln)/max_aln;
363 free_double (tot_res,-1);
364 free_double (max_res,-1);
366 free_double (tot_seq,-1);
367 free_double (max_seq,-1);
372 Alignment * sar_evaluate_output ( Alignment *IN,Constraint_list *CL)
374 int a,b, c,r, s, r1, r2;
391 Residue x: sum of observed extended X.. /sum of possible X..
395 if ( !CL->M)CL->M=read_matrice ("blosum62mt");
397 OUT=copy_aln (IN, OUT);
400 tot_res=declare_double ( IN->nseq, IN->len_aln);
401 max_res=declare_double ( IN->nseq, IN->len_aln);
403 tot_seq=declare_double ( IN->nseq, 1);
404 max_seq=declare_double ( IN->nseq, 1);
405 tot_col=declare_double ( IN->len_aln,1);
406 max_col=declare_double ( IN->len_aln,1);
410 for (a=0; a< IN->len_aln; a++)
412 for (b=0; b< IN->nseq; b++)
414 r1=tolower(IN->seq_al[b][a]);
415 for (c=0; c<IN->nseq; c++)
417 r2=tolower(IN->seq_al[c][a]);
420 if ( is_gap(r1) && is_gap(r2))s=0;
441 for ( a=0; a< IN->nseq; a++)
443 if ( !max_seq[a][0])continue;
444 OUT->score_seq[a]=(max_seq[a][0]*100)/max_seq[a][0];
445 for (b=0; b< IN->len_aln; b++)
448 if ( is_gap(r1) || !max_res[a][b])continue;
449 r1=(tot_res[a][b]*10)/max_res[a][b];
452 (OUT)->seq_al[a][b]=r1+'0';
456 for ( a=0; a< IN->len_aln; a++)
458 r1=(max_col[a][0]==0)?0:((tot_col[a][0]*10)/max_col[a][0]);
460 (OUT)->seq_al[OUT->nseq][a]=r1+'0';
462 sprintf ( OUT->name[IN->nseq], "cons");
463 if (max_aln)OUT->score_aln=(100*tot_aln)/max_aln;
466 free_double (tot_res,-1);
467 free_double (max_res,-1);
469 free_double (tot_seq,-1);
470 free_double (max_seq,-1);
474 Alignment * boxshade_evaluate_output ( Alignment *IN,Constraint_list *CL, int T)
483 Residue x: sum of observed extended X.. /sum of possible X..
487 OUT=copy_aln (IN, OUT);
488 aa=declare_int (26, 2);
490 for ( a=0; a< OUT->len_aln; a++)
492 for ( b=0; b< 26; b++){aa[b][1]=0;aa[b][0]='a'+b;}
493 for ( b=0; b< OUT->nseq; b++)
495 r=tolower(OUT->seq_al[b][a]);
496 if ( !is_gap(r))aa[r-'a'][1]++;
498 sort_int ( aa, 2, 1, 0,25);
499 f=(aa[25][1]*100)/OUT->nseq;
504 bs=((f-T)*10)/(100-T);
507 if (bs==10)bs--;bs+='0';
508 for ( b=0; b< OUT->nseq; b++)
510 r=tolower(OUT->seq_al[b][a]);
511 if (r==br && bs>'1')OUT->seq_al[b][a]=bs;
513 OUT->seq_al[b][a]=bs;
516 sprintf ( OUT->name[IN->nseq], "cons");
521 Alignment * categories_evaluate_output ( Alignment *IN,Constraint_list *CL)
527 float score, nseq2, tot_aln;
530 Residue x: sum of observed extended X.. /sum of possible X..
532 OUT=copy_aln (IN, OUT);
533 aa=vcalloc ( 26, sizeof (int));
534 nseq2=IN->nseq*IN->nseq;
536 for (tot_aln=0, a=0; a< IN->len_aln; a++)
538 for (n=0,b=0; b< IN->nseq; b++)
545 aa[tolower(r)-'a']++;
550 for ( score=0,b=0; b< 26; b++){score+=aa[b]*aa[b];aa[b]=0;}
552 score=(n==0)?0:score/n;
556 (OUT)->seq_al[OUT->nseq][a]='0'+r;
558 OUT->score_aln=(tot_aln/OUT->len_aln)*100;
559 sprintf ( OUT->name[IN->nseq], "cons");
564 Alignment * categories_evaluate_output_old ( Alignment *IN,Constraint_list *CL)
570 float score, nseq2, tot_aln, min=0;
573 Residue x: sum of observed extended X.. /sum of possible X..
575 OUT=copy_aln (IN, OUT);
576 aa=vcalloc ( 26, sizeof (int));
577 nseq2=IN->nseq*IN->nseq;
579 for (tot_aln=0, a=0; a< IN->len_aln; a++)
581 for (ng=0,b=0; b< IN->nseq; b++)
588 aa[tolower(r)-'a']++;
591 for (nc=0, b=0; b<26; b++)
599 score=(2*min)/IN->nseq;
604 (OUT)->seq_al[OUT->nseq][a]='0'+r;
607 OUT->score_aln=(tot_aln/OUT->len_aln)*100;
608 sprintf ( OUT->name[IN->nseq], "cons");
612 Alignment * fork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL);
613 Alignment * nfork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL);
614 Alignment * triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
617 if (!IN || !CL || !CL->residue_index) return IN;
620 if ( get_nproc()==1)return nfork_triplet_coffee_evaluate_output (IN,CL);
621 else if (strstr ( CL->multi_thread, "evaluate"))return fork_triplet_coffee_evaluate_output (IN,CL);
622 else return nfork_triplet_coffee_evaluate_output (IN,CL);
624 Alignment * fork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
629 double score_col=0, score_aln=0, score_res=0, score=0;
630 double max_col=0, max_aln=0, max_res=0;
631 double *max_seq, *score_seq;
634 int s1,r1,s2,r2,w2,s3,r3,w3;
644 OUT=copy_aln (IN, OUT);
645 pos=aln2pos_simple(IN, IN->nseq);
646 sprintf ( OUT->name[IN->nseq], "cons");
648 max_seq=vcalloc ( IN->nseq, sizeof (double));
649 score_seq=vcalloc ( IN->nseq, sizeof (double));
650 lu=declare_int (IN->nseq, IN->len_aln+1);
652 //multi Threading stuff
654 sl=n2splits (njobs,IN->len_aln);
655 pid_tmpfile=vcalloc (njobs, sizeof (char*));
657 for (sjobs=0,j=0; sjobs<njobs; j++)
659 pid_tmpfile[j]=vtmpnam(NULL);
660 if (vvfork (NULL)==0)
662 initiate_vtmpnam(NULL);
663 fp=vfopen (pid_tmpfile[j], "w");
665 for (a=sl[j][0]; a<sl[j][1]; a++)
667 if (j==0)output_completion (stderr,a,sl[0][1],1, "Final Evaluation");
669 for (b=0; b<IN->nseq; b++)
673 if (r1>=0)lu[s1][r1]=1;
675 for (b=0; b<IN->nseq; b++)
677 score_res=max_res=NORM_F;
678 res=NO_COLOR_RESIDUE;
683 fprintf (fp, "%d ", res);
687 for (x=1;x<CL->residue_index[s1][r1][0];x+=ICHUNK)
689 s2=CL->residue_index[s1][r1][x+SEQ2];
690 r2=CL->residue_index[s1][r1][x+R2];
691 w2=CL->residue_index[s1][r1][x+WE];
692 for (y=1; y<CL->residue_index[s2][r2][0];y+=ICHUNK)
694 s3=CL->residue_index[s2][r2][y+SEQ2];
695 r3=CL->residue_index[s2][r2][y+R2];
696 w3=CL->residue_index[s2][r2][y+WE];
714 res=(max_res==0)?NO_COLOR_RESIDUE:((score_res*10)/max_res);
715 res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
716 fprintf ( fp, "%d ", res);
718 for (b=0; b<IN->nseq; b++)
722 if (r1>0)lu[s1][r1]=0;
725 res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
726 res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
727 fprintf (fp, "%d ", res);
728 for (b=0; b<IN->nseq; b++)fprintf (fp, "%f %f ", score_seq[b], max_seq[b]);
730 fprintf (fp, "%f %f ", score_aln, max_aln);
732 myexit (EXIT_SUCCESS);
740 while (sjobs>=0){vwait(NULL); sjobs--;}//wait for all jobs to complete
742 for (j=0; j<njobs; j++)
745 fp=vfopen (pid_tmpfile[j], "r");
746 for (a=sl[j][0];a<sl[j][1]; a++)
748 for (b=0; b<=IN->nseq; b++)//don't forget the consensus
750 fscanf (fp, "%d ", &res);
751 OUT->seq_al[b][a]=res;
753 for (b=0; b<IN->nseq; b++)
755 fscanf (fp, "%f %f", &sseq, &mseq);
760 fscanf (fp, "%f %f", &sseq, &mseq);
764 remove (pid_tmpfile[j]);
766 fprintf ( stderr, "\n");
768 IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
769 for ( a=0; a< OUT->nseq; a++)
771 OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
783 Alignment * nfork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
788 double score_col=0, score_aln=0, score_res=0, score=0;
789 double max_col=0, max_aln=0, max_res=0;
790 double *max_seq, *score_seq;
793 int s1,r1,s2,r2,w2,s3,r3,w3;
798 OUT=copy_aln (IN, OUT);
799 pos=aln2pos_simple(IN, IN->nseq);
800 sprintf ( OUT->name[IN->nseq], "cons");
802 max_seq=vcalloc ( IN->nseq, sizeof (double));
803 score_seq=vcalloc ( IN->nseq, sizeof (double));
804 lu=declare_int (IN->nseq, IN->len_aln+1);
810 for (a=0; a<IN->len_aln; a++)
812 output_completion (stderr,a,IN->len_aln,1, "Final Evaluation");
814 for (b=0; b<IN->nseq; b++)
818 if (r1>=0)lu[s1][r1]=1;
820 for (b=0; b<IN->nseq; b++)
822 score_res=max_res=NORM_F;
823 OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
827 for (x=1;x<CL->residue_index[s1][r1][0];x+=ICHUNK)
829 s2=CL->residue_index[s1][r1][x+SEQ2];
830 r2=CL->residue_index[s1][r1][x+R2];
831 w2=CL->residue_index[s1][r1][x+WE];
832 for (y=1; y<CL->residue_index[s2][r2][0];y+=ICHUNK)
834 s3=CL->residue_index[s2][r2][y+SEQ2];
835 r3=CL->residue_index[s2][r2][y+R2];
836 w3=CL->residue_index[s2][r2][y+WE];
854 res=(max_res==0)?NO_COLOR_RESIDUE:((score_res*10)/max_res);
855 res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
856 OUT->seq_al[b][a]=res;
858 for (b=0; b<IN->nseq; b++)
862 if (r1>0)lu[s1][r1]=0;
865 res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
866 res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
867 OUT->seq_al[IN->nseq][a]=res;
869 fprintf ( stderr, "\n");
871 IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
872 for ( a=0; a< OUT->nseq; a++)
874 OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
883 Alignment * fast_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
885 int a,b, c, m,res, s, s1, s2, r1, r2;
889 double score_col=0, score_aln=0, score_res=0;
890 double max_col, max_aln;
891 double *max_seq, *score_seq;
896 /*NORMALIZE: with the highest scoring pair found in the multiple alignment*/
899 if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
901 OUT=copy_aln (IN, OUT);
902 pos=aln2pos_simple(IN, IN->nseq);
903 pos2=aln2defined_residues (IN, CL);
905 max_seq=vcalloc ( IN->nseq, sizeof (double));
906 score_seq=vcalloc ( IN->nseq, sizeof (double));
910 /*1: Identify the highest scoring pair within the alignment*/
912 for ( m=0, a=0; a< IN->len_aln; a++)
914 for ( b=0; b< IN->nseq; b++)
920 for ( c=0; c< IN->nseq; c++)
924 if ( s1==s2 && !CL->do_self)continue;
926 if ( s1< s2)s=(CL->evaluate_residue_pair)( CL, s1, r1, s2, r2);
927 else s=(CL->evaluate_residue_pair)( CL, s2, r2, s1, r1);
929 s=(s!=UNDEFINED)?s:0;
937 sprintf ( OUT->name[IN->nseq], "cons");
938 for ( max_aln=0,score_aln=0,a=0; a< IN->len_aln; a++)
940 OUT->seq_al[IN->nseq][a]=NO_COLOR_RESIDUE;
941 for ( local_nseq=0,b=0; b<IN->nseq; b++){local_nseq+=(pos[b][a]>0 && pos2[b][a])?1:0;}
942 local_m=m*(local_nseq-1);
944 for ( max_col=0, score_col=0,b=0; b< IN->nseq; b++)
946 OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
950 if (r1<=0 || !pos2[b][a])
955 for ( score_res=0,c=0; c< IN->nseq; c++)
960 if ((s1==s2 && !CL->do_self) || r2<=0 || !pos2[c][a]){continue;}
965 if ( s1< s2)s=(CL->evaluate_residue_pair)( CL, s1, r1, s2, r2);
966 else s=(CL->evaluate_residue_pair)( CL, s2, r2, s1, r1);
967 s=(s!=UNDEFINED)?s:0;
975 res=(local_m==0)?NO_COLOR_RESIDUE:((score_res*10)/local_m);
976 (OUT)->seq_al[b][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
981 res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
982 OUT->seq_al[IN->nseq][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
986 IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
987 for ( a=0; a< OUT->nseq; a++)
989 OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
1000 Alignment * slow_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
1002 int a,b, c,res, s, s1, s2, r1, r2;
1003 Alignment *OUT=NULL;
1005 double max_score_r, score_r, max;
1006 double score_col=0, score_aln=0;
1007 double max_score_col, max_score_aln;
1008 double *max_score_seq, *score_seq;
1009 int ***res_extended_weight;
1014 Residue x: sum of observed extended X.. /sum of possible X..
1020 if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
1023 OUT=copy_aln (IN, OUT);
1024 pos=aln2pos_simple(IN, IN->nseq);
1025 pos2=aln2defined_residues (IN, CL);
1027 max_score_seq=vcalloc ( IN->nseq, sizeof (double));
1028 score_seq=vcalloc ( IN->nseq, sizeof (double));
1029 res_extended_weight=declare_arrayN(3,sizeof(int), (CL->S)->nseq, (CL->S)->max_len+1, 2);
1030 max=(CL->normalise)?(100*CL->normalise)*SCORE_K:100;
1032 for (a=0; a< IN->len_aln; a++)
1034 for ( b=0; b< IN->nseq-1; b++)
1038 for ( c=b+1; c< IN->nseq; c++)
1042 if ( s1==s2 && !CL->do_self)continue;
1043 else if ( r1<=0 || r2<=0) continue;
1046 s=(CL->evaluate_residue_pair)( CL, s1, r1, s2, r2);
1047 res_extended_weight[s1][r1][0]+=s*100;
1048 res_extended_weight[s2][r2][0]+=s*100;
1049 res_extended_weight[s1][r1][1]+=max;
1050 res_extended_weight[s2][r2][1]+=max;
1057 sprintf ( OUT->name[IN->nseq], "cons");
1058 for ( max_score_aln=0,score_aln=0,a=0; a< IN->len_aln; a++)
1060 OUT->seq_al[IN->nseq][a]=NO_COLOR_RESIDUE;
1061 for ( n_res_in_col=0,b=0; b<IN->nseq; b++){n_res_in_col+=(pos[b][a]>0 && pos2[b][a]>0)?1:0;}
1062 for ( max_score_col=0, score_col=0,b=0; b< IN->nseq; b++)
1064 OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
1067 if (r1<=0 || pos2[b][a]<1)continue;
1070 max_score_r =res_extended_weight[s1][r1][1];
1071 score_r =res_extended_weight[s1][r1][0];
1072 if ( max_score_r==0 && n_res_in_col>1)res=0;
1073 else if ( n_res_in_col==1)res=NO_COLOR_RESIDUE;
1074 else res=((score_r*10)/max_score_r);
1077 (OUT)->seq_al[b][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res, 9));
1078 max_score_col+=max_score_r;
1080 max_score_seq[b]+=max_score_r;
1081 score_seq[b]+=score_r;
1082 max_score_aln+=max_score_r;
1085 if ( max_score_col==0 && n_res_in_col>1)res=0;
1086 else if ( n_res_in_col<2)res=NO_COLOR_RESIDUE;
1087 else res=((score_col*10)/max_score_col);
1089 OUT->seq_al[IN->nseq][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
1092 IN->score_aln=OUT->score_aln=(max_score_aln==0)?0:((score_aln*100)/max_score_aln);
1093 for ( a=0; a< OUT->nseq; a++)
1095 OUT->score_seq[a]=(max_score_seq[a]==0)?0:((score_seq[a]*100)/max_score_seq[a]);
1100 vfree ( max_score_seq);
1101 free_arrayN((void*)res_extended_weight, 3);
1105 free_int (pos2, -1);
1109 double genepred2sd (Sequence *S);
1110 double genepred2avg (Sequence *S);
1111 double genepred2zsum2 (Sequence *S);
1112 int * genepred2orf_len (Sequence *S);
1113 double genepred2acc (Sequence *S, Sequence *PS);
1114 double seq_genepred2acc (Sequence *S, Sequence *PS, char *name);
1115 Alignment *coffee_seq_evaluate_output ( Alignment *IN, Constraint_list *CL)
1118 int avg, min_avg, best_cycle;
1119 char **pred_list,**weight_list;
1123 float nsd, best_score,score, acc;
1127 min_ne=CL->ne/10; //avoid side effects due to very relaxed libraries
1129 pred_list=declare_char (ncycle, 100);
1130 weight_list=declare_char(ncycle,100);
1132 fprintf ( stderr, "\nSCORE_MODE: %s\n", CL->genepred_score);
1133 for (a=0; a<ncycle && CL->ne>min_ne; a++)
1135 if (w);free_int (w, -1);
1136 w=list2residue_total_weight (CL);
1137 PS=dnaseq2geneseq (S, w);
1138 nsd=(float)(genepred2sd(PS)/genepred2avg(PS));
1139 if ( strm (CL->genepred_score, "nsd"))score=nsd;
1140 else score=1-seq_genepred2acc (S, PS, CL->genepred_score);
1141 fprintf ( stderr, "Cycle: %2d LIB: %6d AVG LENGTH: %6.2f NSD: %8.2f SCORE: %.2f ",a+1,CL->ne,(float) genepred2avg(PS),nsd,score);
1143 vfree (display_accuracy (genepred_seq2accuracy_counts (S, PS, NULL),stderr));
1144 pred_list[a]=vtmpnam(NULL);
1145 weight_list[a]=vtmpnam(NULL);
1146 output_fasta_seqS(pred_list[a],PS);
1147 output_array_int (weight_list[a],w);
1148 free_sequence (PS, -1);
1150 if (a==0 || score<best_score)
1155 CL=relax_constraint_list_4gp(CL);
1157 if (w)free_int(w,-1);
1158 S=get_fasta_sequence (pred_list[best_cycle], NULL);
1160 fprintf ( stderr, "\nSelected Cycle: %d (SCORE=%.4f)----> ", best_cycle+1, best_score);
1161 vfree (display_accuracy (genepred_seq2accuracy_counts ((CL->S), S, NULL),stderr));
1163 genepred_seq2accuracy_counts4all((CL->S), S);
1165 IN=seq2aln (S,NULL, 1);
1166 IN->score_res=input_array_int (weight_list[best_cycle]);
1170 double seq_genepred2acc (Sequence *S, Sequence *PS, char *name)
1177 if ( strm (name, "best"))return genepred2acc (S, PS);
1179 ref=name_is_in_list (name, S->name, S->nseq, 100);
1180 target=name_is_in_list (name, PS->name, PS->nseq, 100);
1182 if ( target==-1 || ref==-1)
1184 printf_exit (EXIT_FAILURE,stderr, "\nERROR: %s is not a valid sequence", name);
1186 count=genepred2accuracy_counts (S->seq[ref],PS->seq[target],NULL);
1187 acc=counts2accuracy (count);
1197 double genepred2acc (Sequence *S, Sequence *PS)
1202 count=genepred_seq2accuracy_counts (S, PS, NULL);
1203 acc=counts2accuracy (count);
1211 double genepred2sd (Sequence *S)
1213 double sum=0, sum2=0;
1216 len=genepred2orf_len (S);
1218 for (a=0; a<S->nseq; a++)
1220 sum+=(double)len[a];
1221 sum2+=(double)len[a]*(double)len[a];
1223 sum/=(double)S->nseq;
1226 return sqrt (sum2-(sum*sum));
1228 double genepred2avg (Sequence *S)
1233 len=genepred2orf_len (S);
1235 for (a=0; a<S->nseq; a++)avg+=len[a];
1237 return avg/(double)S->nseq;
1239 double genepred2zsum2 (Sequence *S)
1241 double sum=0, sum2=0, zscore=0, zsum2=0;
1246 avg=genepred2avg (S);
1248 len=genepred2orf_len (S);
1249 for (a=0; a<S->nseq; a++)
1251 zscore=((double)len[a]-avg)/sd;
1252 zsum2+=FABS(zscore);
1254 zsum2/=(float)S->nseq;
1258 int *genepred2orf_len (Sequence *S)
1261 len=vcalloc (S->nseq, sizeof (int));
1262 for (a=0; a<S->nseq; a++)
1263 for (b=0; b<S->len[a]; b++)
1264 len[a]+=(isupper(S->seq[a][b]));
1268 Alignment *coffee_seq_evaluate_output_old2 ( Alignment *IN, Constraint_list *CL)
1271 int avg, min_avg, best_cycle;
1272 char **pred_list; FILE *fp;
1276 pred_list=declare_char (ncycle, 100);
1279 //CL=expand_constraint_list_4gp(CL);
1280 min_avg=constraint_list2avg(CL);
1282 for (a=1; a<ncycle && CL->ne>0; a++)
1285 if (w);free_int (w, -1);
1286 w=list2residue_total_weight (CL);
1287 CL=relax_constraint_list_4gp (CL);
1288 fprintf (stderr,"\nRELAX CYCLE: %2d AVG: %5d [%10d] ", a, avg=constraint_list2avg(CL), CL->ne);
1289 for (b=0; b<S->nseq; b++)for (c=0; c<S->len[b]; c++)w[b][c]-=1;
1291 //rescale nuclotide coding weights
1293 PS=dnaseq2geneseq (S, w);
1294 vfree (display_accuracy (genepred_seq2accuracy_counts (S, PS, NULL),stderr));
1297 free_sequence (PS, PS->nseq);
1305 S=get_fasta_sequence (pred_list[best_cycle], NULL);
1306 vfree (display_accuracy (genepred_seq2accuracy_counts ((CL->S), S, NULL),stderr));myexit (0);
1307 for (a=0; a<S->nseq; a++)
1308 HERE (">%s\n%s", S->name[a], S->seq[a]);
1310 return seq2aln (S,NULL, 1);
1315 Alignment * non_extended_t_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
1317 int a,b, c,res, s1, s2, r1, r2;
1318 Alignment *OUT=NULL;
1320 int max_score_r, score_r;
1321 double score_col=0, score_aln=0;
1322 int max_score_col, max_score_aln;
1323 double *max_score_seq, *score_seq;
1325 int **tot_non_extended_weight;
1326 int **res_non_extended_weight;
1328 CLIST_TYPE *entry=NULL;
1332 entry=vcalloc (CL->entry_len+1, CL->el_size);
1333 if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
1335 OUT=copy_aln (IN, OUT);
1336 pos=aln2pos_simple(IN, IN->nseq);
1339 max_score_seq=vcalloc ( IN->nseq, sizeof (double));
1340 score_seq=vcalloc ( IN->nseq, sizeof (double));
1342 tot_non_extended_weight=list2residue_total_weight(CL);
1343 res_non_extended_weight=declare_int ((CL->S)->nseq, (CL->S)->max_len+1);
1345 for (a=0; a< IN->len_aln; a++)
1347 for ( b=0; b< IN->nseq-1; b++)
1351 for ( c=b+1; c< IN->nseq; c++)
1355 if ( s1==s2 && !CL->do_self)continue;
1356 else if ( r1<=0 || r2<=0) continue;
1363 if ((l=main_search_in_list_constraint (entry,&p,4,CL))!=NULL)
1365 res_non_extended_weight[s1][r1]+=l[WE];
1366 res_non_extended_weight[s2][r2]+=l[WE];
1372 if ((l=main_search_in_list_constraint (entry,&p,4,CL))!=NULL)
1374 res_non_extended_weight[s1][r1]+=l[WE];
1375 res_non_extended_weight[s2][r2]+=l[WE];
1377 max_score=MAX(max_score,res_non_extended_weight[s1][r1]);
1378 max_score=MAX(max_score,res_non_extended_weight[s2][r2]);
1385 sprintf ( OUT->name[IN->nseq], "cons");
1386 for ( max_score_aln=0,score_aln=0,a=0; a< IN->len_aln; a++)
1388 OUT->seq_al[IN->nseq][a]=NO_COLOR_RESIDUE;
1389 for ( local_nseq=0,b=0; b<IN->nseq; b++){local_nseq+=(pos[b][a]>0)?1:0;}
1391 for ( max_score_col=0, score_col=0,b=0; b< IN->nseq; b++)
1393 OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
1399 max_score_r =max_score;/*tot_non_extended_weight[s1][r1];*/
1400 score_r=res_non_extended_weight[s1][r1];
1401 res=(max_score_r==0 || local_nseq<2 )?NO_COLOR_RESIDUE:((score_r*10)/max_score_r);
1403 (OUT)->seq_al[b][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res, 9));
1404 max_score_col+=max_score_r;
1406 max_score_seq[b]+=max_score_r;
1407 score_seq[b]+=score_r;
1408 max_score_aln+=max_score_r;
1411 res=(max_score_col==0 || local_nseq<2)?NO_COLOR_RESIDUE:((score_col*10)/max_score_col);
1412 OUT->seq_al[IN->nseq][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
1415 IN->score_aln=OUT->score_aln=(max_score_aln==0)?0:((score_aln*100)/max_score_aln);
1416 for ( a=0; a< OUT->nseq; a++)
1418 OUT->score_seq[a]=(max_score_seq[a]==0)?0:((score_seq[a]*100)/max_score_seq[a]);
1419 OUT->score_seq[a]=(OUT->score_seq[a]>100)?100:OUT->score_seq[a];
1421 OUT->score_aln=(OUT->score_aln>100)?100:OUT->score_aln;
1424 vfree ( max_score_seq);
1426 free_int (tot_non_extended_weight, -1);
1427 free_int (res_non_extended_weight, -1);
1435 /*********************************************************************************************/
1437 /* PROFILE/PRofile Functions */
1439 /*********************************************************************************************/
1440 int channel_profile_profile (int *prf1, int *prf2, Constraint_list *CL);
1442 Profile_cost_func get_profile_mode_function (char *name, Profile_cost_func func)
1446 static Profile_cost_func *flist;
1447 static char **nlist;
1451 /*The first time: initialize the list of pairwse functions*/
1452 /*If func==NULL:REturns a pointer to the function associated with a name*/
1453 /*If name is empty:Prints the name of the function associated with name*/
1457 flist=vcalloc ( 100, sizeof (Pwfunc));
1458 nlist=declare_char (100, 100);
1460 flist[nfunc]=cw_profile_profile;
1461 sprintf (nlist[nfunc], "cw_profile_profile");
1464 flist[nfunc]=muscle_profile_profile;
1465 sprintf (nlist[nfunc], "muscle_profile_profile");
1468 flist[nfunc]=channel_profile_profile;
1469 sprintf (nlist[nfunc], "channel_profile_profile");
1473 for ( a=0; a<nfunc; a++)
1475 if ( (name && strm (nlist[a],name)) || flist[a]==func)
1477 if (name)sprintf (name,"%s", nlist[a]);
1481 fprintf ( stderr, "\n[%s] is an unknown profile_profile function[FATAL:%s]\n",name, PROGRAM);
1486 int generic_evaluate_profile_score (Constraint_list *CL,Alignment *Profile1,int s1, int r1, Alignment *Profile2,int s2, int r2, Profile_cost_func prf_prf)
1493 /*Generic profile function*/
1496 dummy=vcalloc (10, sizeof(int));
1497 dummy[0]=1;/*Number of Amino acid types on colum*/
1498 dummy[1]=5;/*Length of Dummy*/
1499 dummy[3]='\0';/*Amino acid*/
1500 dummy[4]=1; /*Number of occurences*/
1501 dummy[5]=100; /*Frequency in the MSA column*/
1510 prf1=(Profile1)?(Profile1->P)->count2[r1]:NULL;
1511 prf2=(Profile2)?(Profile2->P)->count2[r2]:NULL;
1513 if (!prf1) {prf1=dummy; prf1[3]=(CL->S)->seq[s1][r1];}
1514 else if (!prf2){prf2=dummy; prf2[3]=(CL->S)->seq[s2][r2];}
1516 score=((prf_prf==NULL)?cw_profile_profile:prf_prf) (prf1, prf2, CL);
1523 int cw_profile_profile_count (int *prf1, int *prf2, Constraint_list *CL)
1525 /*see function aln2count2 for prf structure*/
1531 for ( n=0,a=3; a<prf1[1]; a+=3)
1532 for ( b=3; b<prf2[1]; b+=3)
1538 score+=prf1[a+1]*prf2[b+1]*CL->M[res1-'A'][res2-'A'];
1539 n+=prf1[a+1]*prf2[b+1];
1543 score=(score*SCORE_K)/n;
1546 int muscle_profile_profile (int *prf1, int *prf2, Constraint_list *CL)
1548 /*see function aln2count2 for prf structure*/
1551 double score=0, fg1, fg2, fi, fj, m;
1552 static double *exp_lu;
1555 exp_lu=vcalloc ( 10000, sizeof (double));
1557 for ( a=-1000; a<1000; a++)
1558 exp_lu[a]=exp((double)a);
1563 for (a=3; a<prf1[1]; a+=3)
1566 fi=(double)prf1[a+2]/100;
1568 for ( b=3; b<prf2[1]; b+=3)
1571 fj=(double)prf2[b+2]/100;
1572 /*m=exp((double)CL->M[res1-'A'][res2-'A']);*/
1573 m=exp_lu[CL->M[res1-'A'][res2-'A']];
1578 fg1=(double)prf1[2]/100;
1579 fg2=(double)prf2[2]/100;
1580 score=(score==0)?0:log(score)*(1-fg1)*(1-fg2);
1581 score=(score*SCORE_K);
1582 /*if ( score<-100)fprintf ( stderr, "\nSCORE %d %d", (int)score, cw_profile_profile(prf1, prf2, CL));*/
1586 int cw_profile_profile (int *prf1, int *prf2, Constraint_list *CL)
1588 /*see function aln2count2 for prf structure*/
1594 for ( n=0,a=3; a<prf1[1]; a+=3)
1595 for ( b=3; b<prf2[1]; b+=3)
1600 p=prf1[a+1]*prf2[b+1];
1603 score+=p*CL->M[res1-'A'][res2-'A'];
1605 score=(score*SCORE_K)/((double)(n==0)?1:n);
1608 int cw_profile_profile_old (int *prf1, int *prf2, Constraint_list *CL)
1610 /*see function aln2count2 for prf structure*/
1617 for ( n=0,a=3; a<prf1[1]; a+=3)
1618 for ( b=3; b<prf2[1]; b+=3)
1623 p=prf1[a+1]*prf2[b+1];
1626 score+=p*CL->M[res1-'A'][res2-'A'];
1628 score=(score*SCORE_K)/((double)(n==0)?1:n);
1631 int channel_profile_profile ( int *prf1, int *prf2, Constraint_list *CL)
1640 if (prf1[0]!=prf1[0]){fprintf ( stderr, "\nERROR: Inconsistent number of channels [channel_profile_profile::FATAL%s]", PROGRAM);}
1644 for (a=1, n=0; a<=prf1[0]; a++)
1646 if (prf1[a]>0 && prf2[a]>0)
1648 n++;score+=CL->M[prf1[a]-'A'][prf2[a]-'A'];
1655 score=(n==0)?0:(score*SCORE_K)/n;
1661 /*********************************************************************************************/
1663 /* FUNCTIONS FOR GETING THE COST : (Sequences) ->evaluate_residue_pair */
1665 /*********************************************************************************************/
1666 int initialize_scoring_scheme (Constraint_list *CL)
1669 else if (!CL->evaluate_residue_pair)return 0;
1670 else if ( !CL->S) return 0;
1671 else if ( (CL->S)->nseq<2) return 0;
1672 else if ( strlen ((CL->S)->seq[0])==0)return 0;
1673 else if ( strlen ((CL->S)->seq[1])==0)return 0;
1676 //CL->evaluate_residue_pair (CL,13,1,16,1);
1677 CL->evaluate_residue_pair (CL,0,1,1,1);
1683 int evaluate_blast_profile_score (Constraint_list *CL, int s1, int r1, int s2, int r2)
1689 PRF1=(Alignment*)atop(seq2T_value (CL->S, s1, "A", "_RB_"));
1690 PRF2=(Alignment*)atop(seq2T_value (CL->S, s2, "A", "_RB_"));
1692 return generic_evaluate_profile_score (CL,PRF1,s1, r1, PRF2,s2, r2, CL->profile_mode);
1695 int evaluate_aln_profile_score (Constraint_list *CL, int s1, int r1, int s2, int r2)
1698 return generic_evaluate_profile_score (CL,seq2R_template_profile((CL->S),s1),s1, r1, seq2R_template_profile(CL->S,s2),s2, r2, CL->profile_mode);
1702 int evaluate_profile_score (Constraint_list *CL,Alignment *Prf1, int s1, int r1, Alignment *Prf2,int s2, int r2)
1705 return generic_evaluate_profile_score (CL, Prf1, s1,r1,Prf2, s2,r2,CL->profile_mode);
1708 int evaluate_cdna_matrix_score (Constraint_list *CL, int s1, int r1, int s2, int r2)
1717 a1=translate_dna_codon((CL->S)->seq[s1]+r1,'x');
1718 a2=translate_dna_codon((CL->S)->seq[s2]+r2,'x');
1722 if (a1=='x' || a2=='x')return 0;
1723 else return CL->M[a1-'A'][a2-'A']*SCORE_K;
1730 int evaluate_physico_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1734 static float **prop_table;
1737 if (r1<0 || r2<0)return 0;
1740 prop_table= initialise_aa_physico_chemical_property_table(&n_prop);
1741 for (p=0; p< n_prop; p++)max+=100;
1744 a=tolower (( CL->S)->seq[s1][r1]);
1745 b=tolower (( CL->S)->seq[s2][r2]);
1747 for (tot=0,p=0; p< n_prop; p++)
1749 tot+=(double)(prop_table[p][a]-prop_table[p][b])*(prop_table[p][a]-prop_table[p][b]);
1752 tot=(sqrt(tot)/max)*10;
1754 return (int) tot*SCORE_K;
1759 int evaluate_diaa_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1773 m=declare_arrayN(4, sizeof (int), 26, 26, 26, 26);
1774 fp=vfopen ("diaa_mat.mat", "r");
1775 while ((c=fgetc (fp))!=EOF)
1779 buf=vfgets(buf, fp);
1784 sscanf (buf, "%s %s %d %d", k1, k2, &v1, &v2);
1786 m[k1[0]-'a'][k1[1]-'a'][k2[0]-'a'][k2[1]-'a']=v1;
1787 m[k2[0]-'a'][k2[1]-'a'][k1[0]-'a'][k1[1]-'a']=v1;
1791 alp=vcalloc (256, sizeof (int));
1792 for (a=0; a<26; a++)alp[a+'a']=1;
1805 char aa1, aa2, aa3, aa4, u;
1812 aa1=tolower((CL->S)->seq[s1][r1-1]);
1813 aa2=tolower((CL->S)->seq[s1][r1]);
1814 aa3=tolower((CL->S)->seq[s2][r2-1]);
1815 aa4=tolower((CL->S)->seq[s2][r2]);
1816 u=alp[(int)aa1];u+=alp[(int)aa2];u+=alp[(int)aa3];u+=alp[(int)aa4];
1819 s+=m[aa1-'a'][aa2-'a'][aa3-'a'][aa4-'a'];
1824 aa1=tolower((CL->S)->seq[s1][r1]);
1825 aa2=tolower((CL->S)->seq[s1][r1+1]);
1826 aa3=tolower((CL->S)->seq[s2][r2]);
1827 aa4=tolower((CL->S)->seq[s2][r2+1]);
1828 u=alp[(int)aa1];u+=alp[(int)aa2];u+=alp[(int)aa3];u+=alp[(int)aa4];
1831 s+=m[aa1-'a'][aa2-'a'][aa3-'a'][aa4-'a'];
1834 if (n)return (s*SCORE_K)/n;
1838 int evaluate_monoaa_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1852 m=declare_arrayN(2, sizeof (int), 26, 26);
1853 fp=vfopen ("monoaa_mat.mat", "r");
1854 while ((c=fgetc (fp))!=EOF)
1858 buf=vfgets(buf, fp);
1863 sscanf (buf, "%s %s %d %d", k1, k2, &v1, &v2);
1865 m[k1[0]-'a'][k2[0]-'a']=v1;
1866 m[k2[0]-'a'][k1[0]-'a']=v1;
1870 alp=vcalloc (256, sizeof (int));
1871 for (a=0; a<26; a++)alp[a+'a']=1;
1891 aa1=tolower((CL->S)->seq[s1][r1]);
1892 aa3=tolower((CL->S)->seq[s2][r2]);
1893 u=alp[(int)aa1];u+=alp[(int)aa3];
1896 s+=m[aa1-'a'][aa3-'a'];
1901 if (n)return (s*SCORE_K)/n;
1905 int evaluate_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1908 if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
1910 return evaluate_aln_profile_score ( CL, s1,r1, s2, r2);
1918 return CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K;
1923 int *get_curvature ( int s1, Constraint_list *CL);
1924 int evaluate_curvature_score( Constraint_list *CL, int s1, int r1, int s2, int r2)
1931 if (!st) st= vcalloc ((CL->S)->nseq, sizeof (char*));
1934 st[s1]=get_curvature (s1, CL);
1938 st[s2]=get_curvature (s2, CL);
1954 //HERE ("%d", score);
1955 //if (score<0)HERE ("%d", score);
1964 int *get_curvature ( int s1, Constraint_list *CL)
1969 char name [1000], b1[100], b2[100];
1971 sprintf ( name, "%s.curvature", (CL->S)->name[s1]);
1972 array=vcalloc (strlen ((CL->S)->seq[s1]), sizeof (int));
1973 fp=vfopen ( name, "r");
1974 while ( fscanf (fp, "%s %d %c %f\n",b1, &a, &c,&f )==4)
1976 if ( c!=(CL->S)->seq[s1][n]){HERE ("ERROR: %c %c", c,(CL->S)->seq[s1][n] );myexit (0);}
1977 else array[n++]=(int)(float)100*(float)f;
1982 int evaluate_tm_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1989 st= vcalloc ((CL->S)->nseq, sizeof (char*));
1990 RF=atoigetenv ("TM_FACTOR_4_TCOFFEE");
1994 if (!st[s1])st[s1]=seq2T_template_string((CL->S),s1);
1995 if (!st[s2])st[s2]=seq2T_template_string((CL->S),s2);
1997 if (r1<=0 || r2<=0)return 0;
1998 else if (st[s1] && st[s2] && (st[s1][r1-1]==st[s2][r2-1]))F=RF;
1999 else if (st[s1] && st[s2] && (st[s1][r1-1]!=st[s2][r2-1]))F=0;
2002 return (int)(evaluate_matrix_score (CL, s1, r1, s2, r2))+F;
2004 int evaluate_ssp_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
2011 st= vcalloc ((CL->S)->nseq, sizeof (char*));
2012 RF=atoigetenv ("SSP_FACTOR_4_TCOFFEE");
2016 if (!st[s1])st[s1]=seq2T_template_string((CL->S),s1);
2017 if (!st[s2])st[s2]=seq2T_template_string((CL->S),s2);
2019 if (r1<=0 || r2<=0)return 0;
2020 else if (st[s1] && st[s2] && (st[s1][r1-1]==st[s2][r2-1]))F=RF;
2021 else if (st[s1] && st[s2] && (st[s1][r1-1]!=st[s2][r2-1]))F=0;
2024 return (int)(evaluate_matrix_score (CL, s1, r1, s2, r2))+F;
2028 int evaluate_combined_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
2031 function documentation: start
2033 int evaluate_matrix_score ( Constraint_list *CL, int s1, int s2, int r1, int r2)
2035 this function evaluates the score for matching residue S1(r1) wit residue S2(r2)
2038 function documentation: end
2041 if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
2043 return evaluate_aln_profile_score ( CL, s1,r1, s2, r2);
2050 if (r1==0 || r2==0)return CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K;
2053 int A1, A2, B1, B2, a2, b2;
2056 A1=toupper((CL->S)->seq[s1][r1-1]);
2057 A2=toupper((CL->S)->seq[s1][r1]);
2058 B1=toupper((CL->S)->seq[s2][r2-1]);
2059 B2=toupper((CL->S)->seq[s2][r2]);
2063 A1-='A';A2-='A';B1-='A'; B2-='A';a2-='A';b2-='A';
2065 score=CL->M[a2][b2]-FABS((CL->M[A1][A2])-(CL->M[B1][B2]));
2074 int residue_pair_non_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2078 This is the generic Function->works with everything
2080 int residue_pair_non_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2, int field );
2082 Computes the non extended score for aligning residue seq1(r1) Vs seq2(r2)
2084 This function can compare a sequence with itself.
2086 Associated functions: See util constraint list, list extention functions.
2087 function documentation: end
2097 field = CL->weight_field;
2100 if ( r1<=0 || r2<=0)return 0;
2101 else if ( !CL->extend_jit)
2103 if ( !entry) entry=vcalloc (LIST_N_FIELDS , sizeof (int));
2108 if ( r1==r2 && s1==s2) return UNDEFINED;
2109 r=main_search_in_list_constraint( entry,&p,4,CL);
2110 if (r==NULL)return 0;
2111 else return r[field]*SCORE_K;
2114 return UNDEFINED;/*ERROR*/
2121 int residue_pair_extended_list_mixt (Constraint_list *CL, int s1, int r1, int s2, int r2 )
2125 score+= residue_pair_extended_list_quadruplet(CL, s1, r1, s2, r2);
2126 score+= residue_pair_extended_list (CL, s1, r1, s2, r2);
2128 return score*SCORE_K;
2131 int residue_pair_extended_list_quadruplet (Constraint_list *CL, int s1, int r1, int s2, int r2 )
2135 int t_s, t_r, t_w, q_s, q_r, q_w;
2140 /* This measure the quadruplets cost on a pair of residues*/
2144 field=CL->weight_field;
2146 if ( r1<=0 || r2<=0)return 0;
2149 hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
2150 for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
2154 hasch[s1][r1]=100000;
2155 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2157 t_s=CL->residue_index[s1][r1][a+SEQ2];
2158 t_r=CL->residue_index[s1][r1][a+R2];
2159 t_w=CL->residue_index[s1][r1][a+WE];
2160 if ( CL->seq_for_quadruplet[t_s])
2162 for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2164 q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2165 q_r=CL->residue_index[t_s][t_r][b+R2];
2166 q_w=CL->residue_index[t_s][t_r][b+WE];
2167 if (CL-> seq_for_quadruplet[q_s])
2168 hasch[q_s][q_r]=MIN(q_w,t_w);
2175 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2177 t_s=CL->residue_index[s2][r2][a+SEQ2];
2178 t_r=CL->residue_index[s2][r2][a+R2];
2179 t_w=CL->residue_index[s2][r2][a+WE];
2180 if ( CL->seq_for_quadruplet[t_s])
2182 for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2184 q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2185 q_r=CL->residue_index[t_s][t_r][b+R2];
2186 q_w=CL->residue_index[t_s][t_r][b+WE];
2187 if (hasch[q_s][q_r] && CL->seq_for_quadruplet[q_s])
2188 score+=MIN(hasch[q_s][q_r],MIN(q_w,t_w));
2193 score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
2195 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2197 t_s=CL->residue_index[s1][r1][a+SEQ2];
2198 t_r=CL->residue_index[s1][r1][a+R2];
2199 t_w=CL->residue_index[s1][r1][a+WE];
2200 if ( CL->seq_for_quadruplet[t_s])
2202 for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2204 q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2205 q_r=CL->residue_index[t_s][t_r][b+R2];
2211 return (int)(score*SCORE_K);
2215 Constraint_list * R_extension ( Constraint_list *CL, Constraint_list *R);
2216 int residue_pair_extended_list4rna4 (Constraint_list *CL,int s1, int r1, int s2, int r2 )
2222 sprintf ( CL->rna_lib, "%s", seq2rna_lib (CL->S, NULL));
2225 return residue_pair_extended_list4rna2 (CL, s1, r1, s2,r2);
2227 int residue_pair_extended_list4rna1 (Constraint_list *CL, int s1, int r1, int s2, int r2 )
2229 static Constraint_list *R;
2230 if (!R)R=read_rna_lib (CL->S, CL->rna_lib);
2231 return residue_pair_extended_list4rna (CL, R, s1, r1, s2, r2);
2234 int residue_pair_extended_list4rna3 (Constraint_list *CL,int s1, int r1, int s2, int r2 )
2236 static Constraint_list *R;
2239 R=read_rna_lib (CL->S, CL->rna_lib);
2240 rna_lib_extension (CL,R);
2242 return residue_pair_extended_list (CL, s1,r1, s2,r2);
2245 int residue_pair_extended_list4rna2 (Constraint_list *CL,int s1, int r1, int s2, int r2 )
2247 static Constraint_list *R;
2253 R=read_rna_lib (CL->S, CL->rna_lib);
2254 rna_lib_extension (CL,R);
2258 return residue_pair_extended_list4rna (CL, R, s1, r1, s2, r2);
2260 int residue_pair_extended_list4rna ( Constraint_list *CL,Constraint_list *R, int s1, int r1, int s2, int r2 )
2266 int score=0, score2;
2270 if ( r1<0 || r2<0)return 0;
2274 for (a=1; a<R->residue_index[s1][r1][0]; a+=ICHUNK)
2276 list1[n1++]=R->residue_index[s1][r1][a+R2];
2281 for (a=1; a<R->residue_index[s2][r2][0]; a+=ICHUNK)
2283 list2[n2++]=R->residue_index[s2][r2][a+R2];
2287 score=residue_pair_extended_list ( CL, s1,list1[0], s2,list2[0]);
2289 for (score2=0,a=1; a<n1; a++)
2290 for ( b=1; b<n2; b++)
2291 score2=MAX(residue_pair_extended_list ( CL, s1,list1[a], s2,list2[b]), score2);
2293 if ( n1==1 || n2==1);
2294 else score=MAX(score,score2);
2300 int residue_pair_extended_list4rna_ref ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2302 static Constraint_list *R;
2306 int score=0, score2;
2312 list=read_lib_list ( CL->rna_lib, &n);
2314 R=declare_constraint_list ( CL->S,NULL, NULL, 0,NULL, NULL);
2316 for (a=0; a< n; a++)
2318 R=read_constraint_list_file (R, list[a]);
2323 if ( r1<0 || r2<0)return 0;
2327 for (a=1; a<R->residue_index[s1][r1][0]; a+=ICHUNK)
2329 list1[n1++]=R->residue_index[s1][r1][a+R2];
2334 for (a=1; a<R->residue_index[s2][r2][0]; a+=ICHUNK)
2336 list2[n2++]=R->residue_index[s2][r2][a+R2];
2340 score=residue_pair_extended_list ( CL, s1,list1[0], s2,list2[0]);
2342 for (score2=0,a=1; a<n1; a++)
2343 for ( b=1; b<n2; b++)
2344 score2=MAX(residue_pair_extended_list ( CL, s1,list1[a], s2,list2[b]), score2);
2346 if ( n1==1 || n2==1);
2347 else score=MAX(score,score2);
2352 static int ** clean_residue_pair_hasch (int s1, int r1, int s2, int r2,int **hasch, Constraint_list *CL);
2353 int residue_pair_extended_list_raw ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2366 function documentation: start
2368 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2370 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2371 Computes: matrix_score
2375 The extended score depends on the function index_res_constraint_list.
2376 This function can compare a sequence with itself.
2378 Associated functions: See util constraint list, list extention functions.
2380 function documentation: end
2383 field=CL->weight_field;
2386 if ( r1<=0 || r2<=0)return 0;
2387 if ( !hasch || max_len!=(CL->S)->max_len)
2389 max_len=(CL->S)->max_len;
2390 if ( hasch) free_int ( hasch, -1);
2391 hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2396 /* Check matches for R1 in the indexed lib*/
2397 hasch[s1][r1]=FORBIDEN;
2398 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2400 t_s=CL->residue_index[s1][r1][a+SEQ2];
2401 t_r=CL->residue_index[s1][r1][a+R2];
2402 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2405 /*Check Matches for r1 <-> r2 in the indexed lib */
2406 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2408 t_s=CL->residue_index[s2][r2][a+SEQ2];
2409 t_r=CL->residue_index[s2][r2][a+R2];
2412 if (hasch[t_s][t_r])
2414 if (hasch[t_s][t_r]==FORBIDEN)
2416 score+=CL->residue_index[s2][r2][a+field];
2420 delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);
2426 clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2429 int residue_pair_extended_list_4gp ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2443 function documentation: start
2445 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2447 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2448 Computes: matrix_score
2452 The extended score depends on the function index_res_constraint_list.
2453 This function can compare a sequence with itself.
2455 Associated functions: See util constraint list, list extention functions.
2457 function documentation: end
2460 field=CL->weight_field;
2463 if ( r1<=0 || r2<=0)return 0;
2464 if ( !hasch || max_len!=(CL->S)->max_len)
2466 max_len=(CL->S)->max_len;
2467 if ( hasch) free_int ( hasch, -1);
2468 hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2473 /* Check matches for R1 in the indexed lib*/
2474 hasch[s1][r1]=FORBIDEN;
2475 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2477 t_s=CL->residue_index[s1][r1][a+SEQ2];
2478 t_r=CL->residue_index[s1][r1][a+R2];
2479 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2482 /*Check Matches for r1 <-> r2 in the indexed lib */
2483 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2485 t_s=CL->residue_index[s2][r2][a+SEQ2];
2486 t_r=CL->residue_index[s2][r2][a+R2];
2489 if (hasch[t_s][t_r])
2491 if (hasch[t_s][t_r]==FORBIDEN)
2493 score+=((float)CL->residue_index[s2][r2][a+2]/NORM_F);
2497 delta=MAX((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
2503 clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2504 score/=(CL->S)->nseq;
2510 int residue_pair_extended_list_pc ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2512 //old norm does not take into account the number of effective intermediate sequences
2526 function documentation: start
2528 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2530 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2531 Computes: matrix_score
2535 The extended score depends on the function index_res_constraint_list.
2536 This function can compare a sequence with itself.
2538 Associated functions: See util constraint list, list extention functions.
2540 function documentation: end
2543 field=CL->weight_field;
2546 if ( r1<=0 || r2<=0)return 0;
2547 if ( !hasch || max_len!=(CL->S)->max_len)
2549 max_len=(CL->S)->max_len;
2550 if ( hasch) free_int ( hasch, -1);
2551 hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2556 /* Check matches for R1 in the indexed lib*/
2557 hasch[s1][r1]=FORBIDEN;
2558 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2560 t_s=CL->residue_index[s1][r1][a+SEQ2];
2561 t_r=CL->residue_index[s1][r1][a+R2];
2562 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2568 /*Check Matches for r1 <-> r2 in the indexed lib */
2569 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2571 t_s=CL->residue_index[s2][r2][a+SEQ2];
2572 t_r=CL->residue_index[s2][r2][a+R2];
2575 if (hasch[t_s][t_r])
2577 if (hasch[t_s][t_r]==FORBIDEN)
2579 score+=((float)CL->residue_index[s2][r2][a+field]/NORM_F);
2583 delta=MIN((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
2589 clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2592 norm1=MIN(norm1,norm2);
2594 //Old Normailzation: on the number of sequences, useless when not doiang an all against all
2595 //norm1=(CL->S)->nseq;
2597 score=(norm1)?score/norm1:0;
2600 return score*NORM_F;
2602 int residue_pair_extended_list_pc_old ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2604 //old norm does not take into account the number of effective intermediate sequences
2616 function documentation: start
2618 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2620 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2621 Computes: matrix_score
2625 The extended score depends on the function index_res_constraint_list.
2626 This function can compare a sequence with itself.
2628 Associated functions: See util constraint list, list extention functions.
2630 function documentation: end
2633 field=CL->weight_field;
2636 if ( r1<=0 || r2<=0)return 0;
2637 if ( !hasch || max_len!=(CL->S)->max_len)
2639 max_len=(CL->S)->max_len;
2640 if ( hasch) free_int ( hasch, -1);
2641 hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2645 /* Check matches for R1 in the indexed lib*/
2646 hasch[s1][r1]=FORBIDEN;
2647 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2649 t_s=CL->residue_index[s1][r1][a+SEQ2];
2650 t_r=CL->residue_index[s1][r1][a+R2];
2651 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2655 /*Check Matches for r1 <-> r2 in the indexed lib */
2656 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2658 t_s=CL->residue_index[s2][r2][a+SEQ2];
2659 t_r=CL->residue_index[s2][r2][a+R2];
2662 if (hasch[t_s][t_r])
2664 if (hasch[t_s][t_r]==FORBIDEN)
2666 score+=((float)CL->residue_index[s2][r2][a+2]/NORM_F);
2670 delta=MIN((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
2676 clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2677 score/=(CL->S)->nseq;
2678 return score*NORM_F;
2682 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2694 new function: self normalized
2695 function documentation: start
2697 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2699 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2700 Computes: matrix_score
2704 The extended score depends on the function index_res_constraint_list.
2705 This function can compare a sequence with itself.
2707 Associated functions: See util constraint list, list extention functions.
2709 function documentation: end
2712 field=CL->weight_field;
2715 if ( r1<=0 || r2<=0)return 0;
2716 if ( !hasch || max_len!=(CL->S)->max_len)
2718 max_len=(CL->S)->max_len;
2719 if ( hasch) free_int ( hasch, -1);
2720 hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2725 /* Check matches for R1 in the indexed lib*/
2726 hasch[s1][r1]=FORBIDEN;
2727 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2729 t_s=CL->residue_index[s1][r1][a+SEQ2];
2730 t_r=CL->residue_index[s1][r1][a+R2];
2731 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2732 max_score+=CL->residue_index[s1][r1][a+field];
2735 /*Check Matches for r1 <-> r2 in the indexed lib */
2736 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2738 t_s=CL->residue_index[s2][r2][a+SEQ2];
2739 t_r=CL->residue_index[s2][r2][a+R2];
2742 if (hasch[t_s][t_r])
2744 if (hasch[t_s][t_r]==FORBIDEN)
2746 score+=CL->residue_index[s2][r2][a+field];
2747 max_score+=CL->residue_index[s2][r2][a+field];
2752 delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);
2755 max_score-=hasch[t_s][t_r];
2757 max_val=MAX(max_val,delta);
2762 max_score+=CL->residue_index[s2][r2][a+field];
2766 max_score-=hasch[s2][r2];
2767 clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2770 if ( max_score==0)score=0;
2771 else if ( CL->normalise)
2773 score=((score*CL->normalise)/max_score)*SCORE_K;
2774 if (max_val> CL->normalise)
2776 score*=max_val/(double)CL->normalise;
2781 int ** clean_residue_pair_hasch (int s1, int r1, int s2, int r2,int **hasch, Constraint_list *CL)
2784 if ( !hasch) return hasch;
2786 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2788 t_s=CL->residue_index[s1][r1][a+SEQ2];
2789 t_r=CL->residue_index[s1][r1][a+R2];
2792 hasch[s1][r1]=hasch[s2][r2]=0;
2796 int residue_pair_test_function ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2808 function documentation: start
2810 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2812 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2813 Computes: matrix_score
2817 The extended score depends on the function index_res_constraint_list.
2818 This function can compare a sequence with itself.
2820 Associated functions: See util constraint list, list extention functions.
2822 function documentation: end
2826 CL->weight_field=WE;
2827 field=CL->weight_field;
2830 if ( r1<=0 || r2<=0)return 0;
2831 if ( !hasch || max_len!=(CL->S)->max_len)
2833 max_len=(CL->S)->max_len;
2834 if ( hasch) free_int ( hasch, -1);
2835 hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2841 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2843 t_s=CL->residue_index[s1][r1][a+SEQ2];
2844 t_r=CL->residue_index[s1][r1][a+R2];
2845 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2849 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2851 t_s=CL->residue_index[s2][r2][a+SEQ2];
2852 t_r=CL->residue_index[s2][r2][a+R2];
2853 if (hasch[t_s][t_r])
2855 cons1=hasch[t_s][t_r];
2856 cons2=CL->residue_index[s2][r2][a+field];
2857 score +=MIN(cons1,cons2);
2862 score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
2864 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2866 t_s=CL->residue_index[s1][r1][a+SEQ2];
2867 t_r=CL->residue_index[s1][r1][a+R2];
2870 hasch[s1][r1]=hasch[s2][r2]=0;
2873 return (int)(score*SCORE_K);
2876 int residue_pair_relative_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2885 function documentation: start
2887 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2889 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2890 Computes: matrix_score
2894 The extended score depends on the function index_res_constraint_list.
2895 This function can compare a sequence with itself.
2897 Associated functions: See util constraint list, list extention functions.
2899 function documentation: end
2904 field=CL->weight_field;
2907 if ( r1<=0 || r2<=0)return 0;
2908 if ( !hasch || max_len!=(CL->S)->max_len)
2910 max_len=(CL->S)->max_len;
2911 if ( hasch) free_int ( hasch, -1);
2912 hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2917 hasch[s1][r1]=100000;
2918 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2920 t_s=CL->residue_index[s1][r1][a+SEQ2];
2921 t_r=CL->residue_index[s1][r1][a+R2];
2922 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2923 total_score+=CL->residue_index[s1][r1][a+field];
2927 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2929 t_s=CL->residue_index[s2][r2][a+SEQ2];
2930 t_r=CL->residue_index[s2][r2][a+R2];
2931 total_score+=CL->residue_index[s1][r1][a+field];
2932 if (hasch[t_s][t_r])
2934 if (field==WE){score+=2*MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);}
2938 score=((CL->normalise*score)/total_score);
2941 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2943 t_s=CL->residue_index[s1][r1][a+SEQ2];
2944 t_r=CL->residue_index[s1][r1][a+R2];
2947 hasch[s1][r1]=hasch[s2][r2]=0;
2949 return score*SCORE_K;
2951 int residue_pair_extended_list_g_coffee_quadruplet ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2953 int t_s, t_r, t_w, q_s, q_r, q_w;
2959 /* This measure the quadruplets cost on a pair of residues*/
2961 field=CL->weight_field;
2963 if ( r1<=0 || r2<=0)return 0;
2966 hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
2967 for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
2971 hasch[s1][r1]=100000;
2972 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2974 t_s=CL->residue_index[s1][r1][a+SEQ2];
2975 t_r=CL->residue_index[s1][r1][a+R2];
2976 t_w=CL->residue_index[s1][r1][a+field];
2977 if ( CL->seq_for_quadruplet[t_s])
2979 for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2981 q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2982 q_r=CL->residue_index[t_s][t_r][b+R2];
2983 if (CL-> seq_for_quadruplet[q_s])
2986 hasch[q_s][q_r]=MIN(CL->residue_index[t_s][t_r][b+2],t_w);
2993 for (s=0,score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2995 t_s=CL->residue_index[s2][r2][a+SEQ2];
2996 t_r=CL->residue_index[s2][r2][a+R2];
2997 t_w=CL->residue_index[s2][r2][a+field];
2998 if ( CL->seq_for_quadruplet[t_s])
3000 for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
3002 q_s=CL->residue_index[t_s][t_r][b+SEQ2];
3003 q_r=CL->residue_index[t_s][t_r][b+R2];
3004 q_w=CL->residue_index[t_s][t_r][b+field];
3005 if (hasch[q_s][q_r] && CL->seq_for_quadruplet[q_s])
3006 s=MIN(hasch[q_s][q_r],MIN(CL->residue_index[t_s][t_r][b+2],q_w));
3007 score=MAX(score, s);
3012 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3014 t_s=CL->residue_index[s1][r1][a+SEQ2];
3015 t_r=CL->residue_index[s1][r1][a+R2];
3016 t_w=CL->residue_index[s1][r1][a+field];
3017 if ( CL->seq_for_quadruplet[t_s])
3019 for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
3021 q_s=CL->residue_index[t_s][t_r][b+SEQ2];
3022 q_r=CL->residue_index[t_s][t_r][b+R2];
3028 return score*SCORE_K;
3030 int residue_pair_extended_list_g_coffee ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
3038 function documentation: start
3040 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
3042 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
3043 Computes: matrix_score
3047 The extended score depends on the function index_res_constraint_list.
3048 This function can compare a sequence with itself.
3050 Associated functions: See util constraint list, list extention functions.
3052 function documentation: end
3055 field=CL->weight_field;
3057 if ( r1<=0 || r2<=0)return 0;
3060 hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
3061 for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
3066 hasch[s1][r1]=100000;
3067 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3069 t_s=CL->residue_index[s1][r1][a+SEQ2];
3070 t_r=CL->residue_index[s1][r1][a+R2];
3071 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
3075 for (s=0, score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
3077 t_s=CL->residue_index[s2][r2][a+SEQ2];
3078 t_r=CL->residue_index[s2][r2][a+R2];
3080 if (hasch[t_s][t_r])
3083 {s=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+WE]);
3089 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3091 t_s=CL->residue_index[s1][r1][a+SEQ2];
3092 t_r=CL->residue_index[s1][r1][a+R2];
3095 hasch[s1][r1]=hasch[s2][r2]=0;
3097 return score*SCORE_K;
3100 int extend_residue_pair ( Constraint_list *CL, int s1, int r1, int s2, int r2)
3114 This is the generic Function->works with everything
3115 should be gradually phased out
3118 int extend_residue_pair ( Constraint_list *CL, int s1, int r1, int s2, int r2, int field )
3120 Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
3121 Computes: matrix_score
3125 The extended score depends on the function index_res_constraint_list.
3126 This function can compare a sequence with itself.
3128 Associated functions: See util constraint list, list extention functions.
3129 function documentation: end
3133 field=CL->weight_field;
3135 if ( r1<=0 || r2<=0)return 0;
3136 else if ( !CL->residue_index && CL->M)
3138 return evaluate_matrix_score (CL, s1,r1, s2, r2);
3141 else if ( !CL->extend_jit)
3143 if ( !entry) entry=vcalloc (LIST_N_FIELDS , sizeof (int));
3148 r=main_search_in_list_constraint( entry,&p,4,CL);
3149 if (r==NULL)return 0;
3150 else return r[field];
3156 hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
3157 for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
3162 hasch[s1][r1]=100000;
3163 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3165 t_s=CL->residue_index[s1][r1][a+SEQ2];
3166 t_r=CL->residue_index[s1][r1][a+R2];
3167 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+WE];
3171 for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
3173 t_s=CL->residue_index[s2][r2][a+SEQ2];
3174 t_r=CL->residue_index[s2][r2][a+R2];
3175 if (hasch[t_s][t_r])
3177 if (field==WE)score+=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+WE] );
3181 score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
3182 for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3184 t_s=CL->residue_index[s1][r1][a+SEQ2];
3185 t_r=CL->residue_index[s1][r1][a+R2];
3188 hasch[s1][r1]=hasch[s2][r2]=0;
3190 return (int)(score*SCORE_K);
3193 /*********************************************************************************************/
3195 /* FUNCTIONS FOR GETTING THE PW COST : CL->get_dp_cost */
3197 /*********************************************************************************************/
3198 int get_dp_cost_blosum_matrix (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3201 static int **matrix;
3203 if (!matrix) matrix=read_matrice ("blosum62mt");
3204 s1=A->order[list1[0]][0];
3205 s2=A->order[list2[0]][0];
3206 r1=pos1[list1[0]][col1];
3207 r2=pos2[list2[0]][col2];
3209 /*dp cost function: works only with two sequences*/
3211 if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
3212 return evaluate_aln_profile_score ( CL, s1,r1, s2, r2) -CL->nomatch*SCORE_K;
3213 else if (r1>0 && r2>0)
3219 return matrix [(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K -CL->nomatch*SCORE_K;
3223 return -CL->nomatch*SCORE_K ;
3225 int get_dp_cost_pam_matrix (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3228 static int **matrix;
3230 if (!matrix) matrix=read_matrice ("pam250mt");
3231 s1=A->order[list1[0]][0];
3232 s2=A->order[list2[0]][0];
3233 r1=pos1[list1[0]][col1];
3234 r2=pos2[list2[0]][col2];
3236 /*dp cost function: works only with two sequences*/
3239 if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
3240 return evaluate_aln_profile_score ( CL, s1,r1, s2, r2) -CL->nomatch*SCORE_K;
3241 else if (r1>0 && r2>0)
3247 return matrix [(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K -CL->nomatch*SCORE_K;
3251 return -CL->nomatch*SCORE_K ;
3254 int get_dp_cost_pw_matrix (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3258 s1=A->order[list1[0]][0];
3259 s2=A->order[list2[0]][0];
3260 r1=pos1[list1[0]][col1];
3261 r2=pos2[list2[0]][col2];
3263 /*dp cost function: works only with two sequences*/
3264 if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
3265 return evaluate_aln_profile_score ( CL, s1,r1, s2, r2) -CL->nomatch*SCORE_K;
3266 else if (r1>0 && r2>0)
3272 return CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K -CL->nomatch*SCORE_K;
3276 return -CL->nomatch*SCORE_K ;
3279 /*********************************************************************************************/
3281 /* FUNCTIONS FOR GETTING THE COST : CL->get_dp_cost */
3283 /*********************************************************************************************/
3287 int get_cdna_best_frame_dp_cost (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3294 static Alignment *B;
3297 if ( !score)score=declare_int(3, 2);
3307 if (ns1+ns2>2){fprintf ( stderr, "\nERROR: get_cdna_dp_cost mode is only for pair-wise ALN [FATAL]\n");crash("");}
3309 B=copy_aln (A, NULL);
3311 l1=(int)strlen ( A->seq_al[list1[0]]);
3312 for ( b=0; b<l1; b++)
3313 B->seq_al[list1[0]][b]=translate_dna_codon (A->seq_al[list1[0]]+b, 'x');
3314 l2=(int)strlen ( A->seq_al[list2[0]]);
3315 for ( b=0; b<l2; b++)
3316 B->seq_al[list2[0]][b]=translate_dna_codon (A->seq_al[list2[0]]+b, 'x');
3321 for ( a=0; a< 3; a++)score[a][0]=score[a][1]=0;
3322 for ( a=col1-(n*3),b=col2-(n*3); a<col1+(n*3) ; a++, b++)
3324 if ( a<0 || b<0 || a>=l1 || b>=l2)continue;
3326 a1=tolower(B->seq_al[list1[0]][a]);
3327 a2=tolower(B->seq_al[list2[0]][b]);
3329 score[a%3][0]+=(a1=='x' || a2=='x')?0:CL->M[a1-'A'][a2-'A'];
3333 for ( a=0; a< 3; a++)score[a][0]=(score[a][1]>0)?(score[a][0]/score[a][1]):0;
3334 if ( score[0][0]>score[1][0] && score[0][0]>score[2][0])
3336 else if ( score[1][0]>score[0][0] && score[1][0]>score[2][0])
3344 int get_dp_cost_quadruplet ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3349 if ( ns1==1 || ns2==1)
3350 score=slow_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3352 score=fast_get_dp_cost_quadruplet ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3354 return (score==UNDEFINED)?UNDEFINED:(score-SCORE_K*CL->nomatch);
3356 int get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3363 if (A==NULL)return 0;
3365 if (MODE!=2 || MODE==0 || (!CL->residue_index && CL->M) || (!CL->residue_index && CL->T)|| ns1==1 || ns2==1)
3366 score=slow_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3367 else if (MODE==1 || MODE==2)
3368 score=fast_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3374 return (score==UNDEFINED)?UNDEFINED:(score-SCORE_K*CL->nomatch);
3376 int ***make_cw_lu (int **cons, int l, Constraint_list *CL);
3377 int ***make_cw_lu (int **cons, int l, Constraint_list *CL)
3382 lu=declare_arrayN(3, sizeof (int),l,26, 2);
3383 for ( p=0; p<l ; p++)
3385 for (r=0; r<26; r++)
3387 for (a=3; a<cons[p][1]; a+=3)
3389 lu[p][r][0]+=cons[p][a+1]*CL->M[r][cons[p][a]-'A'];
3390 lu[p][r][1]+=cons[p][a+1];
3397 int cw_profile_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3399 static int last_tag;
3400 static int *pr, ***lu;
3402 static int *list[2], ns[2], **cons[2], ref;
3403 int eva_col,ref_col, a, p, r;
3409 if (last_tag!=A->random_tag)
3413 last_tag=A->random_tag;
3414 list[0]=list1;ns[0]=ns1;
3415 list[1]=list2;ns[1]=ns2;
3416 free_int (cons[0],-1);free_int (cons[1],-1);free_arrayN((void*)lu,3);
3417 cons[0]=NULL;cons[1]=NULL;lu=NULL;
3419 n1=sub_aln2nseq_prf (A, ns[0], list[0]);
3420 n2=sub_aln2nseq_prf (A, ns[1], list[1]);
3423 cons[0]=sub_aln2count_mat2 (A, ns[0], list[0]);
3424 cons[1]=sub_aln2count_mat2 (A, ns[1], list[1]);
3425 ref=(ns[0]>ns[1])?0:1;
3426 lu=make_cw_lu(cons[ref],(int)strlen(A->seq_al[list[ref][0]]), CL);
3433 r1=A->seq_al[list1[0]][col1];
3434 r2=A->seq_al[list2[0]][col2];
3435 if ( r1!='-' && r2!='-')
3436 return CL->M[r1-'A'][r2-'A']*SCORE_K -SCORE_K*CL->nomatch;
3438 return -SCORE_K*CL->nomatch;
3442 eva_col= (ref==0)?col2:col1;
3443 ref_col= (ref==0)?col1:col2;
3444 pr=cons[1-ref][eva_col];
3446 for (a=3; a< pr[1]; a+=3)
3451 t1+=lu[ref_col][r-'a'][0]*p;
3452 t2+=lu[ref_col][r-'a'][1]*p;
3454 score=(t2==0)?0:(t1*SCORE_K)/t2;
3455 score -=SCORE_K*CL->nomatch;
3459 int cw_profile_get_dp_cost_old ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3461 static int last_tag;
3462 static int **cons1, **cons2;
3466 if (last_tag!=A->random_tag)
3468 last_tag=A->random_tag;
3469 free_int (cons1,-1);free_int (cons2,-1);
3470 cons1=sub_aln2count_mat2 (A, ns1, list1);
3471 cons2=sub_aln2count_mat2 (A, ns2, list2);
3473 score=cw_profile_profile (cons1[col1], cons2[col2], CL)-SCORE_K*CL->nomatch;
3477 int cw_profile_get_dp_cost_window ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3479 static int last_tag;
3480 static int **cons1, **cons2;
3481 int a, score, window_size=5, n, p1, p2;
3484 if (last_tag!=A->random_tag)
3486 last_tag=A->random_tag;
3487 free_int (cons1,-1);free_int (cons2,-1);
3488 cons1=sub_aln2count_mat2 (A, ns1, list1);
3489 cons2=sub_aln2count_mat2 (A, ns2, list2);
3492 for (n=0,score=0,a=0; a<window_size; a++)
3496 if ( p1<0 || cons1[p1][0]==END_ARRAY)continue;
3497 if ( p2<0 || cons2[p2][0]==END_ARRAY)continue;
3498 score+=cw_profile_profile (cons1[col1], cons2[col2], CL)-SCORE_K*CL->nomatch;
3506 int consensus_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3508 static int last_tag;
3509 static char *seq1, *seq2;
3512 /*Works only with matrix*/
3513 if (last_tag !=A->random_tag)
3515 int ls1, ls2, lenls1, lenls2;
3517 last_tag=A->random_tag;
3518 vfree (seq1);vfree (seq2);
3519 seq1=sub_aln2cons_seq_mat (A, ns1, list1, "blosum62mt");
3520 seq2=sub_aln2cons_seq_mat (A, ns2, list2, "blosum62mt");
3521 ls1=list1[ns1-1];ls2=list2[ns2-1];
3522 lenls1=(int)strlen (A->seq_al[ls1]); lenls2=(int)strlen (A->seq_al[ls2]);
3525 return (CL->M[seq1[col1]-'A'][seq2[col2]-'A']*SCORE_K)-SCORE_K*CL->nomatch;
3528 int fast_get_dp_cost_quadruplet ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3530 /*WARNING: WORKS ONLY WITH List to Extend*/
3531 /*That function does a quruple extension beween two columns by pooling the residues together*/
3539 int s1, rs1, r1, t_r, t_s,t_w, q_r, q_s, q_w, s2, rs2, r2;
3540 int **buf_pos, buf_ns, *buf_list, buf_col;
3542 static int **hasch1;
3543 static int **hasch2;
3545 static int **n_hasch1;
3546 static int **n_hasch2;
3548 static int **is_in_col1;
3549 static int **is_in_col2;
3574 hasch1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3575 hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3576 n_hasch1=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
3577 n_hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3578 is_in_col1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3579 is_in_col2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3582 for ( a=0; a< ns1; a++)
3585 s1=A->order[rs1][0];
3591 is_in_col1[s1][r1]=1;
3592 for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3594 t_s=CL->residue_index[s1][r1][b+SEQ2];
3595 t_r=CL->residue_index[s1][r1][b+R2];
3596 t_w=CL->residue_index[s1][r1][b+WE];
3597 for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3599 q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3600 q_r=CL->residue_index[t_s][t_r][c+R2];
3601 q_w=CL->residue_index[t_s][t_r][c+WE];
3602 hasch1[q_s][q_r]+=MIN(q_w, t_w);
3603 n_hasch1[q_s][q_r]++;
3609 for ( a=0; a< ns2; a++)
3612 s2=A->order[rs2][0];
3618 is_in_col2[s2][r2]=1;
3619 for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3621 t_s=CL->residue_index[s2][r2][b+SEQ2];
3622 t_r=CL->residue_index[s2][r2][b+R2];
3623 t_w=CL->residue_index[s2][r2][b+WE];
3624 for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3626 q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3627 q_r=CL->residue_index[t_s][t_r][c+R2];
3628 q_w=CL->residue_index[t_s][t_r][c+WE];
3629 hasch2[q_s][q_r]+=MIN(t_w, q_w);
3630 n_hasch2[q_s][q_r]++;
3637 for ( a=0; a< ns2; a++)
3640 s2=A->order[rs2][0];
3646 for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3648 t_s=CL->residue_index[s2][r2][b+SEQ2];
3649 t_r=CL->residue_index[s2][r2][b+R2];
3651 for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3653 q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3654 q_r=CL->residue_index[t_s][t_r][c+R2];
3655 if ( hasch2[q_s][q_r] && hasch1[q_s][q_r]&& !(is_in_col1[q_s][q_r] || is_in_col2[q_s][q_r]))
3657 score+=MIN(hasch2[q_s][q_r]*(n_hasch1[q_s][q_r]),hasch1[q_s][q_r]*(n_hasch2[q_s][q_r]));
3659 else if ( hasch2[q_s][q_r] && is_in_col1[q_s][q_r])
3661 score+=hasch2[q_s][q_r]*(n_hasch1[q_s][q_r]+1);
3663 else if (hasch1[q_s][q_r] && is_in_col2[q_s][q_r])
3665 score+=hasch1[q_s][q_r]*(n_hasch2[q_s][q_r]+1);
3668 n_hasch2[q_s][q_r]=0;
3672 is_in_col2[s2][r2]=0;
3677 for ( a=0; a< ns1; a++)
3680 s1=A->order[rs1][0];
3686 is_in_col1[s1][r1]=0;
3688 for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3690 t_s=CL->residue_index[s1][r1][b+SEQ2];
3691 t_r=CL->residue_index[s1][r1][b+R2];
3692 for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3694 q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3695 q_r=CL->residue_index[t_s][t_r][c+R2];
3697 n_hasch1[q_s][q_r]=0;
3704 score=(score*SCORE_K)/((ns1-n_gap1)*(ns2-n_gap2));
3705 score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
3707 return (int)(score-SCORE_K*CL->nomatch);
3711 int fast_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3713 /*WARNING: WORKS ONLY WITH List to Extend*/
3721 int s1, rs1, r1, t_r, t_s, s2, rs2, r2;
3722 static int **hasch1;
3723 static int **hasch2;
3725 static int **n_hasch1;
3726 static int **n_hasch2;
3728 static int **is_in_col1;
3729 static int **is_in_col2;
3737 hasch1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3738 hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3739 n_hasch1=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
3740 n_hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3741 is_in_col1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3742 is_in_col2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3745 for ( a=0; a< ns1; a++)
3748 s1=A->order[rs1][0];
3754 is_in_col1[s1][r1]=1;
3755 for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3757 t_s=CL->residue_index[s1][r1][b+SEQ2];
3758 t_r=CL->residue_index[s1][r1][b+R2];
3759 hasch1[t_s][t_r]+=CL->residue_index[s1][r1][b+WE];
3760 n_hasch1[t_s][t_r]++;
3766 for ( a=0; a< ns2; a++)
3769 s2=A->order[rs2][0];
3775 is_in_col2[s2][r2]=1;
3776 for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3778 t_s=CL->residue_index[s2][r2][b+SEQ2];
3779 t_r=CL->residue_index[s2][r2][b+R2];
3781 hasch2[t_s][t_r]+=CL->residue_index[s2][r2][b+WE];
3782 n_hasch2[t_s][t_r]++;
3790 for ( a=0; a< ns2; a++)
3793 s2=A->order[rs2][0];
3799 for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3801 t_s=CL->residue_index[s2][r2][b+SEQ2];
3802 t_r=CL->residue_index[s2][r2][b+R2];
3804 if ( hasch2[t_s][t_r] && hasch1[t_s][t_r]&& !(is_in_col1[t_s][t_r] || is_in_col2[t_s][t_r]))
3806 score+=MIN(hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]),hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]));
3808 else if ( hasch2[t_s][t_r] && is_in_col1[t_s][t_r])
3810 score+=hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]+1);
3812 else if (hasch1[t_s][t_r] && is_in_col2[t_s][t_r])
3814 score+=hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]+1);
3817 n_hasch2[t_s][t_r]=0;
3820 is_in_col2[s2][r2]=0;
3825 for ( a=0; a< ns1; a++)
3828 s1=A->order[rs1][0];
3834 is_in_col1[s1][r1]=0;
3836 for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3838 t_s=CL->residue_index[s1][r1][b+SEQ2];
3839 t_r=CL->residue_index[s1][r1][b+R2];
3842 n_hasch1[t_s][t_r]=0;
3849 for ( a=0; a< ns1; a++)
3852 s1=A->order[rs1][0];
3858 for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3860 t_s=CL->residue_index[s1][r1][b+SEQ2];
3861 t_r=CL->residue_index[s1][r1][b+R2];
3863 if ( hasch1[t_s][t_r] && hasch2[t_s][t_r]&& !(is_in_col2[t_s][t_r] || is_in_col1[t_s][t_r]))
3865 score+=MIN(hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]),hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]));
3867 else if ( hasch1[t_s][t_r] && is_in_col2[t_s][t_r])
3869 score+=hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]+1);
3871 else if (hasch2[t_s][t_r] && is_in_col1[t_s][t_r])
3873 score+=hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]+1);
3876 n_hasch1[t_s][t_r]=0;
3879 is_in_col1[s1][r1]=0;
3884 for ( a=0; a< ns2; a++)
3887 s2=A->order[rs2][0];
3893 is_in_col2[s2][r2]=0;
3895 for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3897 t_s=CL->residue_index[s2][r2][b+SEQ2];
3898 t_r=CL->residue_index[s2][r2][b+R2];
3901 n_hasch2[t_s][t_r]=0;
3906 score=(score*SCORE_K)/((ns1-n_gap1)*(ns2-n_gap2));
3907 score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
3909 return (int)(score-SCORE_K*CL->nomatch);
3912 int fast_get_dp_cost_2 ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3916 int a, b, s1, s2,r1, r2;
3922 static int last_ns1;
3923 static int last_ns2;
3924 static int last_top1;
3925 static int last_top2;
3926 static int **check_list;
3929 /*New heuristic get dp_cost on a limited number of pairs*/
3930 /*This is the current default*/
3933 if ( last_ns1==ns1 && last_top1==list1[0] && last_ns2==ns2 && last_top2==list2[0]);
3942 if ( check_list) free_int (check_list, -1);
3943 check_list=declare_int ( (CL->S)->nseq*(CL->S)->nseq, 3);
3945 for ( n_pair=0,a=0; a< ns1; a++)
3948 rs1=A->order[s1][0];
3949 for ( b=0; b< ns2; b++, n_pair++)
3952 rs2=A->order[s2][0];
3953 check_list[n_pair][0]=s1;
3954 check_list[n_pair][1]=s2;
3955 check_list[n_pair][2]=(!CL->DM)?0:(CL->DM)->similarity_matrix[rs1][rs2];
3957 sort_int ( check_list, 3, 2, 0, n_pair-1);
3961 if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair;}
3964 for ( a=n_pair-1; a>=0; a--)
3966 s1= check_list[a][0];
3967 rs1=A->order[s1][0];
3970 s2= check_list[a][1];
3971 rs2=A->order[s2][0];
3982 if ((s=(CL->evaluate_residue_pair)(CL, rs1, r1, rs2, r2))!=UNDEFINED) score+=s;
3988 if ( res_res>=CL->max_n_pair && CL->max_n_pair!=0)a=0;
3991 score=(res_res==0)?0:( (score)/res_res);
3992 score=score-SCORE_K*CL->nomatch;
4000 int slow_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4004 int a, b, s1, s2,r1, r2;
4010 static int last_tag;
4013 if (col1<0 || col2<0) return 0;
4014 if ( last_tag !=A->random_tag)
4016 last_tag=A->random_tag;
4019 dummy=vcalloc (10, sizeof(int));
4020 dummy[0]=1;/*Number of Amino acid types on colum*/
4021 dummy[1]=5;/*Length of Dummy*/
4022 dummy[3]='\0';/*Amino acid*/
4023 dummy[4]=1; /*Number of occurences*/
4024 dummy[5]=100; /*Frequency in the MSA column*/
4028 if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair;}
4030 for ( a=0; a< ns1; a++)
4032 for ( b=0; b<ns2; b++)
4036 rs1=A->order[s1][0];
4040 rs2=A->order[s2][0];
4049 if (r1==0 && r2==0)gap_gap++;
4050 else if ( r1<0 || r2<0) gap_res++;
4054 if ((s=(CL->evaluate_residue_pair)(CL, rs1, r1, rs2, r2))!=UNDEFINED) score+=s;
4066 score=(res_res==0)?0:( (score)/res_res);
4068 return score-SCORE_K*CL->nomatch;
4071 int slow_get_dp_cost_pc ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4075 int a, b, s1, s2,r1, r2;
4081 static int last_tag;
4084 if (col1<0 || col2<0) return 0;
4085 if ( last_tag !=A->random_tag)
4087 last_tag=A->random_tag;
4090 dummy=vcalloc (10, sizeof(int));
4091 dummy[0]=1;/*Number of Amino acid types on colum*/
4092 dummy[1]=5;/*Length of Dummy*/
4093 dummy[3]='\0';/*Amino acid*/
4094 dummy[4]=1; /*Number of occurences*/
4095 dummy[5]=100; /*Frequency in the MSA column*/
4099 if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair;}
4101 for ( a=0; a< ns1; a++)
4103 for ( b=0; b<ns2; b++)
4107 rs1=A->order[s1][0];
4111 rs2=A->order[s2][0];
4120 if (r1==0 && r2==0)gap_gap++;
4121 else if ( r1<0 || r2<0) gap_res++;
4125 if ((s=residue_pair_extended_list_pc(CL, rs1, r1, rs2, r2))!=UNDEFINED) score+=s;
4137 score=(res_res==0)?0:( (score)/res_res);
4142 int slow_get_dp_cost_test ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4146 int a, b, s1, s2,r1, r2;
4147 int gap_gap=0, gap_res=0, res_res=0, rs1, rs2;
4149 for ( a=0; a< ns1; a++)
4151 for ( b=0; b<ns2; b++)
4155 rs1=A->order[s1][0];
4159 rs2=A->order[s2][0];
4168 if (r1==0 && r2==0)gap_gap++;
4169 else if ( r1<0 || r2<0) gap_res++;
4173 score+=residue_pair_extended_list_raw (CL, rs1, r1, rs2, r2);
4178 return (int)(score*10)/(ns1*ns2);
4181 int sw_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4190 for ( a=0; a< ns1; a++)
4193 rs1=A->order[s1][0];
4195 if ( r1<=0)continue;
4196 for ( b=0; b< ns2; b++)
4201 rs2=A->order[s2][0];
4207 if (sw_pair_is_defined (CL, rs1, r1, rs2, r2)==UNDEFINED)return UNDEFINED;
4211 return slow_get_dp_cost ( A, pos1, ns1, list1, col1, pos2, ns2, list2,col2, CL);
4222 int get_domain_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2,Constraint_list *CL , int scale , int gop, int gep)
4225 int a, b, s1, s2,r1, r2;
4233 int flag_list_is_aa_sub_mat=0;
4236 /*Needs to be cleanned After Usage*/
4240 if ( entry==NULL) entry=vcalloc (LIST_N_FIELDS , sizeof (int));
4242 for (a=0; a< ns1; a++)
4245 rs1=A->order[s1][0];
4246 for ( b=0; b<ns2; b++)
4249 rs2=A->order[s2][0];
4253 r1=entry[R1]=pos1[s1][col1];
4254 r2=entry[R2]=pos2[s2][col2];
4256 if ( !flag_list_is_aa_sub_mat)
4258 if ( r1==r2 && rs1==rs2)
4263 else if (r1==0 && r2==0)
4267 else if ( r1<=0 || r2<=0)
4271 else if ((r=main_search_in_list_constraint ( entry,&p,4,CL))!=NULL)
4275 if (r[WE]!=UNDEFINED)
4277 score+=(r[WE]*SCORE_K)+scale;
4281 fprintf ( stderr, "**");
4289 score=((res_res+gap_res)==0)?0:score/(res_res+gap_res);
4293 /*********************************************************************************************/
4295 /* FUNCTIONS FOR ANALYSING AL OR MATRIX */
4297 /*********************************************************************************************/
4299 int aln2n_res ( Alignment *A, int start, int end)
4304 for ( a=start; a<end; a++)for ( b=0; b< A->nseq; b++)score+=!is_gap(A->seq_al[b][a]);
4308 float get_gop_scaling_factor ( int **matrix,float id, int l1, int l2)
4310 return id* get_avg_matrix_mm(matrix, AA_ALPHABET);
4313 float get_avg_matrix_mm ( int **matrix, char *alphabet)
4321 l=MIN(20,(int)strlen (alphabet));
4322 for (naa=0, gop=0,a=0; a<l; a++)
4323 for ( b=0; b<l; b++)
4327 gop+=matrix[alphabet[a]-'A'][alphabet[b]-'A'];
4335 float get_avg_matrix_match ( int **matrix, char *alphabet)
4344 l=MIN(20,(int)strlen (alphabet));
4345 for (gop=0,a=0; a<l; a++)
4346 gop+=matrix[alphabet[a]-'A'][alphabet[a]-'A'];
4352 float get_avg_matrix_diff ( int **matrix1,int **matrix2, char *alphabet)
4362 l=MIN(20,(int)strlen (alphabet));
4363 for (naa=0, gop=0,a=0; a<l; a++)
4366 v1=matrix1[alphabet[a]-'A'][alphabet[a]-'A']-matrix2[alphabet[b]-'A'][alphabet[b]-'A'];
4375 float* set_aa_frequencies ()
4379 /*frequencies tqken from psw*/
4381 frequency=vcalloc (100, sizeof (float));
4382 frequency ['x'-'A']=0.0013;
4383 frequency ['a'-'A']=0.0076;
4384 frequency ['c'-'A']=0.0176;
4385 frequency ['d'-'A']=0.0529;
4386 frequency ['e'-'A']=0.0628;
4387 frequency ['f'-'A']=0.0401;
4388 frequency ['g'-'A']=0.0695;
4389 frequency ['h'-'A']=0.0224;
4390 frequency ['i'-'A']=0.0561;
4391 frequency ['k'-'A']=0.0584;
4392 frequency ['l'-'A']=0.0922;
4393 frequency ['m'-'A']=0.0236;
4394 frequency ['n'-'A']=0.0448;
4395 frequency ['p'-'A']=0.0500;
4396 frequency ['q'-'A']=0.0403;
4397 frequency ['r'-'A']=0.0523;
4398 frequency ['s'-'A']=0.0715;
4399 frequency ['t'-'A']=0.0581;
4400 frequency ['v'-'A']=0.0652;
4401 frequency ['w'-'A']=0.0128;
4402 frequency ['y'-'A']=0.0321;
4403 frequency ['X'-'A']=0.0013;
4404 frequency ['A'-'A']=0.0076;
4405 frequency ['C'-'A']=0.0176;
4406 frequency ['D'-'A']=0.0529;
4407 frequency ['E'-'A']=0.0628;
4408 frequency ['F'-'A']=0.0401;
4409 frequency ['G'-'A']=0.0695;
4410 frequency ['H'-'A']=0.0224;
4411 frequency ['I'-'A']=0.0561;
4412 frequency ['J'-'A']=0.0584;
4413 frequency ['L'-'A']=0.0922;
4414 frequency ['M'-'A']=0.0236;
4415 frequency ['N'-'A']=0.0448;
4416 frequency ['P'-'A']=0.0500;
4417 frequency ['Q'-'A']=0.0403;
4418 frequency ['R'-'A']=0.0523;
4419 frequency ['S'-'A']=0.0715;
4420 frequency ['T'-'A']=0.0581;
4421 frequency ['V'-'A']=0.0652;
4422 frequency ['W'-'A']=0.0128;
4423 frequency ['Y'-'A']=0.0321;
4427 float measure_matrix_pos_avg (int **matrix,char *alphabet)
4433 for ( tot=0,a=0; a< 20; a++)
4434 for ( b=0; b<20; b++)
4436 if (matrix[alphabet[a]-'A'][alphabet[b]-'A']>=0){naa++;tot+=matrix[alphabet[a]-'A'][alphabet[b]-'A'];}
4442 float measure_matrix_enthropy (int **matrix,char *alphabet)
4446 double s, p, q, h=0, tq=0;
4449 /*frequencies tqken from psw*/
4451 frequency=set_aa_frequencies ();
4454 lambda=compute_lambda(matrix,alphabet);
4455 fprintf ( stderr, "\nLambda=%f", (float)lambda);
4457 for ( a=0; a< 20; a++)
4458 for ( b=0; b<=a; b++)
4460 s=matrix[alphabet[a]-'A'][alphabet[b]-'A'];
4463 p=frequency[alphabet[a]-'A']*frequency[alphabet[b]-'A'];
4467 q=exp(lambda*s+log(p));
4470 h+=q*log(q/p)*log(2);
4474 fprintf ( stderr,"\ntq=%f\n", (float)tq);
4478 float compute_lambda (int **matrix,char *alphabet)
4482 double lambda, best_lambda=0, delta, best_delta=0, p, tq,s;
4483 static float *frequency;
4485 if ( !frequency)frequency=set_aa_frequencies ();
4487 for ( lambda=0.001; lambda<1; lambda+=0.005)
4490 for ( a=0; a< 20; a++)
4491 for ( b=0; b<20; b++)
4493 p=frequency[alphabet[a]-'A']*frequency[alphabet[b]-'A'];
4494 s=matrix[alphabet[a]-'A'][alphabet[b]-'A'];
4495 tq+=exp(lambda*s+log(p));
4505 if (delta<best_delta)
4512 fprintf ( stderr, "\n%f %f ", lambda, tq);
4515 fprintf ( stderr, "\nRESULT: %f %f ", best_lambda, best_delta);
4516 return (float) best_lambda;
4521 float evaluate_random_match (char *mat, int n, int len,char *alp)
4524 matrix=read_matrice ( mat);
4525 fprintf ( stderr, "Matrix=%15s ", mat);
4526 return evaluate_random_match2 (matrix, n,len,alp);
4530 float evaluate_random_match2 (int **matrix, int n, int len,char *alp)
4532 int a, b, c, d, c1, c2, tot;
4535 float score_random=0;
4545 freq=set_aa_frequencies ();
4546 list=vcalloc ( 10000, sizeof (char));
4549 for (tot=0,c=0,a=0;a<20; a++)
4551 b=freq[alp[a]-'A']*1000;
4553 for (d=0; d<b; d++, c++)
4560 for (a=0; a< len*n; a++)
4564 score_random+=matrix[list[c1]-'A'][list[c2]-'A'];
4565 score_id+=matrix[list[c1]-'A'][list[c1]-'A'];
4567 while (tot_len< len*n)
4572 if ( matrix[list[c1]-'A'][list[c2]-'A']>=0){score_good+=matrix[list[c1]-'A'][list[c2]-'A']; tot_len++;}
4576 score_random=score_random/tot_len;
4577 score_id=score_id/tot_len;
4578 score_good=score_good/tot_len;
4580 fprintf ( stderr, "Random=%8.3f Id=%8.3f Good=%8.3f [%7.2f]\n",score_random, score_id, score_good, tot_len/tot_try);
4582 return score_random;
4584 float compare_two_mat (char *mat1,char*mat2, int n, int len,char *alp)
4586 int **matrix1, **matrix2;
4588 evaluate_random_match (mat1, n, len,alp);
4589 evaluate_random_match (mat2, n, len,alp);
4590 matrix1=read_matrice ( mat1);
4591 matrix2=read_matrice ( mat2);
4592 matrix1=rescale_matrix(matrix1, 10, alp);
4593 matrix2=rescale_matrix(matrix2, 10, alp);
4594 compare_two_mat_array(matrix1,matrix2, n, len,alp);
4599 int ** rescale_two_mat (char *mat1,char*mat2, int n, int len,char *alp)
4602 int **matrix1, **matrix2;
4604 lambda=measure_lambda2 (mat1, mat2, n, len, alp)*10;
4606 fprintf ( stderr, "\nLambda=%.2f", lambda);
4607 matrix2=read_matrice(mat2);
4608 matrix2=neg_matrix2pos_matrix(matrix2);
4609 matrix2=rescale_matrix( matrix2, lambda,"abcdefghiklmnpqrstvwxyz");
4611 matrix1=read_matrice(mat1);
4612 matrix1=neg_matrix2pos_matrix(matrix1);
4613 matrix1=rescale_matrix( matrix1,10,"abcdefghiklmnpqrstvwxyz");
4615 output_matrix_header ( "stdout", matrix2, alp);
4616 evaluate_random_match2(matrix1, 1000, 100, alp);
4617 evaluate_random_match2(matrix2, 1000, 100, alp);
4618 compare_two_mat_array(matrix1,matrix2, n, len,alp);
4622 float measure_lambda2(char *mat1,char*mat2, int n, int len,char *alp)
4627 m1=read_matrice (mat1);
4628 m2=read_matrice (mat2);
4630 m1=neg_matrix2pos_matrix(m1);
4631 m2=neg_matrix2pos_matrix(m2);
4633 f1=measure_matrix_pos_avg( m1, alp);
4634 f2=measure_matrix_pos_avg( m2, alp);
4640 float measure_lambda (char *mat1,char*mat2, int n, int len,char *alp)
4643 int **matrix1, **matrix2, **mat;
4645 float best_quality=0, quality=0, best_lambda=0;
4647 matrix1=read_matrice ( mat1);
4648 matrix2=read_matrice ( mat2);
4649 matrix1=rescale_matrix(matrix1, 10, alp);
4650 matrix2=rescale_matrix(matrix2, 10, alp);
4652 for (c=0, a=0.1; a< 2; a+=0.05)
4654 fprintf ( stderr, "Lambda=%.2f\n", a);
4655 mat=duplicate_int (matrix2,-1,-1);
4656 mat=rescale_matrix(mat, a, alp);
4657 quality=compare_two_mat_array(matrix1,mat, n, len,alp);
4658 quality=MAX((-quality),quality);
4660 if (c==0 || (best_quality>quality))
4663 fprintf ( stderr, "*");
4664 best_quality=quality;
4669 evaluate_random_match2(mat, 1000, 100, alp);
4670 evaluate_random_match2(matrix1, 1000, 100, alp);
4678 float compare_two_mat_array (int **matrix1,int **matrix2, int n, int len,char *alp)
4680 int a, b, c, d, c1, c2, tot;
4683 float delta_random=0;
4684 float delta2_random=0;
4690 float delta2_good=0;
4702 freq=set_aa_frequencies ();
4703 list=vcalloc ( 10000, sizeof (char));
4706 for (tot=0,c=0,a=0;a<20; a++)
4708 b=freq[alp[a]-'A']*1000;
4710 for (d=0; d<b; d++, c++)
4717 for (a=0; a< len*n; a++)
4721 delta=matrix1[list[c1]-'A'][list[c2]-'A']-matrix2[list[c1]-'A'][list[c2]-'A'];
4722 delta_random+=delta;
4723 delta2_random+=MAX(delta,(-delta));
4725 delta=matrix1[list[c1]-'A'][list[c1]-'A']-matrix2[list[c1]-'A'][list[c1]-'A'];
4727 delta2_id+=MAX(delta,(-delta));
4729 while (tot_len< len*n)
4734 if ( matrix1[list[c1]-'A'][list[c2]-'A']>=0 || matrix2[list[c1]-'A'][list[c2]-'A'] )
4736 delta=matrix1[list[c1]-'A'][list[c2]-'A']-matrix2[list[c1]-'A'][list[c2]-'A'];
4738 delta2_good+=MAX(delta,(-delta));
4744 delta_random=delta_random/tot_len;
4745 delta2_random=delta2_random/tot_len;
4748 delta_id=delta_id/tot_len;
4749 delta2_id=delta2_id/tot_len;
4751 delta_good=delta_good/tot_len;
4752 delta2_good=delta2_good/tot_len;
4755 fprintf ( stderr, "\tRand=%8.3f %8.3f\n\tId =%8.3f %8.3f\n\tGood=%8.3f %8.3f\n",delta_random, delta2_random, delta_id,delta2_id, delta_good,delta2_good);
4762 int ** rescale_matrix ( int **matrix, float lambda, char *alp)
4767 for ( a=0; a< 20; a++)
4768 for ( b=0; b< 20; b++)
4770 matrix[alp[a]-'A'][alp[b]-'A']= matrix[alp[a]-'A'][alp[b]-'A']*lambda;
4774 int **mat2inverted_mat (int **matrix, char *alp)
4776 int a, b, min, max, v,l;
4779 l=(int)strlen (alp);
4780 min=max=matrix[alp[0]-'A'][alp[0]-'A'];
4781 for ( a=0; a<l; a++)
4782 for ( b=0; b< l; b++)
4784 v=matrix[alp[a]-'A'][alp[b]-'A'];
4788 for ( a=0; a<l; a++)
4789 for ( b=0; b< l; b++)
4791 v=matrix[alp[a]-'A'][alp[b]-'A'];
4794 c1=tolower(alp[a])-'A';
4795 c2=tolower(alp[b])-'A';
4797 C1=toupper(alp[a])-'A';
4798 C2=toupper(alp[b])-'A';
4799 matrix[C1][C2]=matrix[c1][c2]=matrix[C1][c2]=matrix[c1][C2]=v;
4803 void output_matrix_header ( char *name, int **matrix, char *alp)
4809 nalp=vcalloc ( 1000, sizeof (char));
4813 sprintf ( nalp, "ABCDEFGHIKLMNPQRSTVWXYZ");
4814 l=(int)strlen (nalp);
4816 fp=vfopen ( name, "w");
4817 fprintf (fp, "\nint []={\n");
4821 for ( b=0; b<=a; b++)
4822 fprintf ( fp, "%3d, ",matrix[nalp[a]-'A'][nalp[b]-'A']);
4823 fprintf ( fp, "\n");
4825 fprintf (fp, "};\n");
4831 float ** initialise_aa_physico_chemical_property_table (int *n)
4841 fp=vfopen ("properties.txt", "r");
4842 while ((c=fgetc(fp))!=EOF)n[0]+=(c=='#');
4847 p=declare_float ( n[0], 'z'+1);
4848 fp=vfopen ("properties.txt", "r");
4850 while ((c=fgetc(fp))!=EOF)
4852 if (c=='#'){in=0;while (fgetc(fp)!='\n');}
4855 if ( !in){in=1; ++n[0];}
4856 fscanf (fp, "%f %c %*s",&v1, &c1 );
4857 p[n[0]-1][tolower(c1)]=v1;
4858 while ( (c=fgetc(fp))!='\n');
4863 /*rescale so that Delta max=10*/
4864 for (a=0; a<n[0]; a++)
4867 for (b='a'; b<='z'; b++)
4869 min=(p[a][b]<min)?p[a][b]:min;
4870 max=(p[a][b]>max)?p[a][b]:max;
4872 for (b='a'; b<='z'; b++)
4874 p[a][b]=((p[a][b]-min)/(max-min))*10;
4881 Constraint_list * choose_extension_mode ( char *extend_mode, Constraint_list *CL)
4883 //evaluation_functions: residues start at 1, sequences at 0;
4887 fprintf ( stderr, "\nWarning: CL was not set");
4890 else if ( strm ( extend_mode, "rna0"))
4892 CL->evaluate_residue_pair=residue_pair_extended_list;
4893 CL->get_dp_cost =slow_get_dp_cost;
4895 else if ( strm ( extend_mode, "rna1") || strm (extend_mode, "rna"))
4897 CL->evaluate_residue_pair=residue_pair_extended_list4rna1;
4898 CL->get_dp_cost =slow_get_dp_cost;
4900 else if ( strm ( extend_mode, "rna2"))
4902 CL->evaluate_residue_pair=residue_pair_extended_list4rna2;
4903 CL->get_dp_cost =slow_get_dp_cost;
4905 else if ( strm ( extend_mode, "rna3"))
4907 CL->evaluate_residue_pair=residue_pair_extended_list4rna3;
4908 CL->get_dp_cost =slow_get_dp_cost;
4910 else if ( strm ( extend_mode, "rna4"))
4912 CL->evaluate_residue_pair=residue_pair_extended_list4rna4;
4913 CL->get_dp_cost =slow_get_dp_cost;
4915 else if ( strm ( extend_mode, "pc") && !CL->M)
4917 CL->evaluate_residue_pair=residue_pair_extended_list_pc;
4918 CL->get_dp_cost =slow_get_dp_cost;
4920 else if ( strm ( extend_mode, "triplet") && !CL->M)
4922 CL->evaluate_residue_pair=residue_pair_extended_list;
4923 CL->get_dp_cost =get_dp_cost;
4925 else if ( strm ( extend_mode, "relative_triplet") && !CL->M)
4927 CL->evaluate_residue_pair=residue_pair_relative_extended_list;
4928 CL->get_dp_cost =fast_get_dp_cost_2;
4930 else if ( strm ( extend_mode, "g_coffee") && !CL->M)
4932 CL->evaluate_residue_pair=residue_pair_extended_list_g_coffee;
4933 CL->get_dp_cost =slow_get_dp_cost;
4935 else if ( strm ( extend_mode, "g_coffee_quadruplets") && !CL->M)
4937 CL->evaluate_residue_pair=residue_pair_extended_list_g_coffee_quadruplet;
4938 CL->get_dp_cost =slow_get_dp_cost;
4940 else if ( strm ( extend_mode, "fast_triplet") && !CL->M)
4942 CL->evaluate_residue_pair=residue_pair_extended_list;
4943 CL->get_dp_cost =fast_get_dp_cost;
4946 else if ( strm ( extend_mode, "very_fast_triplet") && !CL->M)
4948 CL->evaluate_residue_pair=residue_pair_extended_list;
4949 CL->get_dp_cost =fast_get_dp_cost_2;
4951 else if ( strm ( extend_mode, "slow_triplet") && !CL->M)
4953 CL->evaluate_residue_pair=residue_pair_extended_list;
4954 CL->get_dp_cost =slow_get_dp_cost;
4956 else if ( strm ( extend_mode, "mixt") && !CL->M)
4958 CL->evaluate_residue_pair=residue_pair_extended_list_mixt;
4959 CL->get_dp_cost=slow_get_dp_cost;
4961 else if ( strm ( extend_mode, "quadruplet") && !CL->M)
4963 CL->evaluate_residue_pair=residue_pair_extended_list_quadruplet;
4964 CL->get_dp_cost =get_dp_cost_quadruplet;
4966 else if ( strm ( extend_mode, "test") && !CL->M)
4968 CL->evaluate_residue_pair=residue_pair_test_function;
4969 CL->get_dp_cost =slow_get_dp_cost_test;
4971 else if ( strm ( extend_mode, "ssp"))
4974 CL->evaluate_residue_pair=evaluate_ssp_matrix_score;
4975 CL->get_dp_cost=slow_get_dp_cost;
4978 else if ( strm ( extend_mode, "tm"))
4981 CL->evaluate_residue_pair=evaluate_tm_matrix_score;
4982 CL->get_dp_cost=slow_get_dp_cost;
4985 else if ( strm ( extend_mode, "matrix"))
4988 CL->evaluate_residue_pair=evaluate_matrix_score;
4989 CL->get_dp_cost=cw_profile_get_dp_cost;
4992 else if ( strm ( extend_mode, "curvature"))
4994 CL->evaluate_residue_pair=evaluate_curvature_score;
4995 CL->get_dp_cost=slow_get_dp_cost;
5000 CL->evaluate_residue_pair=evaluate_matrix_score;
5001 CL->get_dp_cost=cw_profile_get_dp_cost;
5006 fprintf ( stderr, "\nERROR: %s is an unknown extend_mode[FATAL:%s]\n", extend_mode, PROGRAM);
5007 myexit (EXIT_FAILURE);
5012 int ** combine_two_matrices ( int **mat1, int **mat2)
5014 int naa, re1, re2, Re1, Re2, a, b, u, l;
5016 naa=(int)strlen (BLAST_AA_ALPHABET);
5017 for ( a=0; a< naa; a++)
5018 for ( b=0; b< naa; b++)
5020 re1=BLAST_AA_ALPHABET[a];
5021 re2=BLAST_AA_ALPHABET[b];
5022 if (re1=='*' || re2=='*');
5026 Re1=toupper(re1);Re2=toupper(re2);
5027 re1-='A';re2-='A';Re1-='A';Re2-='A';
5031 mat1[re1][re2]=mat2[re1][re2]=l;
5032 mat2[Re1][Re2]=mat2[Re1][Re2]=u;
5038 /* Off the shelves evaluations */
5039 /*********************************************************************************************/
5041 /* OFF THE SHELVES EVALUATION */
5043 /*********************************************************************************************/
5046 int lat_sum_pair (Alignment *A, char *mat)
5048 int a,b,c, tot=0, v1, v2, score;
5051 matrix=read_matrice (mat);
5053 for (a=0; a<A->nseq; a++)
5054 for ( b=0; b<A->nseq; b++)
5056 for (c=1; c<A->len_aln; c++)
5060 r11=A->seq_al[a][c-1];
5061 r12=A->seq_al[a][c];
5062 if (is_gap(r11) || is_gap(r12))continue;
5063 else v1=matrix[r11-'A'][r12-'A'];
5065 r11=A->seq_al[b][c-1];
5066 r12=A->seq_al[b][c];
5067 if (is_gap(r11) || is_gap(r12))continue;
5068 else v2=matrix[r11-'A'][r12-'A'];
5070 score+=(v1-v2)*(v1-v2);
5074 score=(100*score)/tot;
5075 return (float)score;
5080 /* Off the shelves evaluations */
5081 /*********************************************************************************************/
5083 /* OFF THE SHELVES EVALUATION */
5085 /*********************************************************************************************/
5087 int comp_pair ( int len,char *sa, char *sb, int seqA, int seqB,int *tgp_a, int *tgp_b,int gap_op,int gap_ex, int start, int end,int **matrix,int MODE);
5088 int score_gap ( int len, char *sa, char *sb,int seqA, int seqB,int *tgp_a, int *tgp_b, int op, int ex, int start, int end, int MODE);
5089 void evaluate_tgp_decoded_chromosome ( Alignment *A,int **TGP,int start, int end,int MODE);
5090 int gap_type ( char a, char b);
5094 float sum_pair ( Alignment*A,char *mat_name, int gap_op, int gap_ex)
5105 matrix=read_matrice (mat_name);
5106 matrix=mat2inverted_mat (matrix, "acdefghiklmnpqrstvwy");
5112 tgp= declare_int (A->nseq,2);
5114 evaluate_tgp_decoded_chromosome ( A,tgp,start, end,MODE);
5116 for ( a=0; a< A->nseq-1; a++)
5117 for (b=a+1; b<A->nseq; b++)
5119 pscore= comp_pair (A->len_aln,A->seq_al[a], A->seq_al[b],a, b,tgp[a], tgp[b],gap_op,gap_ex, start, end,matrix, MODE);
5121 /*score+=(double)pscore*(int)(PARAM->OFP)->weight[A->order[a][0]][A->order[b][0]];*//*NO WEIGHTS*/
5124 score=score/(A->nseq*A->nseq);
5126 return (float)score;
5129 int comp_pair ( int len,char *sa, char *sb, int seqA, int seqB,int *tgp_a, int *tgp_b,int gap_op,int gap_ex, int start, int end,int **matrix,int MODE)
5136 score+= score_gap (len, sa,sb, seqA, seqB,tgp_a, tgp_b, gap_op,gap_ex, start, end,MODE);
5141 for (a=start; a<=end; a++)
5143 if ( is_gap(sa[a]) || is_gap(sb[a]))
5145 if (is_gap(sa[a]) && is_gap(sb[a]));
5154 score += matrix [sa[a]-'A'][sb[a]-'A'];
5160 int score_gap ( int len, char *sa, char *sb,int seqA, int seqB,int *tgp_a, int *tgp_b, int op, int ex, int start, int end, int MODE)
5167 int right_gap, left_gap;
5177 int sequence_pattern[2][3];
5181 /*op= gor_gap_op ( 0,seqA, seqB, PARAM);
5182 ex= gor_gap_ext ( 0, seqA, seqB, PARAM);*/
5186 for (a=start; a<=end; ++a)
5189 type= gap_type ( sa[a], sb[a]);
5191 if ( type==2 && ga<=gb)
5196 else if (type==1 && ga >=gb)
5218 else if ( (type==flag1) && flag2==1)
5223 else if ( (type!=flag1) && flag2==1)
5232 /*gap_type -/-:0, X/X:-1 X/-:1, -/X:2*/
5233 /*evaluate the pattern of gaps*/
5236 sequence_pattern[0][0]=sequence_pattern[1][0]=0;
5237 for ( a=start; a<=end && continue_loop==1; a++)
5239 left_gap= gap_type ( sa[a], sb[a]);
5244 sequence_pattern[0][0]=sequence_pattern[1][0]=0;
5250 for (b=a; b<=end && continue_loop==1; b++)
5251 {type=gap_type( sa[b], sb[b]);
5254 if ( type!=left_gap && type !=0)
5257 sequence_pattern[2-left_gap][0]= b-a-null_gap;
5258 sequence_pattern [1-(2-left_gap)][0]=0;
5261 if ( continue_loop==1)
5264 sequence_pattern[2-left_gap][0]= b-a-null_gap;
5265 sequence_pattern [1-(2-left_gap)][0]=0;
5271 sequence_pattern[0][2]=sequence_pattern[1][2]=1;
5272 for ( a=start; a<=end; a++)
5274 if ( !is_gap(sa[a]))
5275 sequence_pattern[0][2]=0;
5276 if ( !is_gap(sb[a]))
5277 sequence_pattern[1][2]=0;
5281 sequence_pattern[0][1]=sequence_pattern[1][1]=0;
5282 for ( a=end; a>=start && continue_loop==1; a--)
5284 right_gap= gap_type ( sa[a], sb[a]);
5289 sequence_pattern[0][1]=sequence_pattern[1][1]=0;
5295 for (b=a; b>=start && continue_loop==1; b--)
5296 {type=gap_type( sa[b], sb[b]);
5299 if ( type!=right_gap && type !=0)
5302 sequence_pattern[2-right_gap][1]= a-b-null_gap;
5303 sequence_pattern [1-(2-right_gap)][1]=0;
5306 if ( continue_loop==1)
5309 sequence_pattern[2-right_gap][1]= a-b-null_gap;
5310 sequence_pattern [1-(2-right_gap)][1]=0;
5317 printf ( "\n*****************************************************");
5318 printf ( "\n%c\n%c", sa[start],sb[start]);
5319 printf ( "\n%d %d %d",sequence_pattern[0][0] ,sequence_pattern[0][1], sequence_pattern[0][2]);
5320 printf ( "\n%d %d %d",sequence_pattern[1][0] ,sequence_pattern[1][1], sequence_pattern[1][2]);
5321 printf ( "\n*****************************************************");
5324 /*correct the scoring*/
5329 if ( FABS(tgp_a[0])>1 && (FABS(tgp_a[0])>FABS( tgp_b[0])))
5330 score-= (sequence_pattern[0][0]>0)?op:0;
5331 if ( FABS(tgp_b[0])>1 && (FABS(tgp_b[0])> FABS(tgp_a[0])))
5332 score-= (sequence_pattern[1][0]>0)?op:0;
5334 else if ( MODE ==1 || MODE ==2)
5336 if ( FABS(tgp_a[0])>1 && (FABS(tgp_a[0])>FABS( tgp_b[0])) && (tgp_a[1]!=1 || sequence_pattern[0][2]==0))
5337 score-= (sequence_pattern[0][0]>0)?op:0;
5338 if ( FABS(tgp_b[0])>1 && (FABS(tgp_b[0])> FABS(tgp_a[0])) && (tgp_b[1]!=1 || sequence_pattern[1][2]==0))
5339 score-= (sequence_pattern[1][0]>0)?op:0;
5342 if ( tgp_a[0]>=1 && tgp_a[0]==tgp_b[0])
5343 score -=(sequence_pattern[0][0]>0)?op:0;
5344 if ( tgp_b[0]>=1 && tgp_a[0]==tgp_b[0])
5345 score-= (sequence_pattern[1][0]>0)?op:0;
5348 if ( tgp_a[1]==1 && sequence_pattern[0][2]==0)
5349 score -= ( sequence_pattern[0][1]>0)?op:0;
5350 else if ( tgp_a[1]==1 && sequence_pattern[0][2]==1 && tgp_a[0]<=0)
5351 score -= ( sequence_pattern[0][1]>0)?op:0;
5354 if ( tgp_b[1]==1 && sequence_pattern[1][2]==0)
5355 score -= ( sequence_pattern[1][1]>0)?op:0;
5356 else if ( tgp_b[1]==1 && sequence_pattern[1][2]==1 && tgp_b[0]<=0)
5357 score -= ( sequence_pattern[1][1]>0)?op:0;
5362 score -=sequence_pattern[0][0]*ex;
5364 score -= sequence_pattern[1][0]*ex;
5366 score-=sequence_pattern[0][1]*ex;
5368 score-=sequence_pattern[1][1]*ex;
5378 void evaluate_tgp_decoded_chromosome ( Alignment *A,int **TGP,int start, int end,int MODE)
5385 if (MODE==11 || MODE==13|| MODE==14)
5387 if ( start==0)for ( a=0; a<A->nseq; a++)TGP[a][0]=-1;
5388 else for ( a=0; a<A->nseq; a++)TGP[a][0]=(is_gap(A->seq_al[a][start-1])==1)?0:1;
5390 if ( end==A->len_aln-1)for ( a=0; a<A->nseq; a++)TGP[a][1]=-1;
5391 else for ( a=0; a<A->nseq; a++)TGP[a][1]=(is_gap(A->seq_al[a][start-1])==1)?0:1;
5395 /* 0: in the middle of the alignement
5397 2: q left gap is the continuation of another gap that was open outside the bloc ( don't open it)
5400 for ( a=0; a< A->nseq; a++)
5404 for ( b=0; b< start; b++)
5405 if ( !is_gap(A->seq_al[a][b]))
5409 if (is_gap(A->seq_al[a][start-1]) && TGP[a][0]!=1)
5412 for ( b=(start-1); b>=0 && continue_loop==1; b--)
5413 {TGP[a][0]-= ( is_gap(A->seq_al[a][b])==1)?1:0;
5414 continue_loop= (is_gap(A->seq_al[a][b])==1)?continue_loop:0;
5418 else if (is_gap(A->seq_al[a][start-1]) && TGP[a][0]==1)
5422 for ( b=(start-1); b>=0 && continue_loop==1; b--)
5423 {TGP[a][0]+= ( is_gap(A->seq_al[a][b])==1)?1:0;
5424 continue_loop= (is_gap(A->seq_al[a][b])==1)?continue_loop:0;
5427 for ( b=(A->len_aln-1); b>end; b--)
5428 if ( !is_gap(A->seq_al[a][b]))
5433 int gap_type ( char a, char b)
5435 /*gap_type -/-:0, X/X:-1 X/-:1, -/STAR:2*/
5437 if ( is_gap(a) && is_gap(b))
5439 else if ( !is_gap(a) && !is_gap(b))
5441 else if ( !is_gap(a))
5443 else if ( !is_gap(b))
5449 /******************************COPYRIGHT NOTICE*******************************/
5450 /*© Centro de Regulacio Genomica */
5452 /*Cedric Notredame */
5453 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
5454 /*All rights reserved.*/
5455 /*This file is part of T-COFFEE.*/
5457 /* T-COFFEE is free software; you can redistribute it and/or modify*/
5458 /* it under the terms of the GNU General Public License as published by*/
5459 /* the Free Software Foundation; either version 2 of the License, or*/
5460 /* (at your option) any later version.*/
5462 /* T-COFFEE is distributed in the hope that it will be useful,*/
5463 /* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
5464 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
5465 /* GNU General Public License for more details.*/
5467 /* You should have received a copy of the GNU General Public License*/
5468 /* along with Foobar; if not, write to the Free Software*/
5469 /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
5470 /*............................................... |*/
5471 /* If you need some more information*/
5472 /* cedric.notredame@europe.com*/
5473 /*............................................... |*/
5477 /******************************COPYRIGHT NOTICE*******************************/