IA32 and x64 tcoffee binaries are added
[jabaws.git] / binaries / src / tcoffee / t_coffee_source / reformat_struc.c
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <stdarg.h>
5 #include <string.h>
6 #include <stddef.h>
7 #include <ctype.h>
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
15 #define FATAL "fatal:reformat_struc"
16
17 char * process_repeat (char *aln, char *seq, char *pdb)
18 {
19   char *tf, *file, *name;
20   Alignment *A, *A2;
21   Sequence *S, *P;
22   int r1, r2, is_r1, is_r2, l1=0, l2=0, a;
23   int *pos;
24   FILE *fp;
25   
26   A=main_read_aln (aln, NULL);
27   for (a=0; a<A->nseq; a++)ungap(A->seq_al[a]);
28   
29   S=main_read_seq (seq);
30   P=get_pdb_sequence (pdb);
31   
32
33   A2=align_two_sequences (S->seq[0], P->seq[0], "pam250mt", -10, -1, "myers_miller_pair_wise");
34   
35   pos=vcalloc ( A2->len_aln+1, sizeof (int));
36
37   for (l1=0, l2=0,a=0; a< A2->len_aln; a++)
38     {
39       
40       r1=A2->seq_al[0][a];
41       r2=A2->seq_al[1][a];
42       
43       is_r1=1-is_gap(r1);
44       is_r2=1-is_gap(r2);
45       
46       l1+=is_r1;
47       l2+=is_r2;
48       
49       if (!is_r1);
50       else
51         {
52           pos[l1]=l2;
53         }
54     }
55   tf=vtmpnam(NULL);
56   fp=vfopen (tf, "w");
57   for (a=0; a<A->nseq; a++)
58     {
59       int *coor, b, c;
60       
61       name=A->name[a];
62       file=vtmpnam (NULL);
63       coor=string2num_list2 (name, "-");
64
65       //Check the compatibility between the guide sequence and the coordinates
66       for ( c=0,b=coor[1]-1; b<coor[2]; b++, c++)
67         {
68
69           if (tolower(A->seq_al[a][c])!=tolower(S->seq[0][b]))
70             printf_exit (EXIT_FAILURE, stderr, "Incompatibility between the repeat [%s] and the master Sequence [%s]\n%s",A->name[a], seq, A->seq_al[a]);
71         }
72    
73       printf_system ( "extract_from_pdb %s -coor %d %d -seq_field SEQRES> %s", pdb,pos[coor[1]-1],pos[coor[2]-1], file);
74       fprintf (fp, ">%s _P_ %s\n", name, file);
75       vfree (coor);
76     }
77   vfclose (fp);
78   return tf;
79 }
80   
81   
82   
83 char *     normalize_pdb_file  (char *name_in, char *seq, char *out_file)
84 {
85   char command[1000];
86   Sequence *S;
87   Alignment *A;
88   int a;
89   int start, end, npdb, r1, r2;
90   char name[100];
91
92   
93  if ( !name_in) return NULL;
94  else
95    {
96      sprintf ( name, "%s", name_in);
97    }
98  
99  if ( !is_pdb_file(name))
100     {
101       fprintf(stdout, "\nERROR[normalize_pdb_file]: %s is not a pdb file[FATAL:%s]\n", name, PROGRAM);
102       myexit (EXIT_FAILURE);
103     }
104
105   S=get_pdb_sequence (name);
106   A=align_two_sequences (S->seq[0],seq,"idmat",-3,0, "fasta_pair_wise");
107  
108   
109
110   for (start=-1, end=-1,npdb=0,a=0; a< A->len_aln; a++)
111     {
112       r1=1-is_gap(A->seq_al[0][a]);
113       r2=1-is_gap(A->seq_al[1][a]);
114       
115       npdb+=r1;
116
117       if (r1 && r2 && start==-1)start=npdb;
118       if (r1 && r2)end=npdb;
119     }
120
121   free_aln(A);  
122   free_sequence (S, -1);
123   
124   sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -coor %d %d -nodiagnostic > %s", check_file_exists(name), start, end, out_file);
125   my_system ( command);
126   return out_file;
127   }
128
129 Ca_trace * trim_ca_trace (Ca_trace *T, char *seq )
130 {
131   /*This function:
132     -removes from Ca trace all the residues that are not in the sequence
133     -add in the Ca trace residues unmatched in the structure (it gives them a NULL structure)
134   */
135   Alignment *ALN;
136   Atom *A;  
137   int a,l, s, r, is_r, is_s;
138   int *seq_cache, *struc_cache;
139   int tot_l=0;
140
141   char buf1[10000];
142   char buf2[10000];
143
144   /*  lower_string (T->seq);
145       lower_string (seq);
146   */
147
148
149   sprintf (buf1, "%s", T->seq);
150   sprintf (buf2, "%s", seq);
151   lower_string (buf1);
152   lower_string (buf2);
153   
154
155   if ( strm (buf1,buf2))return T;
156   else
157     {
158      ALN=align_two_sequences (T->seq,seq, "est_idmat",-1, 0,"fasta_pair_wise"); 
159      struc_cache=vcalloc (ALN->len_aln+1, sizeof (int));
160      seq_cache  =vcalloc (ALN->len_aln+1, sizeof (int));
161      
162      for ( r=0, s=0,a=0; a< ALN->len_aln; a++)
163        {
164          is_r=!is_gap(ALN->seq_al[0][a]);
165          is_s=!is_gap(ALN->seq_al[1][a]);
166          
167          r+=is_r;
168          s+=is_s;
169          
170          if ( is_s && is_r)
171            {
172              struc_cache[r-1]=s-1;
173              seq_cache[s-1]=r-1;
174            }
175          else if ( is_s && !is_r)
176            {
177              seq_cache[s-1]=-1;
178              
179            }
180          else if ( !is_s && is_r)
181            {
182              struc_cache[r-1]=-1;
183            }
184        }
185      
186       T->ca=vrealloc ( T->ca, sizeof (Atom*)*(ALN->len_aln+1));
187       T->peptide_chain=vrealloc ( T->peptide_chain, (sizeof (Amino_acid*))*(ALN->len_aln+1));
188       T->seq=vrealloc ( T->seq,   ALN->len_aln+1);
189
190       for ( a=T->len; a< ALN->len_aln; a++)
191         {
192           T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));
193         }
194       
195       
196       /*Read every atom*/
197       for ( a=0; a< T->n_atom; a++)
198         {
199           
200           A=(T->structure[a]);
201           if ( struc_cache[A->res_num-1]==-1)continue;
202           else
203             {
204               /*set the struc residue to its sequence index*/
205               A->res_num=struc_cache[A->res_num-1]+1;
206               if (strm (A->type, "CA")) {T->ca[A->res_num-1]=A;tot_l++;}
207               if ( strm (A->type, "CA"))(T->peptide_chain[A->res_num-1])->CA=A;
208               if ( strm (A->type, "C"))(T->peptide_chain[A->res_num-1] )->C=A;
209               if ( strm (A->type, "CB"))(T->peptide_chain[A->res_num-1])->CB=A;
210               if ( strm (A->type, "N"))(T->peptide_chain[A->res_num-1] )->N=A;
211             }
212         }
213  
214       l=strlen(seq);     
215       for ( a=0;a< l; a++)
216        {
217
218          if ( seq_cache[a]==-1)
219            {
220              tot_l++;
221              T->ca[a]=NULL;
222              
223              if (!T->peptide_chain[a])T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));
224              T->peptide_chain[a]->CA=NULL;
225              T->peptide_chain[a]->C =NULL;
226              T->peptide_chain[a]->CB=NULL;
227              T->peptide_chain[a]->N=NULL;
228              T->seq[a]='x';
229            }
230          else
231            {
232              T->seq[a]=seq[a];
233            }
234        }
235       T->len=ALN->len_aln;
236      
237       /*
238         T->len=tot_l;
239       */
240       
241      free_aln (ALN);
242      vfree(seq_cache);
243      vfree(struc_cache);
244      
245           
246     }
247   return T;
248 }
249
250 Ca_trace * read_ca_trace (char *name, char *seq_field )
251 {
252   char *tp_name=NULL;
253   char command[10000];
254   
255
256   if ( !is_simple_pdb_file (name))
257     {
258       tp_name=vtmpnam (NULL);
259       sprintf ( command, "extract_from_pdb -seq_field %s  -infile %s -atom ALL -chain FIRST -mode simple> %s",seq_field, check_file_exists(name), tp_name);
260       if ( getenv4debug ("DEBUG_EXTRACT_FROM_PDB"))fprintf ( stderr, "\n[DEBUG_EXTRACT_FROM_PDB:read_ca_trace] %s\n", command);
261       my_system ( command);
262     }
263   else
264     tp_name=name;
265   
266   return simple_read_ca_trace (tp_name);
267 }
268
269 Ca_trace * simple_read_ca_trace (char *tp_name )
270     {
271         /*This function reads a pdb file into a Ca_trace structure*/
272          
273          
274
275         int a, c, n;
276         FILE *fp;
277         Atom *A;
278         char res;
279         char *buf;
280         Ca_trace *T=NULL;
281         int res_num=0, last_res_num=0;
282
283         
284         buf=vcalloc ( VERY_LONG_STRING, sizeof (char));
285         n=count_n_line_in_file (tp_name );
286         
287         if ( !T)
288             {
289             T=vcalloc ( 1, sizeof ( Ca_trace));
290             declare_name (T->name);
291             }
292         
293         /*1 Get the complete sequence: replace missing residues with Xs*/
294         for (a=0; a< VERY_LONG_STRING; a++)buf[a]='x';
295         res=res_num=0;
296         
297         fp=vfopen (tp_name, "r");
298         while ( (c=fgetc(fp))!='>');
299         fscanf ( fp, "%*s" );
300         while ( (c=fgetc(fp))!='\n');
301         fscanf ( fp, "%*s" );
302         while ( (c=fgetc(fp))!='\n');   
303         while ((c=fgetc(fp))!=EOF)
304           {
305                ungetc(c, fp);
306
307                fscanf (fp, "%*s %*s %c %*c %d %*f %*f %*f\n",&res,&res_num);
308                if (res)
309                  {
310                   
311                    res=tolower (res);
312                    buf[res_num-1]=res;
313                    last_res_num=res_num;
314                  }
315                res=res_num=0;
316           }
317         buf[last_res_num]='\0';
318         vfclose (fp);
319         /*Sequence Read*/
320         
321
322         T->len=strlen (buf);
323         T->seq=vcalloc ( T->len+1, sizeof (char));
324         buf=lower_string (buf);
325         sprintf ( T->seq, "%s", buf);
326         n+=T->len;
327         T->structure=vcalloc ( n, sizeof (Atom*));
328         for ( a=0; a< n; a++)T->structure[a]=vcalloc ( 1, sizeof (Atom));
329         T->ca=vcalloc ( T->len+1, sizeof ( Atom*));
330         a=0;
331
332         fp=vfopen (tp_name, "r");
333         while ( (c=fgetc(fp))!='>');
334         fscanf ( fp, "%*s" );
335         while ( (c=fgetc(fp))!='\n');
336         fscanf ( fp, "%*s" );
337         while ( (c=fgetc(fp))!='\n');   
338         
339         while ((c=fgetc(fp))!=EOF)
340            {
341                ungetc(c, fp);
342                A=T->structure[a];
343                A->num=a;
344                fscanf (fp, "%*s %s %s %*c %d %f %f %f\n",A->type, A->res,&A->res_num, &A->x, &A->y, &A->z);
345                res=A->res[0];
346                
347                res=tolower (res);
348                
349                T->seq[A->res_num-1]=res;
350
351                if (( strm ( A->type, "CA")) || ( strm ( A->type, "C3'")))
352                            T->ca[A->res_num-1]=A;
353                a++;
354            }
355         T->n_atom=a;
356
357         T->peptide_chain=vcalloc (T->len+1, sizeof (Amino_acid*));
358         for ( a=0; a<=T->len; a++) T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));                                   
359         for ( a=0; a< T->n_atom; a++)
360           {
361             A=T->structure[a];
362
363             if ( strm (A->type, "CA"))(T->peptide_chain[A->res_num-1])->CA=A;
364                 if ( strm (A->type, "C3'"))(T->peptide_chain[A->res_num-1])->CA=A;
365             if ( strm (A->type, "C"))(T->peptide_chain[A->res_num-1] )->C=A;
366             if ( strm (A->type, "CB"))(T->peptide_chain[A->res_num-1])->CB=A;
367             if ( strm (A->type, "N"))(T->peptide_chain[A->res_num-1] )->N=A;
368           }
369             
370
371         vfclose (fp);   
372         vfree (buf);
373         return T;
374     }
375 Ca_trace * hasch_ca_trace    ( Ca_trace *T)
376     {
377       
378       T=hasch_ca_trace_nb (T);
379       T=hasch_ca_trace_bubble (T);
380       T=hasch_ca_trace_transversal (T);
381       return T;
382     }
383 Ca_trace * hasch_ca_trace_transversal ( Ca_trace *TRACE)
384     {
385         /*This function gets the Coordinates of a protein and computes the distance of each Ca to its 
386           
387           given a Ca,
388           Compute the distance between, CA-x and CA+x with x=[1-N_ca]
389           T->nb[a][0]-->Number of distances.
390           T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
391           T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
392         */
393         
394         int a, b, d;
395         float dist;
396         Atom *A, *B;
397                 
398         Struct_nb *T;
399         Pdb_param *PP;
400
401
402         TRACE->Transversal=vcalloc ( 1, sizeof (Struct_nb));
403         
404         T=TRACE->Transversal;
405         PP=TRACE->pdb_param;
406
407         if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
408         if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
409
410         for (d=0,a=0; a< TRACE->len; a++)
411             {
412
413                 for ( b=1; b<=PP->N_ca; b++)
414                     {
415                     if ( (a-b)<0 || (a+b)>=TRACE->len)continue;
416                     A=TRACE->ca[a-b];
417                     B=TRACE->ca[a+b];
418                     dist=get_atomic_distance ( A, B);
419                     
420                     T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
421                     T->nb[a][T->nb[a][0]]=b;
422
423                     T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
424                     T->d_nb[a][T->nb[a][0]]=dist;
425                     
426                     d++;
427                     }
428                 T->max_nb=MAX (T->max_nb, T->nb[a][0]);
429             }
430         return TRACE;
431
432     }                          
433
434 Ca_trace * hasch_ca_trace_nb ( Ca_trace *TRACE)
435     {
436         /*This function gets the Coordinates of a protein and computes the distance of each Ca to its 
437           T->N_ca Ngb.
438           The first Ngb to the left and to the right are excluded
439           Ngd to the left get negative distances
440           Ngb to the right receive positive distances
441           T->nb[a][0]-->Number of ngb.
442           T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
443           T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
444         */
445         
446
447
448         int a, b, d;
449         float dist;
450         Atom *A, *B;
451                 
452         Struct_nb *T;
453         Pdb_param *PP;
454         
455         TRACE->Chain=vcalloc ( 1, sizeof (Struct_nb));
456         
457         T=TRACE->Chain;
458         PP=TRACE->pdb_param;
459
460         if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
461         if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
462
463         for (d=0,a=0; a< TRACE->len; a++)
464             {
465                 for ( b=MAX(0,a-PP->N_ca); b< MIN( a+PP->N_ca, TRACE->len); b++)
466                         {
467                         if (FABS(a-b)<2)continue;
468                         A=TRACE->ca[a];
469                         B=TRACE->ca[b];
470                         if ( !A || !B)continue;
471                         dist=get_atomic_distance ( A, B);
472                         if (b<a) dist=-dist;
473                 
474
475                         T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
476                         T->nb[a][T->nb[a][0]]=b;
477                         
478                         T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
479                         T->d_nb[a][T->nb[a][0]]=dist;
480                         d++;
481                         }
482                     
483                 T->max_nb=MAX (T->max_nb, T->nb[a][0]);
484             }
485         return TRACE;
486
487     }   
488 Ca_trace * hasch_ca_trace_bubble ( Ca_trace *TRACE)
489     {
490       
491         int a, b;
492         float dist;
493         Atom *A, *B;
494         float **list;
495         Pdb_param *PP;
496         Struct_nb *T;
497
498         PP=TRACE->pdb_param;
499         TRACE->Bubble=vcalloc ( 1, sizeof (Struct_nb));
500         T=TRACE->Bubble;
501         
502                 
503         
504         if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
505         if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
506         list=declare_float ( TRACE->n_atom, 3);
507                 
508
509         for (a=0; a< TRACE->len; a++)
510             {
511                 for ( b=0; b< TRACE->len; b++)
512                     {
513                         A=TRACE->ca[a];
514                         B=TRACE->ca[b];
515                         if ( !A || !B)continue;
516                         dist=get_atomic_distance ( A, B);
517                         
518                         if ( dist<PP->maximum_distance && FABS((A->res_num-B->res_num))>2)
519                            {
520                              T->nb[a][0]++;                                                    
521                              T->nb[a]=vrealloc ( T->nb[a], (T->nb[a][0]+1)*sizeof (int));
522                              T->nb[a][T->nb[a][0]]=(TRACE->ca[b])->num;
523                              
524                              T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
525                              T->d_nb[a][T->nb[a][0]]= ((a<b)?dist:-dist);
526                            }                                            
527                     }
528                 T->max_nb=MAX (T->max_nb, T->nb[a][0]);
529         
530             }
531         
532         for ( a=0; a< TRACE->len; a++)
533             {
534             for ( b=0; b< T->nb[a][0]; b++)
535                 {
536                     list[b][0]=T->nb[a][b+1];
537                     list[b][1]=T->d_nb[a][b+1];
538                     list[b][2]=(TRACE->structure[T->nb[a][b+1]])->res_num;
539                 }
540             
541             sort_float ( list, 3,2, 0, T->nb[a][0]-1);
542             for ( b=0; b< T->nb[a][0]; b++)
543                 {
544                     T->nb[a][b+1]=list[b][0];
545                     T->d_nb[a][b+1]=list[b][1];
546                 }
547             }
548
549         free_float ( list, -1);
550         return TRACE;
551
552     }
553
554        
555 float ** measure_ca_distances(Ca_trace *T)
556    {
557        int a, b;
558        Atom *A, *B;
559        float **dist;
560        
561        dist=declare_float ( T->len, T->len);
562     
563        for (a=0; a< T->len-1; a++)
564             {
565                 for ( b=a+1; b< T->len; b++)
566                     {
567                       
568                         A=T->ca[a];
569                         B=T->ca[b];
570                         dist[a][b]=dist[b][a]=get_atomic_distance ( A,      B);
571                         
572                     }
573             }
574        return dist;
575    }
576 float    get_atomic_distance ( Atom *A, Atom*B)
577    {
578        float dx, dy, dz, d;
579        
580        if ( !A || !B)
581          {
582            return UNDEFINED;
583          }
584
585        
586        dx=A->x - B->x;
587        dy=A->y - B->y;
588        dz=A->z - B->z;
589        
590
591        d=(float) sqrt ( (double) ( dx*dx +dy*dy +dz*dz));
592        return d;
593    }
594
595 char * map_contacts ( char  *file1, char *file2, float T)
596 {
597   
598   Ca_trace *ST1, *ST2;
599   int *contact_list;
600   int a;
601   
602
603   ST1=read_ca_trace (file1, "ATOM");
604   ST2=read_ca_trace (file2, "ATOM");
605
606   contact_list=identify_contacts (ST1, ST2, T);
607   for ( a=0; a<ST1->len; a++)
608     {
609       ST1->seq[a]=(contact_list[a]==1)?toupper(ST1->seq[a]):tolower(ST1->seq[a]);
610     }
611
612   return ST1->seq;
613 }
614
615 float ** print_contacts ( char  *file1, char *file2, float T)
616 {
617   
618   Ca_trace *ST1, *ST2;
619   float **dist, d;
620   int a, b;
621   Atom *A, *B;
622   char *list;
623   int *list1=NULL, *list2=NULL;
624   int *cache1, *cache2;
625   
626   
627    
628   if ((list=strstr (file1, "_R_"))!=NULL)
629       {
630         list[0]='\0';
631         list+=3;
632         list1=string2num_list2(list, "_");
633       }
634   
635   if ((list=strstr (file2, "_R_"))!=NULL)
636       {
637         list[0]='\0';
638         list+=3;
639         list2=string2num_list2(list, "_");
640       } 
641
642   fprintf ( stdout, "\n#%s (%s) struc2contacts_01", PROGRAM, VERSION);
643   fprintf ( stdout, "\nStructure %s vs %s", file1, file2);
644   ST1=read_ca_trace (file1, "SEQRES");
645   ST2=read_ca_trace (file2, "SEQRES");
646   
647   
648   cache1=vcalloc (ST1->len+1, sizeof (int));
649   cache2=vcalloc (ST2->len+1, sizeof (int));
650   
651   if (list1)for ( a=1; a<list1[0]; a++)cache1[list1[a]]=1;
652   else for ( a=1; a<ST1->len; a++)cache1[a]=1;
653   
654   if (list2)for ( a=1; a<list2[0]; a++)cache2[list2[a]]=1;
655   else for ( a=1; a<ST2->len; a++)cache2[a]=1;
656   
657   
658   dist=declare_float (ST1->len+1,ST2->len+1);
659   vfree (list1); vfree(list2);
660   
661   for ( a=0; a< ST1->n_atom; a++)
662     {
663       A=ST1->structure[a];
664       if ( !cache1[A->res_num])continue;
665       for ( b=0; b<ST2->n_atom; b++)
666             {
667               
668               B=ST2->structure[b];
669               if( !cache2[B->res_num])continue;
670             
671               d=get_atomic_distance (A,B);
672               
673               if (dist[A->res_num][B->res_num]==0 || dist[A->res_num][B->res_num]>d)dist[A->res_num][B->res_num]=d;
674             }
675     }
676
677   for ( a=1; a<=ST1->len; a++)
678     {
679       A=ST1->ca[a-1];
680       if ( !A || !cache1[A->res_num])continue;
681       for ( b=1; b<= ST2->len; b++)
682         {
683           B=ST2->ca[b-1];
684           if( !B || !cache2[B->res_num])continue;
685           if(dist[a][b]!=0)fprintf ( stdout, "\nResidue %3d [%s] vs %3d [%s] %9.4f Angstrom",A->res_num,A->res,B->res_num,B->res,dist[a][b]);
686         }
687     }
688   fprintf ( stdout, "\n");
689   vfree (cache1);vfree(cache2);
690   free_float (dist, -1);
691   return NULL;
692 }
693   
694
695 int * identify_contacts (Ca_trace *ST1,Ca_trace *ST2, float T) 
696 {
697   int a, b;
698   float d;
699   int *result;
700   
701   
702     
703   result=vcalloc ( ST1->len+1, sizeof (int));
704
705   
706   for ( a=0; a< ST1->n_atom; a++)
707     for ( b=0; b<ST2->n_atom; b++)
708   
709       {
710         
711         d=get_atomic_distance (ST1->structure[a], ST2->structure[b]);
712         if (d<T)
713           {
714             /*Atom *A, *B;
715             A=ST1->structure[a]; B=ST2->structure[b];
716             fprintf ( stderr, "\n%d %s %s Vs %d %s %s: %f", A->res_num, A->res, A->type, B->res_num, B->res, B->type, d); */
717             result[(ST1->structure[a])->res_num-1]=1;
718           }
719       }
720   return result;
721 }
722
723 Sequence *seq2contacts ( Sequence *S, float T)
724 {
725   int a;
726   Sequence *NS;
727   
728  
729   
730   NS=duplicate_sequence (S);
731   for ( a=0; a< S->nseq; a++)
732     {
733       NS->seq[a]=string2contacts ( S->seq[a], S->name[a], S->seq_comment[a], T);
734     }
735   
736   return NS;
737 }
738
739 char *string2contacts (char *seq,char *name, char *comment, float T)
740 {
741   char **nlist;
742   char *r;
743   char *result;
744   int a, b, n;
745   char *struc_name;
746   static char *struc_file;
747   static char *ligand_file;
748   Alignment *A;
749   char r0, r1;
750   
751   int l, ln;
752   char command[1000];
753   /*>seq__struc Ligand1 Chain1 Ligand 2 cahin2
754     Chain: index or ANY if unknown
755     Ligand: name of pdb file
756   */
757   
758   if ( !struc_file)
759     {
760       struc_file=vtmpnam  (NULL);
761       ligand_file=vtmpnam (NULL);
762     }
763   
764
765   
766   result=vcalloc ( strlen (seq)+1, sizeof (char));
767   for ( a=0; a< strlen (seq); a++)result[a]='0';
768   
769   nlist=string2list (comment);
770   if ( !nlist)return result;
771   else
772     n=atoi(nlist[0]);
773
774   struc_name=strstr(name, "_S_");
775   if (!struc_name && !is_pdb_struc (name))
776     {
777       
778       return result;
779     }
780   else if ( !struc_name && is_pdb_struc (name))
781     {
782       struc_name=name;
783     }
784   else
785     {
786       struc_name+=3;
787       if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
788         {
789          sprintf (command, "cp %s %s", name, struc_file);
790         }
791       else
792         {
793           sprintf ( command, "extract_from_pdb -infile %s -atom ALL -mode simple -force >%s",name, struc_file);
794         }
795       my_system (command);
796     }
797   
798   
799   
800   for ( a=1, ln=1;a<n; ln++)
801         {
802           fprintf ( stderr, "\n\tProcess: Structure %s Against Ligand %s T=%.2f Angstrom",struc_name, nlist[a], T);
803          
804           
805           if (check_file_exists ( nlist[a]) && is_simple_pdb_file(nlist[a]))
806               {         
807                 sprintf (command, "cp %s %s", nlist[a], ligand_file);
808                 a++;
809               }
810           else if ( check_file_exists ( nlist[a]))
811             {
812               sprintf (command, "extract_from_pdb -infile %s -ligand ALL -atom ALL -mode simple -force >%s", nlist[a], ligand_file);
813               a++;
814             }
815           else
816             {
817               sprintf ( command, "extract_from_pdb -infile %s -chain %s -ligand %s -ligand_only -atom ALL -mode simple -force >%s",struc_name, nlist[a+1],nlist[a], ligand_file);
818               a+=2;           
819             }
820           my_system (command);
821           
822           if ( T>0)
823             {
824               r=map_contacts (struc_file,ligand_file,T);
825               
826               toggle_case_in_align_two_sequences (KEEP_CASE);
827               A=align_two_sequences (seq,r,"pam250mt", -10, -1, "myers_miller_pair_wise");     
828               toggle_case_in_align_two_sequences (CHANGE_CASE);
829               
830               
831               for ( l=0,b=0; b< A->len_aln; b++)
832                 {
833                   r0=A->seq_al[0][b];r1=A->seq_al[1][b];
834                   if (!is_gap(r0))
835                     {
836                       if (isupper(r1))result[l]=(result[l]!='0')?'9':'0'+ln;
837                       l++;
838                     }
839                 }
840               
841               free_aln (A);
842               fprintf ( stderr, " [DONE]");
843             }
844           else if ( T==0)
845             {
846               print_contacts( struc_file,ligand_file,T);
847             }
848           
849         }
850       fprintf ( stderr, "\n");
851       
852       return result;
853 }
854
855
856 char **struclist2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
857 {
858   char *list, **pdb;
859   int a;
860   char **R, *tmpf;
861   
862   tmpf=vtmpnam (NULL);
863   
864   list=strstr ( comment, "_P_")+3;
865   if ( !strstr (comment, "_P_"))return NULL;
866   else
867     {
868       pdb=string2list ( strstr ( comment, "_P_")+3);
869     }
870   
871   for (a=1; a<atoi(pdb[0]); a++)
872     {
873       char com[1000];
874       sprintf ( com, "_P_ %s", pdb[a]);
875       R=struc2nb (name, seq, com, Threshold, atom, tmpf);
876     }
877   if (output==NULL|| strstr (output, "fasta"))
878     {
879       for (a=0; a< strlen (seq); a++)
880         if (isupper (R[a][a]))fprintf (stdout, ">%s R=%d T=%.2f %s\n%s\n", name, a+1, Threshold, comment, R[a]);
881     }
882   else
883     {
884       FILE *fp;
885       char c;
886       
887       fp=vfopen (tmpf, "r");
888       while ( (c=fgetc(fp))!=EOF)fprintf (stdout, "%c", c);
889       vfclose (fp);
890     }
891   return NULL;
892 }
893
894 char **struc2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
895 {
896   char *struc_file;
897   char *struc_name;
898   Ca_trace *T;
899   Atom *A1, *A2;
900   int a, b;
901   short **hasch;
902   FILE *fp;
903   float d;
904   char command[10000];
905   static char **R;
906   
907   struc_file=vtmpnam (NULL);
908   declare_name (struc_name);
909   
910   sscanf ( strstr(comment, "_P_"), "_P_ %s", struc_name);
911   //struc_name=strstr(name, "_S_");
912   
913   if (!R)
914     {
915       int l;
916       l=strlen (seq);
917       R=declare_char (l+1, l+1);
918       for ( a=0; a<l; a++)
919         {
920           sprintf ( R[a], "%s", seq);
921           lower_string (R[a]);
922         }
923     }
924   if ( atom) atom=translate_string (atom, "_", " ");
925   if ( !struc_name) return NULL;
926   
927   else 
928     {
929       //struc_name+=3;
930       if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
931         {
932           sprintf (command, "cp %s %s", name, struc_file);
933         }
934       else
935         {
936           sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",struc_name,(atom==NULL)?"ALL":atom, struc_file);
937         }
938   
939       my_system (command);
940     }
941   T=read_ca_trace (struc_file, "ATOM");
942   hasch=declare_short (T->len, T->len);
943   
944   if (!R)
945     {
946       int l;
947       l=strlen (seq);
948       R=declare_char (l+1, l+1);
949       for ( a=0; a<l; a++)
950         {
951           sprintf ( R[a], "%s", seq);
952           lower_string (R[a]);
953         }
954     }
955   for ( a=0; a< T->n_atom; a++)
956     for ( b=0; b< T->n_atom; b++)
957       {
958         A1=T->structure[a];A2=T->structure[b];
959         d=get_atomic_distance (A1, A2);
960
961         if ( d<Threshold)
962           {
963             hasch[A1->res_num-1][A2->res_num-1]=1;
964           }
965       }
966   fp=vfopen (output, "a");
967   fprintf ( fp, "#Distance_map_format_01\n#Sequence %s with T= %.2f", struc_name, Threshold);
968   for ( a=0; a<T->len; a++)
969     {
970       int c;
971       c=change_residue_coordinate ( T->seq,seq,a);
972       
973       if ( c!=-1 && seq)
974         {
975           char r1, r2;
976           r1=(T->seq)[a];
977           r2=seq[c];
978           r1=tolower(r1);r2=tolower(r2);
979           if ( r1!=r2) continue;
980           R[c][c]=toupper (R[c][c]);
981           fprintf (fp, "\n%s Residue %d ",struc_name,c+1);
982           for ( b=0; b<T->len; b++)
983             {
984               int d;
985               char r3, r4;
986               r3=(T->seq)[a];r4=seq[c];
987               r3=tolower(r3);r4=tolower(r4);
988               if ( r3!=r4) continue;
989               
990               d=change_residue_coordinate (T->seq,seq,b);
991               
992               if ( hasch[a][b] && d!=-1)
993                 {
994                   fprintf (fp, "%d ",d+1);
995                   R[c][d]=toupper (R[c][d]);
996                 }
997             }
998           fprintf (fp, ";");
999         }
1000     }
1001   free_short (hasch, -1);
1002   fprintf (fp, "\n");
1003   return R;
1004 }
1005
1006 short **seq2nb (char *seq, char *pdb, float Threshold, char *atom)
1007 {
1008   char *struc_file;
1009  
1010   Ca_trace *T;
1011   Atom *A1, *A2;
1012   int a, b;
1013   short **hasch;
1014   short **result; //Contains the result for every residue of seq
1015   float d;
1016   char command[10000];
1017
1018   // get a clean pdb file
1019   struc_file=vtmpnam (NULL);
1020   if ( check_file_exists (struc_file) && is_simple_pdb_file(struc_file))
1021     {
1022       sprintf (command, "cp %s %s", pdb, struc_file);
1023     }
1024   else
1025     {
1026       sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",pdb,(atom==NULL)?"ALL":atom, struc_file);
1027     }
1028   my_system (command);
1029   
1030   //Read and hasch the PDB file
1031   T=read_ca_trace (struc_file, "ATOM");
1032   hasch=declare_short (T->len, T->len);
1033   result=declare_short (strlen (seq)+1, strlen (seq)+1);
1034   for ( a=0; a< T->n_atom; a++)
1035     for ( b=0; b< T->n_atom; b++)
1036       {
1037         A1=T->structure[a];A2=T->structure[b];
1038         d=get_atomic_distance (A1, A2);
1039
1040         if ( d<Threshold)
1041           {
1042             hasch[A1->res_num-1][A2->res_num-1]=1;
1043           }
1044       }
1045   
1046   for ( a=0; a<T->len; a++)
1047     {
1048       int c;
1049       c=change_residue_coordinate ( T->seq,seq,a);
1050       if ( c!=-1)
1051         {
1052           for ( b=0; b<T->len; b++)
1053             {
1054               int d;
1055               d=change_residue_coordinate (T->seq,seq,b);
1056               
1057               if ( hasch[a][b] && d!=-1)
1058                 {
1059                   result[c][result[c][0]++]=d;
1060                 }
1061             }
1062         }
1063     }
1064   free_short (hasch, -1);
1065   return result;
1066 }
1067 /******************************COPYRIGHT NOTICE*******************************/
1068 /*© Centro de Regulacio Genomica */
1069 /*and */
1070 /*Cedric Notredame */
1071 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
1072 /*All rights reserved.*/
1073 /*This file is part of T-COFFEE.*/
1074 /**/
1075 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1076 /*    it under the terms of the GNU General Public License as published by*/
1077 /*    the Free Software Foundation; either version 2 of the License, or*/
1078 /*    (at your option) any later version.*/
1079 /**/
1080 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1081 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1082 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1083 /*    GNU General Public License for more details.*/
1084 /**/
1085 /*    You should have received a copy of the GNU General Public License*/
1086 /*    along with Foobar; if not, write to the Free Software*/
1087 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1088 /*...............................................                                                                                      |*/
1089 /*  If you need some more information*/
1090 /*  cedric.notredame@europe.com*/
1091 /*...............................................                                                                                                                                     |*/
1092 /**/
1093 /**/
1094 /*      */
1095 /******************************COPYRIGHT NOTICE*******************************/