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*******************************/