Fix core WST file
[jabaws.git] / binaries / src / tcoffee / t_coffee_source / aln_compare.c
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <stdarg.h>
5 #include <string.h>
6 #include <ctype.h>
7
8 #include "io_lib_header.h"
9 #include "util_lib_header.h"
10 #include "dp_lib_header.h"
11 #include "define_header.h"
12
13   
14 int aln_compare ( int argc, char *argv[])
15     {
16     int a, b, c,f;
17         
18         int ik =0;
19     Alignment *A, *B;
20     Sequence *SA, *SB, *TOT_SEQ=NULL;
21     Sequence *defined_residueA;
22     Sequence *defined_residueB;
23     char **seq_list;
24     int n_seq_file;
25     
26     
27     Sequence  *S=NULL;
28     Structure *ST=NULL;
29     Result *R=NULL;
30     char *buf1;
31     char *buf2;
32
33 /*PARAMETERS*/
34     char ***grep_list;
35     int n_greps;
36     char compare_mode[STRING];
37     char sim_aln[STRING];
38     char *alignment1_file;
39     char *alignment2_file;
40     
41     char *pep1_file;
42     char *pep2_file;
43     
44     char io_format[STRING];
45
46     int n_structure;
47     char **struct_file;
48     char **struct_format;
49     int *n_symbol;
50     char ***symbol_list;
51     int pep_compare;
52     int aln_compare;
53     int count;
54     int output_aln;
55     int output_aln_threshold;
56     char output_aln_file[LONG_STRING];
57     char output_aln_format[LONG_STRING];
58     char output_aln_modif[LONG_STRING];
59 /*LIST VARIABLES*/
60     Constraint_list *CL_A;
61     Constraint_list *CL_B;
62
63     int pos_in_clist;
64     Sequence *CLS;
65
66
67 /*Column Comparison Variables*/
68  int **posA;
69  int **posB;
70  int **seq_cache;
71  int is_same;
72  int *correct_column;
73 /*RESULTS_VARIABLES*/
74     int **tot_count;
75     int **pos_count;
76     int ***pw_tot_count;
77     int ***pw_pos_count;
78     int *glob;
79     int **pw_glob;
80     int s1, r1, s2, r2;
81 /*IO VARIABLES*/
82     int n_categories;
83     char ***category;
84     char *category_list;
85     int *n_sub_categories;
86     char sep_l;
87     char sep_r;
88     char io_file[STRING];
89     FILE *fp;
90     int **aln_output_count;
91     int **aln_output_tot;
92     
93 /*Sims VARIABLES*/
94     float **sim;
95     float **sim_param;
96     char sim_matrix[STRING];
97     
98     int sim_n_categories;
99     char ***sim_category;
100     char *sim_category_list;
101     int *sim_n_sub_categories;
102     
103     
104     if ( argc==1|| strm6 ( argv[1], "h", "-h", "help", "-help", "-man", "?"))
105       {
106         output_informations();
107       }
108     
109     argv=standard_initialisation (argv, &argc);
110 /*Declarations and Initializations*/
111     alignment1_file=vcalloc ( LONG_STRING, sizeof (char));
112     alignment2_file=vcalloc ( LONG_STRING, sizeof (char));
113     
114     pep1_file=vcalloc ( LONG_STRING, sizeof (char));
115     pep2_file=vcalloc ( LONG_STRING, sizeof (char));
116     
117
118     sprintf (compare_mode, "sp");
119     count=0;
120     pep_compare=0;
121     aln_compare=1;
122
123     
124     grep_list=vcalloc ( STRING, sizeof (char**));
125     for ( a=0; a< STRING; a++)grep_list[a]=declare_char (3, STRING);
126     n_greps=0;
127
128     n_structure=0;
129     struct_file=declare_char ( MAX_N_STRUC, LONG_STRING);
130     struct_format=declare_char (MAX_N_STRUC, STRING);
131     
132     n_symbol=vcalloc ( MAX_N_STRUC, sizeof (int));
133     symbol_list=vcalloc (MAX_N_STRUC, sizeof (char**));
134     for ( a=0; a< MAX_N_STRUC; a++)symbol_list[a]=declare_char ( 100, 100);
135     
136
137    
138     
139     n_categories=1;
140     category=vcalloc ( MAX_N_CATEGORIES, sizeof (char**));
141     for ( a=0; a< MAX_N_CATEGORIES; a++)category[a]=declare_char(100, STRING);
142     n_sub_categories=vcalloc ( 100, sizeof (int));
143     category_list=vcalloc ( LONG_STRING, sizeof (char));
144     sprintf ( category_list, "[*][*]=[ALL]");
145     
146
147     sim_n_categories=1;
148     sim_category=vcalloc ( MAX_N_CATEGORIES, sizeof (char**));
149     for ( a=0; a< MAX_N_CATEGORIES; a++)sim_category[a]=declare_char(100, STRING);
150     sim_n_sub_categories=vcalloc ( 100, sizeof (int));
151     sim_category_list=vcalloc ( LONG_STRING, sizeof (char));
152     sprintf ( sim_category_list, "[*][*]=[ALL]");
153     sprintf ( sim_aln, "al1");
154     sim_matrix[0]='\0';
155     sprintf ( sim_category_list, "[*][*]=[ALL]");
156     
157     sprintf ( io_format, "ht");
158     sprintf ( io_file, "stdout");
159     sep_l='[';
160     sep_r=']';
161
162     output_aln=0;
163     output_aln_threshold=100;
164     sprintf ( output_aln_file, "stdout");
165     sprintf ( output_aln_format, "clustalw");
166     sprintf ( output_aln_modif, "lower");
167 /*END OF INITIALIZATION*/
168     
169    
170
171     
172 /*PARAMETERS INPUT*/
173
174     for ( a=1; a< argc; a++)
175         {
176         if      (strcmp ( argv[a], "-f")==0)
177                 {
178                 sprintf  (io_file,"%s", argv[++a]);
179                 }
180         else if ( strcmp ( argv[a], "-sim_aln")==0)
181                 {
182                 sprintf  (sim_aln,"%s", argv[++a]);
183                 }
184         else if ( strcmp ( argv[a], "-sim_matrix")==0)
185                 {
186                 sprintf  (sim_matrix,"%s", argv[++a]);
187                 }
188         else if ( strm ( argv[a], "-compare_mode"))
189                 {
190                 sprintf ( compare_mode, "%s", argv[++a]);
191                 }
192         else if ( strcmp ( argv[a], "-sim_cat")==0)
193                 {       
194                 if ( argv[++a][0]!='[')
195                    {
196                    if ( strcmp ( argv[a], "3d_ali")==0)sprintf ( sim_category_list, "[b][b]+[h][h]=[struc]");
197                    else
198                        {
199                        fprintf ( stderr, "\n%s: Unknown category for distance measure", argv[a]);
200                        }
201                        
202                    }
203                 else
204                    {
205                    sprintf ( sim_category_list, "%s", argv[a]);
206                    }
207                 }
208         else if ( strcmp ( argv[a], "-grep_value")==0)
209                 {
210                 sprintf ( grep_list[n_greps][0], "%s", argv[++a]);
211                 sprintf ( grep_list[n_greps][1], "%s", argv[++a]);
212                 n_greps++;
213              
214                 }
215         else if ( strcmp ( argv[a], "-al1")==0)
216                 {
217                 
218                 sprintf ( alignment1_file, "%s", argv[++a]);
219                 
220                 }       
221         else if ( strcmp ( argv[a], "-al2")==0)
222                 {
223                 
224                 sprintf ( alignment2_file, "%s", argv[++a]);
225                 
226                 }
227         else if (  strcmp ( argv[a], "-pep1")==0)
228                 {
229                 pep_compare=1;
230                 sprintf ( pep1_file, "%s", argv[++a]);
231                 }
232         else if (  strcmp ( argv[a], "-pep2")==0)
233                 {
234                 pep_compare=1;
235                 sprintf ( pep2_file, "%s", argv[++a]);
236                 }
237         else if (  strcmp ( argv[a], "-pep")==0)
238                 {
239                 pep_compare=1;
240                 }
241         else if (  strcmp ( argv[a], "-count")==0)
242                 {
243                 count=1;
244                 }
245         else if (  strcmp ( argv[a], "-output_aln")==0)
246                 {
247                 output_aln=1;
248                 }
249         else if (  strcmp ( argv[a], "-output_aln_threshold")==0)
250                 {
251                 output_aln_threshold=atoi(argv[++a]);
252                 }
253         else if (  strcmp ( argv[a], "-output_aln_file")==0)
254                 {
255                 sprintf (output_aln_file,"%s",argv[++a]);
256                 }
257         else if (  strcmp ( argv[a], "-output_aln_format")==0)
258                 {
259                 sprintf (output_aln_format,"%s",argv[++a]);
260                 }
261         else if (  strcmp ( argv[a], "-output_aln_modif")==0)
262                 {
263                 sprintf (output_aln_modif,"%s",argv[++a]);
264                 }
265         else if ( strcmp ( argv[a], "-st")==0)
266                 {
267                 sprintf ( struct_file  [n_structure], "%s", argv[++a]);
268                 if (!NEXT_ARG_IS_FLAG && is_a_struc_format (argv[a+1]))
269                     sprintf ( struct_format[n_structure], "%s", argv[++a]);
270                 else
271                     sprintf ( struct_format[n_structure], "%s", "pep");
272                 
273                 if ( !NEXT_ARG_IS_FLAG && strcmp ( argv[a+1], "conv")==0)
274                    {
275                    a++;
276                    while(!NEXT_ARG_IS_FLAG)
277                          { 
278                          sprintf ( symbol_list[n_structure][n_symbol[n_structure]], "%s", argv[++a]);
279                          n_symbol[n_structure]++; 
280                          }
281                    }
282                 else if (!NEXT_ARG_IS_FLAG)
283                          {
284                          symbol_list[n_structure]=make_symbols ( argv[++a], &n_symbol[n_structure]);
285                          }
286
287                 else 
288                          {
289                          symbol_list[n_structure]=make_symbols ( "any", &n_symbol[n_structure]);
290                          }
291                 
292                 n_structure++;
293                 
294                 }
295         else if ( strcmp  (argv[a], "-sep")==0)
296                 {
297                 if ( !NEXT_ARG_IS_FLAG)
298                     get_separating_char ( argv[++a][0], &sep_l, &sep_r);
299                 else
300                     sep_l=sep_r=' ';
301                 }
302         else if ( strncmp ( argv[a], "-io_format",5)==0)
303                 {
304                 sprintf ( io_format, "%s", argv[++a]);
305                 }
306         else if ( strcmp ( argv[a], "-io_cat")==0)
307                 {
308                 if ( argv[++a][0]!='[')
309                    {
310                    if ( strcmp ( argv[a], "3d_ali")==0)sprintf ( category_list, "[b][b]+[h][h]=[struc];[*][*]=[tot]");
311                    }
312                 else
313                    {
314                    sprintf ( category_list, "%s", argv[a]);
315                    }
316                 }
317         else
318                 {
319                 fprintf ( stdout, "\nOPTION %s UNKNOWN[FATAL]\n", argv[a]);
320                 myexit (EXIT_FAILURE);
321                 }
322         }
323  
324 /*PARAMETER PROCESSING*/
325
326 if ( pep_compare==1 || count==1)aln_compare=0;
327 if ( aln_compare==1)pep_compare=0; 
328
329 /*READ THE TOTAL SEQUENCES*/
330        seq_list=declare_char ( 100,STRING);
331        n_seq_file=0;
332
333         if (  alignment1_file[0] && !check_file_exists ( alignment1_file))
334          {
335            fprintf (stderr, "\nERROR: %s DOES NOT EXIST[FATAL:%s]\n",  alignment1_file, PROGRAM);
336            myexit(EXIT_FAILURE);
337          } 
338         if (  alignment2_file[0] && !check_file_exists ( alignment2_file))
339          {
340            fprintf (stderr, "\nERROR: %s DOES NOT EXIST[FATAL:%s]\n",  alignment2_file, PROGRAM);
341            myexit(EXIT_FAILURE);
342          }
343         if (  pep1_file[0] && !check_file_exists ( pep1_file))
344          {
345            fprintf (stderr, "\nERROR: %s DOES NOT EXIST[FATAL:%s]\n",  pep1_file, PROGRAM);
346            myexit(EXIT_FAILURE);
347          }
348         if (  pep2_file[0] && !check_file_exists ( pep2_file))
349          {
350            fprintf (stderr, "\nERROR: %s DOES NOT EXIST[FATAL:%s]\n",  pep2_file, PROGRAM);
351            myexit(EXIT_FAILURE);
352          }
353         
354        if ( alignment1_file[0])sprintf ( seq_list[n_seq_file++], "A%s", alignment1_file);
355        if ( alignment2_file[0])sprintf ( seq_list[n_seq_file++], "A%s", alignment2_file);
356        if ( pep1_file[0])sprintf ( seq_list[n_seq_file++], "S%s", pep1_file);
357        if ( pep2_file[0])sprintf ( seq_list[n_seq_file++], "S%s", pep2_file);
358
359      
360
361        TOT_SEQ=read_seq_in_n_list ( seq_list, n_seq_file, NULL, NULL);
362              
363        A=declare_aln (TOT_SEQ);
364        B=declare_aln (TOT_SEQ);
365
366
367 /*1 COMPARISON OF THE SEQUENCES*/
368     if ( pep_compare==1 || count==1)
369        {
370        f=0; 
371
372        if ( pep1_file[0]!='\0')SA=main_read_seq (pep1_file);
373        else if (alignment1_file[0]!='\0')  
374            {           
375            main_read_aln ( alignment1_file, A);
376            SA=aln2seq ( A);
377            }
378        else 
379            {
380            main_read_aln ("stdin", A);
381            sprintf ( alignment1_file, "stdin");
382            SA=aln2seq ( A);
383            }
384        if ( pep2_file[0]!='\0')SB=main_read_seq (pep2_file);
385        else if (alignment2_file[0]!='\0') 
386            {
387            main_read_aln ( alignment2_file, B);
388            
389            SB=aln2seq (B);
390            }
391        else 
392            {
393            SB=SA;
394            sprintf ( alignment2_file, "%s", alignment1_file );
395            }
396        buf1=(pep1_file[0]!='\0')?pep1_file: alignment1_file;
397        buf2=(pep2_file[0]!='\0')?pep2_file: alignment2_file;
398        /*Output of the Results*/    
399
400        fp=vfopen ( io_file, "w");      
401        
402
403        if ( count==1)
404           {
405           fprintf (fp, "Number of seq: %d %d\n", SA->nseq,SA->nseq);
406           for ( a=0; a< SA->nseq; a++) fprintf (fp, "%-15s %d %d\n", SA->name[a], (int)strlen (SA->seq[a]), (int)strlen (SA->seq[a]));
407           }
408               
409        if (SA->nseq!=SB->nseq)
410            {
411            
412            fprintf ( fp, "DIFFERENCE TYPE 1: Different number of sequences %3d/%3d",SA->nseq,SB->nseq);
413            f=1;
414            }
415
416        trim_seq ( SA, SB);
417        for ( a=0; a< SA->nseq; a++)
418            {
419            lower_string (SA->seq[a]);
420            lower_string (SB->seq[a]);
421            ungap ( SA->seq[a]);
422            ungap ( SB->seq[a]);
423
424            if ( strcmp ( SA->seq[a], SB->seq[a])!=0) 
425                {
426                fprintf ( fp, "DIFFERENCE TYPE 2: %s is different in the 2 files\n", SA->name[a]);
427                f=1;
428                }
429            }
430        for ( a=0; a< SA->nseq; a++)
431            {
432            lower_string (SA->seq[a]);
433            lower_string (SB->seq[a]);
434            ungap ( SA->seq[a]);
435            ungap ( SB->seq[a]);
436
437            if ( strlen ( SA->seq[a])!= strlen (SB->seq[a]))
438                {
439                  fprintf ( fp, "DIFFERENCE TYPE 3: %s has != length in the 2 files (%d-%d)\n", SA->name[a],(int)strlen ( SA->seq[a]), (int)strlen (SB->seq[a]));
440                f=1;
441                }
442            }
443        if ( f==1)
444            {
445            fprintf ( fp, "\nDIFFERENCES found between:\n\t%s\n\t%s\n**********\n\n",buf1, buf2);
446            }
447        fclose (fp);
448        }
449
450 /*2 COMPARISON OF THE ALIGNMENTS*/                       
451     else if ( aln_compare==1)
452        {
453
454        n_categories=parse_category_list ( category_list, category, n_sub_categories);
455        sim_n_categories=parse_category_list ( sim_category_list, sim_category, sim_n_sub_categories);
456        
457        main_read_aln ( alignment1_file, A);
458        main_read_aln ( alignment2_file, B);
459        CLS=trim_aln_seq ( A, B);
460        
461        
462        defined_residueA=get_defined_residues (A);
463        defined_residueB=get_defined_residues (B);
464        
465
466        A=thread_defined_residues_on_aln(A, defined_residueA);
467        A=thread_defined_residues_on_aln(A, defined_residueB);
468        B=thread_defined_residues_on_aln(B, defined_residueA);
469        B=thread_defined_residues_on_aln(B, defined_residueB);
470       
471      
472        CL_A=declare_constraint_list ( CLS, NULL, NULL, 0, NULL, NULL);
473        CL_B=declare_constraint_list ( CLS, NULL, NULL, 0, NULL, NULL);
474        
475       
476        CL_A=aln2constraint_list (A,CL_A, "sim");
477        CL_B=aln2constraint_list (B,CL_B, "sim");
478        
479        
480        
481        glob=vcalloc ( A->nseq+1, sizeof (int));
482        pw_glob=declare_int ( A->nseq+1, A->nseq+1);
483
484
485        if ( strm( compare_mode, "sp"))
486             {
487               int *entry;
488               while (entry=extract_entry (CL_A))
489                 {
490                   s1=entry[SEQ1];
491                   s2=entry[SEQ2];
492                   glob[0]++;
493                   glob[s1+1]++;
494                   glob[s2+1]++;
495                   pw_glob[s1][s2]++;
496                   pw_glob[s2][s1]++;
497                   entry[MISC]=1;
498                   
499                   if ((main_search_in_list_constraint (entry,&pos_in_clist,4,CL_B))!=NULL)
500                     add_entry2list (entry,CL_A);
501                   
502                 }
503             }
504        else if ( strm( compare_mode, "column") )
505             {
506               int *entry;
507               posA=aln2pos_simple_2(A);
508               posB=aln2pos_simple_2(B);
509               seq_cache=declare_int ( A->nseq, A->len_aln+1);
510               for ( a=0; a< A->len_aln; a++)
511                 for ( b=0; b<B->len_aln; b++)
512                   {
513                     is_same=compare_pos_column(posA, a, posB, b, A->nseq);    
514                     if(is_same)
515                     {
516                       for (c=0; c< A->nseq;c++)
517                           {
518                             if(posA[c][a]>0)
519                             {
520                                 seq_cache[c][posA[c][a]]=1;
521                             }
522                           }
523                     }
524                   }
525               
526               while (entry=extract_entry (CL_A))
527                 {
528                   s1=entry[SEQ1];
529                   s2=entry[SEQ2];
530                   r1=entry[R1];
531                   entry[MISC]=0;
532                   glob[0]++;
533                   glob[s1+1]++;
534                   glob[s2+1]++;
535                   pw_glob[s1][s2]++;
536                   pw_glob[s2][s1]++;
537                   
538                   if (seq_cache[s1][r1])
539                   {
540                      entry[MISC]=1;
541                      add_entry2list (entry, CL_A);
542                   }
543                 }
544               free_int (posA, -1);
545               free_int (posB, -1);
546               free_int (seq_cache, -1);      
547             }   
548         else if ( strm( compare_mode, "tc") )
549         {
550               correct_column = vcalloc(A->len_aln+1, sizeof (int));
551               posA=aln2pos_simple_2(A);
552               posB=aln2pos_simple_2(B);
553               for ( a=0; a< A->len_aln; a++)
554                 for ( b=0; b<B->len_aln; b++)
555                   {
556                     is_same = compare_pos_column(posA, a, posB, b, A->nseq);
557                     if(is_same)
558                       correct_column[a] = is_same;
559                   }
560               free_int (posB, -1);
561         }
562         
563        for ( a=0; a< n_structure; a++)
564            {
565                ST=read_structure (struct_file[a],struct_format[a], A,B,ST,n_symbol[a], symbol_list[a]); 
566            }
567        
568       /*RESULT ARRAY DECLARATION*/    
569     
570        tot_count=declare_int (n_categories+1, A->nseq+1);
571        pos_count=declare_int (n_categories+1, A->nseq+1);
572        pw_tot_count=vcalloc ( A->nseq, sizeof (int**));
573        for ( a=0; a< A->nseq; a++)pw_tot_count[a]=declare_int ( A->nseq, n_categories);
574
575     
576        pw_pos_count=vcalloc ( A->nseq, sizeof (int**));
577        for ( a=0; a< A->nseq; a++)pw_pos_count[a]=declare_int ( A->nseq, n_categories);
578  
579      /*COMPARISON MODULE*/
580         if (strm( compare_mode, "tc") )
581         {
582           int column_is_structure;
583           int pair_is_structure;
584
585           for ( a=0; a< n_categories; a++)
586            {     
587              for(b = 0; b < A->len_aln; b++)
588              {
589                column_is_structure = 1;
590                pair_is_structure = 0;
591                for (s1=0; s1< A->nseq-1; s1++)
592                {
593                  for(s2=s1+1; s2 < A->nseq; s2++)
594                  {
595                    if ((posA[s1][b]>0) && (posA[s2][b]>0))
596                    {  
597                      pair_is_structure=is_in_struct_category(s1,s2,posA[s1][b],posA[s2][b],ST,category[a], n_sub_categories[a]);     
598                      column_is_structure = (column_is_structure && pair_is_structure);
599                      if(pair_is_structure)
600                      {
601                         tot_count[a][s1+1]++;
602                         tot_count[a][s2+1]++;
603                         pw_tot_count[s1][s2][a]++;
604                         pw_tot_count[s2][s1][a]++;
605                         if (correct_column[b])
606                           {
607                             pw_pos_count[s1][s2][a]++;
608                             pw_pos_count[s2][s1][a]++;
609                             pos_count[a][s1+1]++;
610                             pos_count[a][s2+1]++;
611                           }
612                      }
613                    }
614                  }
615                }               
616                 if(column_is_structure && pair_is_structure)
617                 {
618                  tot_count[a][0]++;
619                  if(correct_column[b]) 
620                  {
621                    pos_count[a][0]++;
622                  }
623                 }
624              }
625            }
626            free_int (posA, -1);
627         }
628         else
629         {
630         /*COMPARISON MODULE*/
631           int *entry;
632           for ( a=0; a< n_categories; a++)
633               {
634                 while (entry=extract_entry (CL_A))
635                   {
636                     s1=entry[SEQ1];
637                     s2=entry[SEQ2];
638                     r1=entry[R1];
639                     r2=entry[R2];
640                     c= entry[MISC];
641                     if ( is_in_struct_category ( s1, s2, r1, r2, ST, category[a], n_sub_categories[a]))
642                       {
643                             tot_count[a][0]++;
644                             tot_count[a][s1+1]++;
645                             tot_count[a][s2+1]++;
646                             pw_tot_count[s1][s2][a]++;
647                             pw_tot_count[s2][s1][a]++;
648                             if ( c==1)
649                               {
650                                 pw_pos_count[s1][s2][a]++;
651                                 pw_pos_count[s2][s1][a]++;
652                                 pos_count[a][0]++;
653                                 pos_count[a][s1+1]++;
654                                 pos_count[a][s2+1]++;
655                               }
656                         }
657                       }
658                   }
659         }
660 //      printf("pos=%d,tot=%d\n", pos_count[0][0], tot_count[0][0]);
661     /*Measure of Aligned Sequences Similarity*/
662     
663        sim=get_aln_compare_sim ((strcmp (sim_aln, "al1")==0)?A:B, ST,sim_category[0], sim_n_sub_categories[0], sim_matrix);
664        sim_param=analyse_sim ((strcmp (sim_aln, "al1")==0)?A:B, sim);
665
666    /*Fill the Result_structure*/
667        R=vcalloc ( 1, sizeof (Result));
668     
669        R->grep_list=grep_list;
670        R->n_greps=n_greps;
671        R->A=A;
672        R->B=B;
673
674        R->S=S;
675        R->ST=ST;
676        R->sim_aln=sim_aln;
677        R->alignment1_file=alignment1_file;
678        R->alignment2_file=alignment2_file;
679        R->io_format=io_format;
680        R->n_structure=n_structure;
681        R->struct_file=struct_file;
682        R->struct_format=struct_format;
683        R->n_symbol=n_symbol;
684        R->symbol_list=symbol_list;
685
686        
687        R->tot_count=tot_count;
688        R->pos_count=pos_count;
689        R->pw_tot_count=pw_tot_count;
690        R->pw_pos_count=pw_pos_count;
691        R->glob=glob;
692        R->pw_glob=pw_glob;
693        R->n_categories=n_categories;
694        R->category=category;
695        R->category_list=category_list;
696        R->n_sub_categories=n_sub_categories;
697        R->sim=sim;
698        R->sim_param=sim_param;
699        R->sim_matrix=sim_matrix;
700        R->sim_n_categories=sim_n_categories;
701        R->sim_category=sim_category;
702        R->sim_category_list=sim_category_list;
703        R->sim_n_sub_categories=sim_n_sub_categories; 
704        R->sep_r=sep_r;
705        R->sep_l=sep_l;
706
707       /*Output of the Results*/    
708
709        fp=vfopen ( io_file, "w");      
710        fp=output_format (io_format, fp, R);
711        vfclose ( fp); 
712
713
714      /*Rewriting of Alignment A*/
715        if ( output_aln)
716          {
717            int *entry;
718            A->residue_case=2;
719            aln_output_tot  =declare_int ( A->nseq, A->len_aln+1);
720            aln_output_count=declare_int ( A->nseq, A->len_aln+1);
721            
722            while (entry=extract_entry (CL_A))
723              {
724                aln_output_tot[entry[SEQ1]][entry[R1]]++;
725                aln_output_tot[entry[SEQ2]][entry[R2]]++;
726                
727                aln_output_count[entry[SEQ1]][entry[R1]]+=entry[MISC];
728                aln_output_count[entry[SEQ2]][entry[R2]]+=entry[MISC];
729              }
730           
731            for ( a=0; a< A->nseq; a++)
732              {
733                
734                for (c=0, b=0; b< A->len_aln; b++)
735                  {
736                    if ( !is_gap(A->seq_al[a][b]))
737                      {
738                        c++;
739                        if ( aln_output_tot[a][c] && ((aln_output_count[a][c]*100)/aln_output_tot[a][c])<output_aln_threshold)
740                          { 
741                            if ( strm (output_aln_modif, "lower"))
742                                 A->seq_al[a][b]=tolower(A->seq_al[a][b]);
743                            else
744                              A->seq_al[a][b]=output_aln_modif[0];
745                          }
746                        else A->seq_al[a][b]=toupper(A->seq_al[a][b]);
747                      }
748                      
749                  }
750              }
751            A->score_aln=(int)(R->tot_count[0][0]==0)?0:((R->pos_count[0][0]*100)/(float)R->tot_count[0][0]);
752
753            output_format_aln (output_aln_format,A,NULL,output_aln_file);
754            
755            free_int ( aln_output_tot, -1);
756            free_int ( aln_output_count, -1  );
757          }
758        }
759     return EXIT_SUCCESS;
760     }   
761 /************************************************************************************/
762 /*                                                                                  */
763 /*                            OUTPUT                                                */
764 /*                                                                                  */
765 /*                                                                                  */
766 /************************************************************************************/
767 FILE *output_format (char *iof,FILE *fp, Result *R)
768       {
769       int a;
770       int l;
771       
772       /*
773       H: files Header
774       h: basic header;
775       s: sequence results
776       t: total results (global);
777       p: pairwise_results;
778       */
779       l=strlen ( iof);
780
781       for ( a=0; a< l; a++)
782           {
783           if ( iof[a]=='H')fp=output_large_header (fp,R);
784           else if ( iof[a]=='h')fp=output_header (fp,R);
785           else if ( iof[a]=='t')fp=output_total_results (fp, R);
786           else if ( iof[a]=='s')fp=output_sequence_results (fp,R);
787           else if ( iof[a]=='p')fp=output_pair_wise_sequence_results(fp,R);
788           }
789       return fp;
790       }
791
792 FILE *output_pair_wise_sequence_results (FILE *fp,  Result *R)
793      {
794      int a,c,d;
795      
796      
797      for ( c=0; c<(R->A)->nseq-1; c++)
798          {
799          for ( d=c+1; d< (R->A)->nseq; d++)    
800              {
801              fprintf (fp, "%-10s %-10s%s",(R->A)->name[c],(R->A)->name[d],SSPACE);
802              fprintf (fp, "%5.1f%s", R->sim[c][d], SSPACE);
803              
804              for (a=0; a< R->n_categories; a++)
805                  {
806                  fprintf ( fp, "%5.1f ",(R->pw_tot_count[c][d][a]==0)?0:((float)(R->pw_pos_count[c][d][a]*100)/(float)R->pw_tot_count[c][d][a]));
807                  fprintf ( fp, "%c%5.1f%c%s",(R->sep_l),(R->pw_glob[c][d]==0)?0:((float)(R->pw_tot_count[c][d][a]*100)/(float)R->pw_glob[c][d]),(R->sep_r),SSPACE);
808                  }
809              fprintf ( fp, "%c%5d%c\n",(R->sep_l), R->pw_glob[c][d],(R->sep_r));
810              }
811          }
812          
813     return fp;
814     }
815 FILE *output_sequence_results (FILE *fp,  Result *R)
816      {
817      int a,c;
818     
819      for ( c=1; c<=R->A->nseq; c++)
820          {
821          fprintf (fp, "%-10s %-10s%s",(R->A)->name[c-1], "..",SSPACE);
822          fprintf (fp, "%5.1f%s", R->sim_param[c-1][0],SSPACE);
823          for (a=0; a< R->n_categories; a++)
824                  {
825                  fprintf ( fp, "%5.1f ",(R->tot_count[a][c]==0)?0:((float)(R->pos_count[a][c]*100)/(float)R->tot_count[a][c]));
826                  fprintf ( fp, "%c%5.1f%c%s",(R->sep_l),(R->glob[c]==0)?0:((float)(R->tot_count[a][c]*100)/(float)R->glob[c]),(R->sep_r), SSPACE);
827                  }
828          fprintf ( fp, "%c%5d%c\n",(R->sep_l), R->glob[c],(R->sep_r));
829          }
830     return fp;
831     }
832          
833
834 FILE *output_total_results (FILE *fp,  Result *R)
835      {
836      int a;
837      
838      
839      fprintf ( fp, "%-13s %-7d%s",extract_suffixe (R->alignment1_file),(R->A)->nseq, SSPACE);
840      fprintf (fp, "%5.1f%s", R->sim_param[(R->A)->nseq][0], SSPACE);
841      for (a=0; a< R->n_categories; a++)
842          {
843          fprintf ( fp, "%5.1f ",(R->tot_count[a][0]==0)?0:((float)(R->pos_count[a][0]*100)/(float)R->tot_count[a][0]));
844          fprintf ( fp, "%c%5.1f%c%s",(R->sep_l),(R->glob[0]==0)?0:((float)(R->tot_count[a][0]*100)/(float)R->glob[0]),(R->sep_r), SSPACE);
845          }
846      fprintf ( fp, "%c%5d%c\n",(R->sep_l),  R->glob[0],(R->sep_r));
847      return fp;
848      }
849 FILE *output_header (FILE *fp, Result *R)
850      {
851      int a;
852
853      
854      fprintf ( fp, "%s\n",generate_string ( R->n_categories*(13+strlen(SSPACE))+31+2*strlen(SSPACE),'*'));
855  
856      fprintf ( fp, "%-10s %-10s %s%-3s%s", "seq1", "seq2",SSPACE,"Sim",SSPACE);
857      
858      for ( a=0; a< R->n_categories; a++)
859          fprintf ( fp, "%-12s%s ",R->category[a][0], SSPACE);
860      fprintf (fp, "%-5s", "Tot");
861      fprintf (fp, "\n");
862      return fp;
863      }
864 FILE *output_large_header ( FILE *fp, Result *R)
865      {
866      int a, b;
867
868      fprintf ( fp, "AL1: %s\n", R->alignment1_file);
869      fprintf ( fp, "AL2: %s\n", R->alignment2_file);
870      for ( a=0; a< R->n_structure; a++)
871           {
872            fprintf (fp, "ST %d: %s [%s]/[", a, R->struct_file[a], R->struct_format[a]);
873            for ( b=0; b< R->n_symbol[a]; b++)fprintf (fp, "%s ", R->symbol_list[a][b]);
874            fprintf ( fp, "]\n");
875           }
876      return (fp);
877      }
878
879 void get_separating_char ( char s, char *l, char *r)
880     {
881     if ( s=='{' || s=='}')
882        {
883        l[0]='{';
884        r[0]='}';
885        return;
886        }
887     else if ( s==']' || s=='[')
888        {
889        l[0]='[';
890        r[0]=']';
891        return;
892        }
893     else if ( s==')' || s=='(')
894        {
895        l[0]='(';
896        r[0]=')';
897        return;
898        }
899     else
900        {
901        l[0]=s;
902        r[0]=s;
903        return;
904        }
905     }
906
907 /************************************************************************************/
908 /*                                                                                  */
909 /*                            SIM MEASURE                                           */
910 /*                                                                                  */
911 /*                                                                                  */
912 /************************************************************************************/
913 float **get_aln_compare_sim ( Alignment *A, Structure *ST, char **cat, int n_cat, char *matrix)
914       {
915       int a, b, c;
916
917       float **sim;
918       char cr1, cr2;
919       int r1, r2;
920       int p1, p2;
921       float pos,tot;
922
923       sim=declare_float ( A->nseq, A->nseq);
924
925       for ( a=0; a<A->nseq-1; a++)
926           {
927           for (b=a+1; b< A->nseq; b++)
928               {
929               for ( r1=0, r2=0,tot=0, pos=0,c=0; c< A->len_aln; c++)
930                   {
931                   p1=is_gap(A->seq_al[a][c]);
932                   p2=is_gap(A->seq_al[b][c]);
933                   
934                   r1+=1-p1;
935                   r2+=1-p2;
936                   cr1=A->seq_al[a][c];
937                   cr2=A->seq_al[b][c];
938                   if (!p1 && !p2)
939                       {
940                       if (is_in_struct_category (a, b, r1, r2, ST, cat, n_cat))
941                           {
942                           tot++;
943                           pos+=is_in_same_group_aa ( cr1, cr2, 0, NULL, matrix);
944                           }
945                       }
946                   }
947               sim[a][b]=sim[b][a]=(tot==0)?0:((pos*100)/tot);
948               }
949           }
950       return sim;
951       }
952 float **analyse_sim ( Alignment *A, float **sim)
953         {
954         int a,b,c;
955
956         float **an, d;
957         
958         
959         an=declare_float ( A->nseq+1,2);
960         
961         for (d=0, a=0; a< A->nseq;a++)
962                 {
963                 for ( b=0; b< A->nseq; b++)
964                         {
965                         if ( b!=a)
966                                 {
967                                 an[a][0]+=sim[a][b];
968                                 an[A->nseq][0]+=sim[a][b];
969                                 d++;
970                                 }
971                         }
972                 an[a][0]=((float)an[a][0]/(float)(A->nseq-1));
973                 }
974         
975         
976         an[A->nseq][0]=an[A->nseq][0]/d;
977                 
978         for ( d=0,a=0; a< A->nseq; a++)
979                 {
980                 for ( b=0; b< A->nseq; b++)
981                         {
982                         if ( b!=a)
983                                 {
984                                 c=an[a][0]-sim[a][b];
985                                 an[a][1]+=(c>0)?c:-c;
986                                 an[A->nseq][1]+=(c>0)?c:-c;
987                                 d++;
988                                 }
989                         }
990                 an[a][1]=((float)an[a][1]/(float)(A->nseq-1));
991                 }
992         
993         an[A->nseq][1]=an[A->nseq][1]/d;
994         
995         return an;
996         }
997                 
998
999
1000
1001
1002
1003
1004
1005 /************************************************************************************/
1006 /*                                                                                  */
1007 /*                            STRUC ANALYSE                                         */
1008 /*                                                                                  */
1009 /*                                                                                  */
1010 /************************************************************************************/
1011 int is_in_struct_category ( int s1, int s2, int r1, int r2, Structure *ST, char **cat, int n_cat)
1012     {
1013     int a;
1014     static char *struc_r1;
1015     static char *struc_r2;
1016     char first[STRING];
1017     char second[STRING];
1018     static int **r;
1019     
1020     
1021
1022     
1023     if ( ST==NULL)return 1;
1024     if ( struc_r1!=NULL)
1025        {
1026        vfree ( struc_r1);
1027        vfree ( struc_r2);
1028        }
1029
1030     if ( r==NULL)r=declare_int (2, 2);
1031     else r[0][0]=r[1][1]=r[1][0]=r[0][1]=0;
1032     
1033     
1034     struc_r1=get_structure_residue ( s1, r1, ST);
1035     struc_r2=get_structure_residue ( s2, r2, ST);
1036
1037     
1038     for ( a=1; a< n_cat; a+=2)
1039         {
1040         sprintf ( first, "%s", cat[a]);
1041         sprintf ( second,"%s", cat[a+1]);
1042         r[0][0]=struc_matches_pattern ( struc_r1, first);
1043         r[0][1]=struc_matches_pattern ( struc_r2, first);
1044         r[1][0]=struc_matches_pattern ( struc_r1, second);
1045         r[1][1]=struc_matches_pattern ( struc_r2, second);
1046         
1047         if ( (r[0][0]&&r[1][1])||(r[1][0]&&r[0][1]))return 1;
1048         }       
1049    return 0;
1050    }
1051
1052 char * get_structure_residue (int seq, int res, Structure *ST)
1053     {
1054     int a;
1055     char *s;
1056     
1057     s=vcalloc ( ST->n_fields+1, sizeof (char));
1058     for (a=0; a< ST->n_fields; a++)
1059         s[a]=ST->struc[seq][res][a];
1060     s[a]='\0';
1061     return s;
1062     }
1063
1064 int struc_matches_pattern ( char *struc, char *pattern)
1065     {
1066     char p[STRING];
1067     char *y;
1068     
1069     int b,l;
1070     
1071
1072     sprintf ( p, "%s", pattern);
1073     
1074     if ( strcmp (p, "*")==0)return 1;
1075     else
1076         {
1077         l=strlen ( struc);
1078         y=strtok ( p, ".");
1079
1080         for ( b=0; b<l; b++)
1081                 {
1082
1083                 if ( strcmp ( y, "*")==0);                 
1084                 else if ( strchr(y, struc[b])==NULL)return 0;
1085                 y=strtok ( NULL, ".");
1086                 
1087                 if ( y==NULL && b<(l-1))
1088                    {
1089                    crash ( "\nA Structure is Missing[FATAL]");
1090                    }
1091                 }
1092         }
1093
1094     return 1;
1095     }
1096     
1097
1098 Structure * read_structure (char *fname, char *format, Alignment *A,Alignment *B, Structure *ST, int n_symbols, char **symbol_table)        
1099         {
1100                
1101         Alignment *StrucAln;
1102         Sequence *SEQ=NULL, *SA=NULL;
1103         
1104         
1105
1106
1107         if (ST==NULL)ST=declare_structure (A->nseq, A->seq_al);
1108         else
1109             ST=extend_structure ( ST);
1110         
1111         StrucAln=declare_Alignment(NULL); 
1112         if (strm ( format, "pep"))
1113                         {
1114                         SEQ=main_read_seq ( fname);
1115                         StrucAln=seq2aln(SEQ,StrucAln,0);
1116                         }
1117         
1118         else if ( strcmp ( format, "aln")==0)
1119                         {
1120                         StrucAln=main_read_aln (fname, StrucAln);
1121                         }
1122
1123         reorder_aln(StrucAln, A->name, A->nseq);
1124
1125         SA=aln2seq(StrucAln);
1126         string_array_convert (SA->seq, SA->nseq, n_symbols, symbol_table);
1127         seq2struc (SA, ST);
1128         
1129         free_aln(StrucAln);
1130         if (SEQ)free_sequence (SEQ, SEQ->nseq);
1131         
1132         return ST;
1133         }
1134
1135 int parse_category_list ( char *category_list_in, char ***category, int*n_sub_categories)
1136     {
1137     int n,a;
1138     char *category_list;
1139     char *y,*z;
1140     category_list=vcalloc ( strlen(category_list_in)+1, sizeof (char));
1141     sprintf (category_list, "%s", category_list_in);
1142
1143     n=0;
1144     z=category_list;
1145     while ((y=strtok(z, ";"))!=NULL)
1146         {
1147         sprintf ( category[n++][2], "%s", y);
1148         z=NULL;
1149         }
1150
1151     
1152     for ( a=0; a< n; a++)
1153         {
1154          sprintf (category_list,"%s",strtok(category[a][2], "="));
1155          sprintf (category[a][0],"%s",strtok (NULL, "="));
1156          sprintf ( category[a][++n_sub_categories[a]],"%s", strtok(category_list, "[]+"));
1157          while ( ((y=strtok(NULL, "[]+"))!=NULL))
1158                 {
1159                 if ( strcmp (y, "#")==0)y=category[a][n_sub_categories[a]];
1160                 sprintf ( category[a][++n_sub_categories[a]],"%s",y);
1161                 }
1162         }
1163     return n;
1164     }
1165
1166
1167 int is_a_struc_format (char *format)
1168     {
1169     if ( strcmp ( format, "pep")==0)return 1;
1170     if ( strcmp ( format, "aln")==0)return 1;
1171     return 0;
1172     }
1173 /************************************************************************************/
1174 /*                                                                                  */
1175 /*                            Informations                                          */
1176 /*                                                                                  */
1177 /*                                                                                  */
1178 /************************************************************************************/
1179 void output_informations ()
1180 {
1181 fprintf ( stderr, "\nPROGRAM: %s (%s)\n",PROGRAM,VERSION);
1182 fprintf ( stderr, "******INPUT***************************");
1183 fprintf ( stderr, "\n-al1 al1_file");
1184 fprintf ( stderr, "\n-al2 al2_file");
1185 fprintf ( stderr, "\n-compare_mode [sp] or column");
1186 fprintf ( stderr, "\n-pep (compare only the sequences");
1187 fprintf ( stderr, "\n-count");
1188 fprintf ( stderr, "\n-pep1 pep1_file");
1189 fprintf ( stderr, "\n-pep1 pep2_file");
1190 fprintf ( stderr, "\n-st  str_file st_format conversion");
1191 fprintf ( stderr, "\n   **st_format:  aln, pep");
1192 fprintf ( stderr, "\n   **conversion: 3d_ali,  conv abcZ #X");
1193 fprintf ( stderr, "\nNOTE: Several structures in a row are possible"); 
1194
1195 fprintf ( stderr, "\n\n****DISTANCE MEASURE*****************");
1196 fprintf ( stderr, "\n-sim_cat category_format  or category_name");
1197 fprintf ( stderr, "\n       **category_format: [*][*]=[tot]");
1198 fprintf ( stderr, "\n       **category_name  : 3d_ali (<=>[h][h]+[e][e]=[Struc]");
1199 fprintf ( stderr, "\n-sim_matrix matrix_name");
1200 fprintf ( stderr, "\n          **matrix_name: idmat,pam250mt..");
1201 fprintf ( stderr, "\n-sim_aln al1 or al2");
1202 fprintf ( stderr, "\n\n****COMPARISON***********************");
1203 fprintf ( stderr, "\n-io_cat category_format or category_name");
1204 fprintf ( stderr, "\n      **category_format: [*][*]=[tot]");
1205 fprintf ( stderr, "\n      **category_name  : 3d_ali(<=>[h][h]+[e][e]=[Struc];[*][*]=[Tot]");
1206 fprintf ( stderr, "\n\nNOTE: if two structures:");
1207 fprintf ( stderr, "\n[he.123][#]=[he123VShe123];[beh.*][he.2345]=[other]");
1208 fprintf ( stderr, "\n\n****OUTPUT****************************");
1209 fprintf ( stderr, "\n-f stdout");
1210 fprintf ( stderr, "\n   stderr");
1211 fprintf ( stderr, "\n   file_name");
1212 fprintf ( stderr, "\n-io_format hts");
1213 fprintf ( stderr, "\n           H ->large Header");
1214 fprintf ( stderr, "\n           h ->small Header");
1215 fprintf ( stderr, "\n           t->global (average)results");
1216 fprintf ( stderr, "\n           s ->average results for each sequence");
1217 fprintf ( stderr, "\n           p ->results for each pair of sequences");
1218 fprintf ( stderr, "\n-output_aln           Outputs al1 with conserved bits in Upper");
1219 fprintf ( stderr, "\n-output_aln_threshold [100]");
1220 fprintf ( stderr, "\n-output_aln_file      [stdout]");
1221 fprintf ( stderr, "\n-output_aln_format    [clustalw]");
1222 fprintf ( stderr, "\n-output_aln_modif     [lower]");
1223 fprintf ( stderr, "\n");
1224 myexit (EXIT_SUCCESS);
1225 }
1226
1227 /*********************************COPYRIGHT NOTICE**********************************/
1228 /*� Centro de Regulacio Genomica */
1229 /*and */
1230 /*Cedric Notredame */
1231 /*Tue Jul  1 10:00:54 WEST 2008. */
1232 /*All rights reserved.*/
1233 /*This file is part of T-COFFEE.*/
1234 /**/
1235 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1236 /*    it under the terms of the GNU General Public License as published by*/
1237 /*    the Free Software Foundation; either version 2 of the License, or*/
1238 /*    (at your option) any later version.*/
1239 /**/
1240 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1241 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1242 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1243 /*    GNU General Public License for more details.*/
1244 /**/
1245 /*    You should have received a copy of the GNU General Public License*/
1246 /*    along with Foobar; if not, write to the Free Software*/
1247 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1248 /*...............................................                                                                                      |*/
1249 /*  If you need some more information*/
1250 /*  cedric.notredame@europe.com*/
1251 /*...............................................                                                                                                                                     |*/
1252 /**/
1253 /**/
1254 /*      */
1255 /*********************************COPYRIGHT NOTICE**********************************/
1256 /******************************COPYRIGHT NOTICE*******************************/
1257 /*© Centro de Regulacio Genomica */
1258 /*and */
1259 /*Cedric Notredame */
1260 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
1261 /*All rights reserved.*/
1262 /*This file is part of T-COFFEE.*/
1263 /**/
1264 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1265 /*    it under the terms of the GNU General Public License as published by*/
1266 /*    the Free Software Foundation; either version 2 of the License, or*/
1267 /*    (at your option) any later version.*/
1268 /**/
1269 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1270 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1271 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1272 /*    GNU General Public License for more details.*/
1273 /**/
1274 /*    You should have received a copy of the GNU General Public License*/
1275 /*    along with Foobar; if not, write to the Free Software*/
1276 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1277 /*...............................................                                                                                      |*/
1278 /*  If you need some more information*/
1279 /*  cedric.notredame@europe.com*/
1280 /*...............................................                                                                                                                                     |*/
1281 /**/
1282 /**/
1283 /*      */
1284 /******************************COPYRIGHT NOTICE*******************************/