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")) T->ca[A->res_num-1]=A;
 352                a++;
 353            }
 354         T->n_atom=a;
 355
 356         T->peptide_chain=vcalloc (T->len+1, sizeof (Amino_acid*));
 357         for ( a=0; a<=T->len; a++) T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));
 358         for ( a=0; a< T->n_atom; a++)
 359           {
 360             A=T->structure[a];
 361
 362             if ( strm (A->type, "CA"))(T->peptide_chain[A->res_num-1])->CA=A;
 363             if ( strm (A->type, "C"))(T->peptide_chain[A->res_num-1] )->C=A;
 364             if ( strm (A->type, "CB"))(T->peptide_chain[A->res_num-1])->CB=A;
 365             if ( strm (A->type, "N"))(T->peptide_chain[A->res_num-1] )->N=A;
 366           }
 367
 368
 369         vfclose (fp);
 370         vfree (buf);
 371         return T;
 372     }
 373 Ca_trace * hasch_ca_trace    ( Ca_trace *T)
 374     {
 375
 376       T=hasch_ca_trace_nb (T);
 377       T=hasch_ca_trace_bubble (T);
 378       T=hasch_ca_trace_transversal (T);
 379       return T;
 380     }
 381 Ca_trace * hasch_ca_trace_transversal ( Ca_trace *TRACE)
 382     {
 383         /*This function gets the Coordinates of a protein and computes the distance of each Ca to its
 384
 385           given a Ca,
 386           Compute the distance between, CA-x and CA+x with x=[1-N_ca]
 387           T->nb[a][0]-->Number of distances.
 388           T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
 389           T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
 390         */
 391
 392         int a, b, d;
 393         float dist;
 394         Atom *A, *B;
 395
 396         Struct_nb *T;
 397         Pdb_param *PP;
 398
 399
 400         TRACE->Transversal=vcalloc ( 1, sizeof (Struct_nb));
 401
 402         T=TRACE->Transversal;
 403         PP=TRACE->pdb_param;
 404
 405         if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
 406         if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
 407
 408         for (d=0,a=0; a< TRACE->len; a++)
 409             {
 410
 411                 for ( b=1; b<=PP->N_ca; b++)
 412                     {
 413                     if ( (a-b)<0 || (a+b)>=TRACE->len)continue;
 414                     A=TRACE->ca[a-b];
 415                     B=TRACE->ca[a+b];
 416                     dist=get_atomic_distance ( A, B);
 417
 418                     T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
 419                     T->nb[a][T->nb[a][0]]=b;
 420
 421                     T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
 422                     T->d_nb[a][T->nb[a][0]]=dist;
 423
 424                     d++;
 425                     }
 426                 T->max_nb=MAX (T->max_nb, T->nb[a][0]);
 427             }
 428         return TRACE;
 429
 430     }
 431
 432 Ca_trace * hasch_ca_trace_nb ( Ca_trace *TRACE)
 433     {
 434         /*This function gets the Coordinates of a protein and computes the distance of each Ca to its
 435           T->N_ca Ngb.
 436           The first Ngb to the left and to the right are excluded
 437           Ngd to the left get negative distances
 438           Ngb to the right receive positive distances
 439           T->nb[a][0]-->Number of ngb.
 440           T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
 441           T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
 442         */
 443
 444
 445
 446         int a, b, d;
 447         float dist;
 448         Atom *A, *B;
 449
 450         Struct_nb *T;
 451         Pdb_param *PP;
 452
 453         TRACE->Chain=vcalloc ( 1, sizeof (Struct_nb));
 454
 455         T=TRACE->Chain;
 456         PP=TRACE->pdb_param;
 457
 458         if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
 459         if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
 460
 461         for (d=0,a=0; a< TRACE->len; a++)
 462             {
 463                 for ( b=MAX(0,a-PP->N_ca); b< MIN( a+PP->N_ca, TRACE->len); b++)
 464                         {
 465                         if (FABS(a-b)<2)continue;
 466                         A=TRACE->ca[a];
 467                         B=TRACE->ca[b];
 468                         if ( !A || !B)continue;
 469                         dist=get_atomic_distance ( A, B);
 470                         if (b<a) dist=-dist;
 471
 472
 473                         T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
 474                         T->nb[a][T->nb[a][0]]=b;
 475
 476                         T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
 477                         T->d_nb[a][T->nb[a][0]]=dist;
 478                         d++;
 479                         }
 480
 481                 T->max_nb=MAX (T->max_nb, T->nb[a][0]);
 482             }
 483         return TRACE;
 484
 485     }
 486 Ca_trace * hasch_ca_trace_bubble ( Ca_trace *TRACE)
 487     {
 488
 489         int a, b;
 490         float dist;
 491         Atom *A, *B;
 492         float **list;
 493         Pdb_param *PP;
 494         Struct_nb *T;
 495
 496         PP=TRACE->pdb_param;
 497         TRACE->Bubble=vcalloc ( 1, sizeof (Struct_nb));
 498         T=TRACE->Bubble;
 499
 500
 501
 502         if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
 503         if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
 504         list=declare_float ( TRACE->n_atom, 3);
 505
 506
 507         for (a=0; a< TRACE->len; a++)
 508             {
 509                 for ( b=0; b< TRACE->len; b++)
 510                     {
 511                         A=TRACE->ca[a];
 512                         B=TRACE->ca[b];
 513                         if ( !A || !B)continue;
 514                         dist=get_atomic_distance ( A, B);
 515
 516                         if ( dist<PP->maximum_distance && FABS((A->res_num-B->res_num))>2)
 517                            {
 518                              T->nb[a][0]++;
 519                              T->nb[a]=vrealloc ( T->nb[a], (T->nb[a][0]+1)*sizeof (int));
 520                              T->nb[a][T->nb[a][0]]=(TRACE->ca[b])->num;
 521
 522                              T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
 523                              T->d_nb[a][T->nb[a][0]]= ((a<b)?dist:-dist);
 524                            }
 525                     }
 526                 T->max_nb=MAX (T->max_nb, T->nb[a][0]);
 527
 528             }
 529
 530         for ( a=0; a< TRACE->len; a++)
 531             {
 532             for ( b=0; b< T->nb[a][0]; b++)
 533                 {
 534                     list[b][0]=T->nb[a][b+1];
 535                     list[b][1]=T->d_nb[a][b+1];
 536                     list[b][2]=(TRACE->structure[T->nb[a][b+1]])->res_num;
 537                 }
 538
 539             sort_float ( list, 3,2, 0, T->nb[a][0]-1);
 540             for ( b=0; b< T->nb[a][0]; b++)
 541                 {
 542                     T->nb[a][b+1]=list[b][0];
 543                     T->d_nb[a][b+1]=list[b][1];
 544                 }
 545             }
 546
 547         free_float ( list, -1);
 548         return TRACE;
 549
 550     }
 551
 552
 553 float ** measure_ca_distances(Ca_trace *T)
 554    {
 555        int a, b;
 556        Atom *A, *B;
 557        float **dist;
 558
 559        dist=declare_float ( T->len, T->len);
 560
 561        for (a=0; a< T->len-1; a++)
 562             {
 563                 for ( b=a+1; b< T->len; b++)
 564                     {
 565
 566                         A=T->ca[a];
 567                         B=T->ca[b];
 568                         dist[a][b]=dist[b][a]=get_atomic_distance ( A,      B);
 569
 570                     }
 571             }
 572        return dist;
 573    }
 574 float    get_atomic_distance ( Atom *A, Atom*B)
 575    {
 576        float dx, dy, dz, d;
 577
 578        if ( !A || !B)
 579          {
 580            return UNDEFINED;
 581          }
 582
 583
 584        dx=A->x - B->x;
 585        dy=A->y - B->y;
 586        dz=A->z - B->z;
 587
 588
 589        d=(float) sqrt ( (double) ( dx*dx +dy*dy +dz*dz));
 590        return d;
 591    }
 592
 593 char * map_contacts ( char  *file1, char *file2, float T)
 594 {
 595
 596   Ca_trace *ST1, *ST2;
 597   int *contact_list;
 598   int a;
 599
 600
 601   ST1=read_ca_trace (file1, "ATOM");
 602   ST2=read_ca_trace (file2, "ATOM");
 603
 604   contact_list=identify_contacts (ST1, ST2, T);
 605   for ( a=0; a<ST1->len; a++)
 606     {
 607       ST1->seq[a]=(contact_list[a]==1)?toupper(ST1->seq[a]):tolower(ST1->seq[a]);
 608     }
 609
 610   return ST1->seq;
 611 }
 612
 613 float ** print_contacts ( char  *file1, char *file2, float T)
 614 {
 615
 616   Ca_trace *ST1, *ST2;
 617   float **dist, d;
 618   int a, b;
 619   Atom *A, *B;
 620   char *list;
 621   int *list1=NULL, *list2=NULL;
 622   int *cache1, *cache2;
 623
 624
 625
 626   if ((list=strstr (file1, "_R_"))!=NULL)
 627       {
 628         list[0]='\0';
 629         list+=3;
 630         list1=string2num_list2(list, "_");
 631       }
 632
 633   if ((list=strstr (file2, "_R_"))!=NULL)
 634       {
 635         list[0]='\0';
 636         list+=3;
 637         list2=string2num_list2(list, "_");
 638       }
 639
 640   fprintf ( stdout, "\n#%s (%s) struc2contacts_01", PROGRAM, VERSION);
 641   fprintf ( stdout, "\nStructure %s vs %s", file1, file2);
 642   ST1=read_ca_trace (file1, "SEQRES");
 643   ST2=read_ca_trace (file2, "SEQRES");
 644
 645
 646   cache1=vcalloc (ST1->len+1, sizeof (int));
 647   cache2=vcalloc (ST2->len+1, sizeof (int));
 648
 649   if (list1)for ( a=1; a<list1[0]; a++)cache1[list1[a]]=1;
 650   else for ( a=1; a<ST1->len; a++)cache1[a]=1;
 651
 652   if (list2)for ( a=1; a<list2[0]; a++)cache2[list2[a]]=1;
 653   else for ( a=1; a<ST2->len; a++)cache2[a]=1;
 654
 655
 656   dist=declare_float (ST1->len+1,ST2->len+1);
 657   vfree (list1); vfree(list2);
 658
 659   for ( a=0; a< ST1->n_atom; a++)
 660     {
 661       A=ST1->structure[a];
 662       if ( !cache1[A->res_num])continue;
 663       for ( b=0; b<ST2->n_atom; b++)
 664             {
 665
 666               B=ST2->structure[b];
 667               if( !cache2[B->res_num])continue;
 668
 669               d=get_atomic_distance (A,B);
 670
 671               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;
 672             }
 673     }
 674
 675   for ( a=1; a<=ST1->len; a++)
 676     {
 677       A=ST1->ca[a-1];
 678       if ( !A || !cache1[A->res_num])continue;
 679       for ( b=1; b<= ST2->len; b++)
 680         {
 681           B=ST2->ca[b-1];
 682           if( !B || !cache2[B->res_num])continue;
 683           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]);
 684         }
 685     }
 686   fprintf ( stdout, "\n");
 687   vfree (cache1);vfree(cache2);
 688   free_float (dist, -1);
 689   return NULL;
 690 }
 691
 692
 693 int * identify_contacts (Ca_trace *ST1,Ca_trace *ST2, float T)
 694 {
 695   int a, b;
 696   float d;
 697   int *result;
 698
 699
 700
 701   result=vcalloc ( ST1->len+1, sizeof (int));
 702
 703
 704   for ( a=0; a< ST1->n_atom; a++)
 705     for ( b=0; b<ST2->n_atom; b++)
 706
 707       {
 708
 709         d=get_atomic_distance (ST1->structure[a], ST2->structure[b]);
 710         if (d<T)
 711           {
 712             /*Atom *A, *B;
 713             A=ST1->structure[a]; B=ST2->structure[b];
 714             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); */
 715             result[(ST1->structure[a])->res_num-1]=1;
 716           }
 717       }
 718   return result;
 719 }
 720
 721 Sequence *seq2contacts ( Sequence *S, float T)
 722 {
 723   int a;
 724   Sequence *NS;
 725
 726
 727
 728   NS=duplicate_sequence (S);
 729   for ( a=0; a< S->nseq; a++)
 730     {
 731       NS->seq[a]=string2contacts ( S->seq[a], S->name[a], S->seq_comment[a], T);
 732     }
 733
 734   return NS;
 735 }
 736
 737 char *string2contacts (char *seq,char *name, char *comment, float T)
 738 {
 739   char **nlist;
 740   char *r;
 741   char *result;
 742   int a, b, n;
 743   char *struc_name;
 744   static char *struc_file;
 745   static char *ligand_file;
 746   Alignment *A;
 747   char r0, r1;
 748
 749   int l, ln;
 750   char command[1000];
 751   /*>seq__struc Ligand1 Chain1 Ligand 2 cahin2
 752     Chain: index or ANY if unknown
 753     Ligand: name of pdb file
 754   */
 755
 756   if ( !struc_file)
 757     {
 758       struc_file=vtmpnam  (NULL);
 759       ligand_file=vtmpnam (NULL);
 760     }
 761
 762
 763
 764   result=vcalloc ( strlen (seq)+1, sizeof (char));
 765   for ( a=0; a< strlen (seq); a++)result[a]='0';
 766
 767   nlist=string2list (comment);
 768   if ( !nlist)return result;
 769   else
 770     n=atoi(nlist[0]);
 771
 772   struc_name=strstr(name, "_S_");
 773   if (!struc_name && !is_pdb_struc (name))
 774     {
 775
 776       return result;
 777     }
 778   else if ( !struc_name && is_pdb_struc (name))
 779     {
 780       struc_name=name;
 781     }
 782   else
 783     {
 784       struc_name+=3;
 785       if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
 786         {
 787          sprintf (command, "cp %s %s", name, struc_file);
 788         }
 789       else
 790         {
 791           sprintf ( command, "extract_from_pdb -infile %s -atom ALL -mode simple -force >%s",name, struc_file);
 792         }
 793       my_system (command);
 794     }
 795
 796
 797
 798   for ( a=1, ln=1;a<n; ln++)
 799         {
 800           fprintf ( stderr, "\n\tProcess: Structure %s Against Ligand %s T=%.2f Angstrom",struc_name, nlist[a], T);
 801
 802
 803           if (check_file_exists ( nlist[a]) && is_simple_pdb_file(nlist[a]))
 804               {
 805                 sprintf (command, "cp %s %s", nlist[a], ligand_file);
 806                 a++;
 807               }
 808           else if ( check_file_exists ( nlist[a]))
 809             {
 810               sprintf (command, "extract_from_pdb -infile %s -ligand ALL -atom ALL -mode simple -force >%s", nlist[a], ligand_file);
 811               a++;
 812             }
 813           else
 814             {
 815               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);
 816               a+=2;
 817             }
 818           my_system (command);
 819
 820           if ( T>0)
 821             {
 822               r=map_contacts (struc_file,ligand_file,T);
 823
 824               toggle_case_in_align_two_sequences (KEEP_CASE);
 825               A=align_two_sequences (seq,r,"pam250mt", -10, -1, "myers_miller_pair_wise");
 826               toggle_case_in_align_two_sequences (CHANGE_CASE);
 827
 828
 829               for ( l=0,b=0; b< A->len_aln; b++)
 830                 {
 831                   r0=A->seq_al[0][b];r1=A->seq_al[1][b];
 832                   if (!is_gap(r0))
 833                     {
 834                       if (isupper(r1))result[l]=(result[l]!='0')?'9':'0'+ln;
 835                       l++;
 836                     }
 837                 }
 838
 839               free_aln (A);
 840               fprintf ( stderr, " [DONE]");
 841             }
 842           else if ( T==0)
 843             {
 844               print_contacts( struc_file,ligand_file,T);
 845             }
 846
 847         }
 848       fprintf ( stderr, "\n");
 849
 850       return result;
 851 }
 852
 853
 854 char **struclist2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
 855 {
 856   char *list, **pdb;
 857   int a;
 858   char **R, *tmpf;
 859
 860   tmpf=vtmpnam (NULL);
 861
 862   list=strstr ( comment, "_P_")+3;
 863   if ( !strstr (comment, "_P_"))return NULL;
 864   else
 865     {
 866       pdb=string2list ( strstr ( comment, "_P_")+3);
 867     }
 868
 869   for (a=1; a<atoi(pdb[0]); a++)
 870     {
 871       char com[1000];
 872       sprintf ( com, "_P_ %s", pdb[a]);
 873       R=struc2nb (name, seq, com, Threshold, atom, tmpf);
 874     }
 875   if (output==NULL|| strstr (output, "fasta"))
 876     {
 877       for (a=0; a< strlen (seq); a++)
 878         if (isupper (R[a][a]))fprintf (stdout, ">%s R=%d T=%.2f %s\n%s\n", name, a+1, Threshold, comment, R[a]);
 879     }
 880   else
 881     {
 882       FILE *fp;
 883       char c;
 884
 885       fp=vfopen (tmpf, "r");
 886       while ( (c=fgetc(fp))!=EOF)fprintf (stdout, "%c", c);
 887       vfclose (fp);
 888     }
 889   return NULL;
 890 }
 891
 892 char **struc2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
 893 {
 894   char *struc_file;
 895   char *struc_name;
 896   Ca_trace *T;
 897   Atom *A1, *A2;
 898   int a, b;
 899   short **hasch;
 900   FILE *fp;
 901   float d;
 902   char command[10000];
 903   static char **R;
 904
 905   struc_file=vtmpnam (NULL);
 906   declare_name (struc_name);
 907
 908   sscanf ( strstr(comment, "_P_"), "_P_ %s", struc_name);
 909   //struc_name=strstr(name, "_S_");
 910
 911   if (!R)
 912     {
 913       int l;
 914       l=strlen (seq);
 915       R=declare_char (l+1, l+1);
 916       for ( a=0; a<l; a++)
 917         {
 918           sprintf ( R[a], "%s", seq);
 919           lower_string (R[a]);
 920         }
 921     }
 922   if ( atom) atom=translate_string (atom, "_", " ");
 923   if ( !struc_name) return NULL;
 924
 925   else
 926     {
 927       //struc_name+=3;
 928       if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
 929         {
 930           sprintf (command, "cp %s %s", name, struc_file);
 931         }
 932       else
 933         {
 934           sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",struc_name,(atom==NULL)?"ALL":atom, struc_file);
 935         }
 936
 937       my_system (command);
 938     }
 939   T=read_ca_trace (struc_file, "ATOM");
 940   hasch=declare_short (T->len, T->len);
 941
 942   if (!R)
 943     {
 944       int l;
 945       l=strlen (seq);
 946       R=declare_char (l+1, l+1);
 947       for ( a=0; a<l; a++)
 948         {
 949           sprintf ( R[a], "%s", seq);
 950           lower_string (R[a]);
 951         }
 952     }
 953   for ( a=0; a< T->n_atom; a++)
 954     for ( b=0; b< T->n_atom; b++)
 955       {
 956         A1=T->structure[a];A2=T->structure[b];
 957         d=get_atomic_distance (A1, A2);
 958
 959         if ( d<Threshold)
 960           {
 961             hasch[A1->res_num-1][A2->res_num-1]=1;
 962           }
 963       }
 964   fp=vfopen (output, "a");
 965   fprintf ( fp, "#Distance_map_format_01\n#Sequence %s with T= %.2f", struc_name, Threshold);
 966   for ( a=0; a<T->len; a++)
 967     {
 968       int c;
 969       c=change_residue_coordinate ( T->seq,seq,a);
 970
 971       if ( c!=-1 && seq)
 972         {
 973           char r1, r2;
 974           r1=(T->seq)[a];
 975           r2=seq[c];
 976           r1=tolower(r1);r2=tolower(r2);
 977           if ( r1!=r2) continue;
 978           R[c][c]=toupper (R[c][c]);
 979           fprintf (fp, "\n%s Residue %d ",struc_name,c+1);
 980           for ( b=0; b<T->len; b++)
 981             {
 982               int d;
 983               char r3, r4;
 984               r3=(T->seq)[a];r4=seq[c];
 985               r3=tolower(r3);r4=tolower(r4);
 986               if ( r3!=r4) continue;
 987
 988               d=change_residue_coordinate (T->seq,seq,b);
 989
 990               if ( hasch[a][b] && d!=-1)
 991                 {
 992                   fprintf (fp, "%d ",d+1);
 993                   R[c][d]=toupper (R[c][d]);
 994                 }
 995             }
 996           fprintf (fp, ";");
 997         }
 998     }
 999   free_short (hasch, -1);
1000   fprintf (fp, "\n");
1001   return R;
1002 }
1003
1004 short **seq2nb (char *seq, char *pdb, float Threshold, char *atom)
1005 {
1006   char *struc_file;
1007   char *struc_name;
1008   Ca_trace *T;
1009   Atom *A1, *A2;
1010   int a, b;
1011   short **hasch;
1012   short **result; //Contains the result for every residue of seq
1013   float d;
1014   char command[10000];
1015
1016   // get a clean pdb file
1017   struc_file=vtmpnam (NULL);
1018   if ( check_file_exists (struc_file) && is_simple_pdb_file(struc_name))
1019     {
1020       sprintf (command, "cp %s %s", pdb, struc_file);
1021     }
1022   else
1023     {
1024       sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",pdb,(atom==NULL)?"ALL":atom, struc_file);
1025     }
1026   my_system (command);
1027
1028   //Read and hasch the PDB file
1029   T=read_ca_trace (struc_file, "ATOM");
1030   hasch=declare_short (T->len, T->len);
1031   result=declare_short (strlen (seq)+1, strlen (seq)+1);
1032   for ( a=0; a< T->n_atom; a++)
1033     for ( b=0; b< T->n_atom; b++)
1034       {
1035         A1=T->structure[a];A2=T->structure[b];
1036         d=get_atomic_distance (A1, A2);
1037
1038         if ( d<Threshold)
1039           {
1040             hasch[A1->res_num-1][A2->res_num-1]=1;
1041           }
1042       }
1043
1044   for ( a=0; a<T->len; a++)
1045     {
1046       int c;
1047       c=change_residue_coordinate ( T->seq,seq,a);
1048       if ( c!=-1)
1049         {
1050           for ( b=0; b<T->len; b++)
1051             {
1052               int d;
1053               d=change_residue_coordinate (T->seq,seq,b);
1054
1055               if ( hasch[a][b] && d!=-1)
1056                 {
1057                   result[c][result[c][0]++]=d;
1058                 }
1059             }
1060         }
1061     }
1062   free_short (hasch, -1);
1063   return result;
1064 }
1065 /*********************************COPYRIGHT NOTICE**********************************/
1066 /*© Centro de Regulacio Genomica */
1067 /*and */
1068 /*Cedric Notredame */
1069 /*Tue Oct 27 10:12:26 WEST 2009. */
1070 /*All rights reserved.*/
1071 /*This file is part of T-COFFEE.*/
1072 /**/
1073 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1074 /*    it under the terms of the GNU General Public License as published by*/
1075 /*    the Free Software Foundation; either version 2 of the License, or*/
1076 /*    (at your option) any later version.*/
1077 /**/
1078 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1079 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1080 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1081 /*    GNU General Public License for more details.*/
1082 /**/
1083 /*    You should have received a copy of the GNU General Public License*/
1084 /*    along with Foobar; if not, write to the Free Software*/
1085 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1086 /*...............................................                                                                                      |*/
1087 /*  If you need some more information*/
1088 /*  cedric.notredame@europe.com*/
1089 /*...............................................                                                                                                                                     |*/
1090 /**/
1091 /**/
1092 /*      */
1093 /*********************************COPYRIGHT NOTICE**********************************/