Logo Search packages:      
Sourcecode: t-coffee version File versions  Download package

reformat_struc.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <string.h>
#include <stddef.h>

#include "io_lib_header.h"
#include "util_lib_header.h"
#include "dp_lib_header.h" 
#include "define_header.h"



#define FATAL "fatal:reformat_struc"


char *     normalize_pdb_file  (char *name_in, char *seq, char *out_file)
{
  char command[1000];
  Sequence *S;
  Alignment *A;
  int a;
  int start, end, npdb, r1, r2;
  char name[100];

  
 if ( !name_in) return NULL;
 else
   {
     sprintf ( name, "%s", name_in);
   }
 
 if ( !is_pdb_file(name))
    {
      fprintf(stderr, "\nERROR[normalize_pdb_file]: %s is not a pdb file[FATAL:%s]\n", name, PROGRAM);
    }

  S=get_pdb_sequence (name);
  A=align_two_sequences (S->seq[0],seq,"idmat",-3,0, "fasta_pair_wise");
 
  

  for (start=-1, end=-1,npdb=0,a=0; a< A->len_aln; a++)
    {
      r1=1-is_gap(A->seq_al[0][a]);
      r2=1-is_gap(A->seq_al[1][a]);
      
      npdb+=r1;

      if (r1 && r2 && start==-1)start=npdb;
      if (r1 && r2)end=npdb;
    }

  free_aln(A);  
  free_sequence (S, -1);
  
  sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -coor %d %d -nodiagnostic > %s", check_file_exists(name), start, end, out_file);
  my_system ( command);
  return out_file;
  }

Ca_trace * trim_ca_trace (Ca_trace *T, char *seq )
{
  /*This function:
    -removes from Ca trace all the residues that are not in the sequence
    -add in the Ca trace residues unmatched in the structure (it gives them a NULL structure)
  */
  Alignment *ALN;
  Atom *A;  
  int a,l, s, r, is_r, is_s;
  int *seq_cache, *struc_cache;
  int tot_l=0;

  char buf1[10000];
  char buf2[10000];

  /*  lower_string (T->seq);
      lower_string (seq);
  */


  sprintf (buf1, "%s", T->seq);
  sprintf (buf2, "%s", seq);
  lower_string (buf1);
  lower_string (buf2);
  

  if ( strm (buf1,buf2))return T;
  else
    {
     ALN=align_two_sequences (T->seq,seq, "est_idmat",-1, 0,"fasta_pair_wise"); 
     struc_cache=vcalloc (ALN->len_aln+1, sizeof (int));
     seq_cache  =vcalloc (ALN->len_aln+1, sizeof (int));
     
     for ( r=0, s=0,a=0; a< ALN->len_aln; a++)
       {
       is_r=!is_gap(ALN->seq_al[0][a]);
       is_s=!is_gap(ALN->seq_al[1][a]);
       
       r+=is_r;
       s+=is_s;
       
       if ( is_s && is_r)
         {
           struc_cache[r-1]=s-1;
           seq_cache[s-1]=r-1;
         }
       else if ( is_s && !is_r)
         {
           seq_cache[s-1]=-1;
           
         }
       else if ( !is_s && is_r)
         {
           struc_cache[r-1]=-1;
         }
       }
     
      T->ca=vrealloc ( T->ca, sizeof (Atom*)*(ALN->len_aln+1));
      T->peptide_chain=vrealloc ( T->peptide_chain, (sizeof (Amino_acid*))*(ALN->len_aln+1));
      T->seq=vrealloc ( T->seq,   ALN->len_aln+1);

      for ( a=T->len; a< ALN->len_aln; a++)
      {
        T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));
      }
      
      
      /*Read every atom*/
      for ( a=0; a< T->n_atom; a++)
      {
        
        A=(T->structure[a]);
        if ( struc_cache[A->res_num-1]==-1)continue;
        else
          {
            /*set the struc residue to its sequence index*/
            A->res_num=struc_cache[A->res_num-1]+1;
            if (strm (A->type, "CA")) {T->ca[A->res_num-1]=A;tot_l++;}
            if ( strm (A->type, "CA"))(T->peptide_chain[A->res_num-1])->CA=A;
            if ( strm (A->type, "C"))(T->peptide_chain[A->res_num-1] )->C=A;
            if ( strm (A->type, "CB"))(T->peptide_chain[A->res_num-1])->CB=A;
            if ( strm (A->type, "N"))(T->peptide_chain[A->res_num-1] )->N=A;
          }
      }
 
      l=strlen(seq);     
      for ( a=0;a< l; a++)
       {

       if ( seq_cache[a]==-1)
         {
           tot_l++;
           T->ca[a]=NULL;
           
           if (!T->peptide_chain[a])T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));
           T->peptide_chain[a]->CA=NULL;
           T->peptide_chain[a]->C =NULL;
           T->peptide_chain[a]->CB=NULL;
           T->peptide_chain[a]->N=NULL;
           T->seq[a]='x';
         }
       else
         {
           T->seq[a]=seq[a];
         }
       }
      T->len=ALN->len_aln;
     
      /*
      T->len=tot_l;
      */
      
     free_aln (ALN);
     vfree(seq_cache);
     vfree(struc_cache);
     
        
    }
  return T;
}

Ca_trace * read_ca_trace (char *name, char *seq_field )
{
  char *tp_name=NULL;
  char command[10000];
  
  if ( !is_simple_pdb_file (name))
    {
      tp_name=vtmpnam (NULL);
      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);
      my_system ( command);
    }
  else
    tp_name=name;
  
  return simple_read_ca_trace (tp_name);
}

Ca_trace * simple_read_ca_trace (char *tp_name )
    {
      /*This function reads a pdb file into a Ca_trace structure*/
       
       

      int a, c, n;
      FILE *fp;
      Atom *A;
      char res;
      char *buf;
      Ca_trace *T=NULL;
      int res_num, last_res_num;

      
      buf=vcalloc ( VERY_LONG_STRING, sizeof (char));
      n=count_n_line_in_file (tp_name );
      
      if ( !T)
          {
          T=vcalloc ( 1, sizeof ( Ca_trace));
          declare_name (T->name);
          }
      
      /*1 Get the complete sequence: replace missing residues with Xs*/
      for (a=0; a< VERY_LONG_STRING; a++)buf[a]='x';
      res=res_num=0;
      
      fp=vfopen (tp_name, "r");
        while ( (c=fgetc(fp))!='>');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');     
      while ((c=fgetc(fp))!=EOF)
        {
             ungetc(c, fp);

             fscanf (fp, "%*s %*s %c %*c %d %*f %*f %*f\n",&res,&res_num);
             if (res)
             {
              
               res=tolower (res);
               buf[res_num-1]=res;
               last_res_num=res_num;
             }
             res=res_num=0;
        }
      buf[last_res_num]='\0';
      vfclose (fp);
      /*Sequence Read*/
      

      T->len=strlen (buf);
      T->seq=vcalloc ( T->len+1, sizeof (char));
      buf=lower_string (buf);
      sprintf ( T->seq, "%s", buf);
      n+=T->len;
      T->structure=vcalloc ( n, sizeof (Atom*));
      for ( a=0; a< n; a++)T->structure[a]=vcalloc ( 1, sizeof (Atom));
      T->ca=vcalloc ( T->len+1, sizeof ( Atom*));
      a=0;

      fp=vfopen (tp_name, "r");
        while ( (c=fgetc(fp))!='>');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');     
      
      while ((c=fgetc(fp))!=EOF)
         {
             ungetc(c, fp);
             A=T->structure[a];
             A->num=a;
             fscanf (fp, "%*s %s %s %*c %d %f %f %f\n",A->type, A->res,&A->res_num, &A->x, &A->y, &A->z);
             res=A->res[0];
             
             res=tolower (res);
             
             T->seq[A->res_num-1]=res;

             if ( strm ( A->type, "CA")) T->ca[A->res_num-1]=A;
             a++;
         }
      T->n_atom=a;

      T->peptide_chain=vcalloc (T->len+1, sizeof (Amino_acid*));
      for ( a=0; a<=T->len; a++) T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));                         
      for ( a=0; a< T->n_atom; a++)
        {
          A=T->structure[a];

          if ( strm (A->type, "CA"))(T->peptide_chain[A->res_num-1])->CA=A;
          if ( strm (A->type, "C"))(T->peptide_chain[A->res_num-1] )->C=A;
          if ( strm (A->type, "CB"))(T->peptide_chain[A->res_num-1])->CB=A;
          if ( strm (A->type, "N"))(T->peptide_chain[A->res_num-1] )->N=A;
        }
          

      vfclose (fp);     
      vfree (buf);
      return T;
    }
Ca_trace * hasch_ca_trace    ( Ca_trace *T)
    {
      
      T=hasch_ca_trace_nb (T);
      T=hasch_ca_trace_bubble (T);
      T=hasch_ca_trace_transversal (T);
      return T;
    }
Ca_trace * hasch_ca_trace_transversal ( Ca_trace *TRACE)
    {
      /*This function gets the Coordinates of a protein and computes the distance of each Ca to its 
        
        given a Ca,
        Compute the distance between, CA-x and CA+x with x=[1-N_ca]
        T->nb[a][0]-->Number of distances.
        T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
        T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
      */
      
      int a, b, d;
      float dist;
      Atom *A, *B;
            
      Struct_nb *T;
      Pdb_param *PP;


      TRACE->Transversal=vcalloc ( 1, sizeof (Struct_nb));
      
      T=TRACE->Transversal;
      PP=TRACE->pdb_param;

      if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
      if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);

      for (d=0,a=0; a< TRACE->len; a++)
          {

              for ( b=1; b<=PP->N_ca; b++)
                {
                if ( (a-b)<0 || (a+b)>=TRACE->len)continue;
                A=TRACE->ca[a-b];
                B=TRACE->ca[a+b];
                dist=get_atomic_distance ( A, B);
                
                T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
                T->nb[a][T->nb[a][0]]=b;

                T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
                T->d_nb[a][T->nb[a][0]]=dist;
                
                d++;
                }
            T->max_nb=MAX (T->max_nb, T->nb[a][0]);
          }
      return TRACE;

    }                    

Ca_trace * hasch_ca_trace_nb ( Ca_trace *TRACE)
    {
      /*This function gets the Coordinates of a protein and computes the distance of each Ca to its 
        T->N_ca Ngb.
        The first Ngb to the left and to the right are excluded
        Ngd to the left get negative distances
        Ngb to the right receive positive distances
        T->nb[a][0]-->Number of ngb.
        T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
        T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
      */
      


      int a, b, d;
      float dist;
      Atom *A, *B;
            
      Struct_nb *T;
      Pdb_param *PP;
      
      TRACE->Chain=vcalloc ( 1, sizeof (Struct_nb));
      
      T=TRACE->Chain;
      PP=TRACE->pdb_param;

      if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
      if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);

      for (d=0,a=0; a< TRACE->len; a++)
          {
            for ( b=MAX(0,a-PP->N_ca); b< MIN( a+PP->N_ca, TRACE->len); b++)
                    {
                  if (FABS(a-b)<2)continue;
                  A=TRACE->ca[a];
                  B=TRACE->ca[b];
                  if ( !A || !B)continue;
                  dist=get_atomic_distance ( A, B);
                  if (b<a) dist=-dist;
            

                  T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
                  T->nb[a][T->nb[a][0]]=b;
                  
                  T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
                  T->d_nb[a][T->nb[a][0]]=dist;
                  d++;
                  }
                
            T->max_nb=MAX (T->max_nb, T->nb[a][0]);
          }
      return TRACE;

    } 
Ca_trace * hasch_ca_trace_bubble ( Ca_trace *TRACE)
    {
      
      int a, b;
      float dist;
      Atom *A, *B;
      float **list;
      Pdb_param *PP;
      Struct_nb *T;

      PP=TRACE->pdb_param;
      TRACE->Bubble=vcalloc ( 1, sizeof (Struct_nb));
      T=TRACE->Bubble;
      
            
      
      if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
      if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
      list=declare_float ( TRACE->n_atom, 3);
            

      for (a=0; a< TRACE->len; a++)
          {
            for ( b=0; b< TRACE->len; b++)
                {
                  A=TRACE->ca[a];
                  B=TRACE->ca[b];
                  if ( !A || !B)continue;
                  dist=get_atomic_distance ( A, B);
                  
                  if ( dist<PP->maximum_distance && FABS((A->res_num-B->res_num))>2)
                     {
                       T->nb[a][0]++;                                                
                       T->nb[a]=vrealloc ( T->nb[a], (T->nb[a][0]+1)*sizeof (int));
                       T->nb[a][T->nb[a][0]]=(TRACE->ca[b])->num;
                       
                       T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
                       T->d_nb[a][T->nb[a][0]]= ((a<b)?dist:-dist);
                     }                                
                }
            T->max_nb=MAX (T->max_nb, T->nb[a][0]);
      
          }
      
      for ( a=0; a< TRACE->len; a++)
          {
          for ( b=0; b< T->nb[a][0]; b++)
              {
                list[b][0]=T->nb[a][b+1];
                list[b][1]=T->d_nb[a][b+1];
                list[b][2]=(TRACE->structure[T->nb[a][b+1]])->res_num;
            }
          
          sort_float ( list, 3,2, 0, T->nb[a][0]-1);
          for ( b=0; b< T->nb[a][0]; b++)
              {
                T->nb[a][b+1]=list[b][0];
                T->d_nb[a][b+1]=list[b][1];
            }
          }

      free_float ( list, -1);
      return TRACE;

    }

       
float ** measure_ca_distances(Ca_trace *T)
   {
       int a, b;
       Atom *A, *B;
       float **dist;
       
       dist=declare_float ( T->len, T->len);
    
       for (a=0; a< T->len-1; a++)
          {
            for ( b=a+1; b< T->len; b++)
                {
                  
                  A=T->ca[a];
                  B=T->ca[b];
                  dist[a][b]=dist[b][a]=get_atomic_distance ( A,      B);
                  
                }
          }
       return dist;
   }
float    get_atomic_distance ( Atom *A, Atom*B)
   {
       float dx, dy, dz, d;
       
       if ( !A || !B)
       {
         return UNDEFINED;
       }

       
       dx=A->x - B->x;
       dy=A->y - B->y;
       dz=A->z - B->z;
       

       d=(float) sqrt ( (double) ( dx*dx +dy*dy +dz*dz));
       return d;
   }

char * map_contacts ( char  *file1, char *file2, float T)
{
  
  Ca_trace *ST1, *ST2;
  int *contact_list;
  int a;
  

  ST1=read_ca_trace (file1, "ATOM");
  ST2=read_ca_trace (file2, "ATOM");

  contact_list=identify_contacts (ST1, ST2, T);
  for ( a=0; a<ST1->len; a++)
    {
      ST1->seq[a]=(contact_list[a]==1)?toupper(ST1->seq[a]):tolower(ST1->seq[a]);
    }

  return ST1->seq;
}

float ** print_contacts ( char  *file1, char *file2, float T)
{
  
  Ca_trace *ST1, *ST2;
  float **dist, d;
  int a, b;
  Atom *A, *B;
  char *list;
  int *list1=NULL, *list2=NULL;
  int *cache1, *cache2;
  
  
   
  if ((list=strstr (file1, "_R_"))!=NULL)
      {
      list[0]='\0';
      list+=3;
      list1=string2num_list2(list, "_");
      }
  
  if ((list=strstr (file2, "_R_"))!=NULL)
      {
      list[0]='\0';
      list+=3;
      list2=string2num_list2(list, "_");
      }     

  fprintf ( stdout, "\n#%s (%s) struc2contacts_01", PROGRAM, VERSION);
  fprintf ( stdout, "\nStructure %s vs %s", file1, file2);
  ST1=read_ca_trace (file1, "SEQRES");
  ST2=read_ca_trace (file2, "SEQRES");
  
  
  cache1=vcalloc (ST1->len+1, sizeof (int));
  cache2=vcalloc (ST2->len+1, sizeof (int));
  
  if (list1)for ( a=1; a<list1[0]; a++)cache1[list1[a]]=1;
  else for ( a=1; a<ST1->len; a++)cache1[a]=1;
  
  if (list2)for ( a=1; a<list2[0]; a++)cache2[list2[a]]=1;
  else for ( a=1; a<ST2->len; a++)cache2[a]=1;
  
  
  dist=declare_float (ST1->len+1,ST2->len+1);
  vfree (list1); vfree(list2);
  
  for ( a=0; a< ST1->n_atom; a++)
    {
      A=ST1->structure[a];
      if ( !cache1[A->res_num])continue;
      for ( b=0; b<ST2->n_atom; b++)
          {
            
            B=ST2->structure[b];
            if( !cache2[B->res_num])continue;
          
            d=get_atomic_distance (A,B);
            
            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;
          }
    }

  for ( a=1; a<=ST1->len; a++)
    {
      A=ST1->ca[a-1];
      if ( !A || !cache1[A->res_num])continue;
      for ( b=1; b<= ST2->len; b++)
      {
        B=ST2->ca[b-1];
        if( !B || !cache2[B->res_num])continue;
        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]);
      }
    }
  fprintf ( stdout, "\n");
  vfree (cache1);vfree(cache2);
  free_float (dist, -1);
  return NULL;
}
  

int * identify_contacts (Ca_trace *ST1,Ca_trace *ST2, float T) 
{
  int a, b;
  float d;
  int *result;
  
  
    
  result=vcalloc ( ST1->len+1, sizeof (int));

  
  for ( a=0; a< ST1->n_atom; a++)
    for ( b=0; b<ST2->n_atom; b++)
  
      {
      
      d=get_atomic_distance (ST1->structure[a], ST2->structure[b]);
      if (d<T)
        {
          /*Atom *A, *B;
          A=ST1->structure[a]; B=ST2->structure[b];
          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); */
          result[(ST1->structure[a])->res_num-1]=1;
        }
      }
  return result;
}

Sequence *seq2contacts ( Sequence *S, float T)
{
  int a;
  Sequence *NS;
  
 
  
  NS=duplicate_sequence (S);
  for ( a=0; a< S->nseq; a++)
    {
      NS->seq[a]=string2contacts ( S->seq[a], S->name[a], S->seq_comment[a], T);
    }
  
  return NS;
}

char *string2contacts (char *seq,char *name, char *comment, float T)
{
  char **nlist;
  char *r;
  char *result;
  int a, b, n;
  char *struc_name;
  static char *struc_file;
  static char *ligand_file;
  Alignment *A;
  char r0, r1;
  
  int l, ln;
  char command[1000];
  /*>seq__struc Ligand1 Chain1 Ligand 2 cahin2
    Chain: index or ANY if unknown
    Ligand: name of pdb file
  */
  
  if ( !struc_file)
    {
      struc_file=vtmpnam  (NULL);
      ligand_file=vtmpnam (NULL);
    }
  

  
  result=vcalloc ( strlen (seq)+1, sizeof (char));
  for ( a=0; a< strlen (seq); a++)result[a]='0';
  
  nlist=string2list (comment);
  if ( !nlist)return result;
  else
    n=atoi(nlist[0]);

  struc_name=strstr(name, "_S_");
  if (!struc_name && !is_pdb_struc (name))
    {
      
      return result;
    }
  else if ( !struc_name && is_pdb_struc (name))
    {
      struc_name=name;
    }
  else
    {
      struc_name+=3;
      if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
      {
       sprintf (command, "cp %s %s", name, struc_file);
      }
      else
      {
        sprintf ( command, "extract_from_pdb -infile %s -atom ALL -mode simple -force >%s",name, struc_file);
      }
      my_system (command);
    }
  
  
  
  for ( a=1, ln=1;a<n; ln++)
      {
        fprintf ( stderr, "\n\tProcess: Structure %s Against Ligand %s T=%.2f Angstrom",struc_name, nlist[a], T);
       
        
        if (check_file_exists ( nlist[a]) && is_simple_pdb_file(nlist[a]))
            {           
            sprintf (command, "cp %s %s", nlist[a], ligand_file);
            a++;
            }
        else if ( check_file_exists ( nlist[a]))
          {
            sprintf (command, "extract_from_pdb -infile %s -ligand ALL -atom ALL -mode simple -force >%s", nlist[a], ligand_file);
            a++;
          }
        else
          {
            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);
            a+=2;       
          }
        my_system (command);
        
        if ( T>0)
          {
            r=map_contacts (struc_file,ligand_file,T);
            
            toggle_case_in_align_two_sequences (KEEP_CASE);
            A=align_two_sequences (seq,r,"pam250mt", -10, -1, "myers_miller_pair_wise");     
            toggle_case_in_align_two_sequences (CHANGE_CASE);
            
            
            for ( l=0,b=0; b< A->len_aln; b++)
            {
              r0=A->seq_al[0][b];r1=A->seq_al[1][b];
              if (!is_gap(r0))
                {
                  if (isupper(r1))result[l]=(result[l]!='0')?'9':'0'+ln;
                  l++;
                }
            }
            
            free_aln (A);
            fprintf ( stderr, " [DONE]");
          }
        else if ( T==0)
          {
            print_contacts( struc_file,ligand_file,T);
          }
        
      }
      fprintf ( stderr, "\n");
      
      return result;
}

short **struc2nb (char *name,char *seq, char *comment, float Threshold, char *atom)
{
  char *struc_file;
  char *struc_name;
  Ca_trace *T;
  Atom *A1, *A2;
  int a, b;
  short **hasch;
  float d;
  char command[10000];
  
  struc_file=vtmpnam (NULL);
  struc_name=strstr(name, "_S_");
  if ( atom) atom=translate_string (atom, "_", " ");
  if ( !struc_name) return NULL;
  
  else 
    {
      struc_name+=3;
      if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
      {
        sprintf (command, "cp %s %s", name, struc_file);
      }
      else
      {
        sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",name,(atom==NULL)?"ALL":atom, struc_file);
      }
      my_system (command);
    }
  T=read_ca_trace (struc_file, "ATOM");
  hasch=declare_short (T->len, T->len);

  for ( a=0; a< T->n_atom; a++)
    for ( b=0; b< T->n_atom; b++)
      {
      A1=T->structure[a];A2=T->structure[b];
      d=get_atomic_distance (A1, A2);

      if ( d<Threshold)
        {
          hasch[A1->res_num-1][A2->res_num-1]=1;
        }
      }
  
  fprintf ( stdout, "#Distance_map_format_01\n#Sequence %s with T= %.2f", struc_name, Threshold);
  for ( a=0; a<T->len; a++)
    {
      int c;
      c=change_residue_coordinate ( T->seq,seq,a);
      if ( c!=-1)
      {
        fprintf ( stdout, "\n%s Residue %d ",struc_name,c+1);
        for ( b=0; b<T->len; b++)
          {
            int d;
            d=change_residue_coordinate (T->seq,seq,b);
            if ( hasch[a][b] && d!=-1)
            fprintf (stdout, "%d ",d+1);
          }
        fprintf (stdout, ";");
      }
    }
  fprintf (stdout, "\n");
  return hasch;
}
/*********************************COPYRIGHT NOTICE**********************************/
/* Centre National de la Recherche Scientifique (CNRS) */
/*and */
/*Cedric Notredame */
/*Tue May 10 12:08:44     2005. */
/*All rights reserved.*/
/*This file is part of T-COFFEE.*/
/**/
/*    T-COFFEE is free software; you can redistribute it and/or modify*/
/*    it under the terms of the GNU General Public License as published by*/
/*    the Free Software Foundation; either version 2 of the License, or*/
/*    (at your option) any later version.*/
/**/
/*    T-COFFEE is distributed in the hope that it will be useful,*/
/*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
/*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
/*    GNU General Public License for more details.*/
/**/
/*    You should have received a copy of the GNU General Public License*/
/*    along with Foobar; if not, write to the Free Software*/
/*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
/*...............................................                                                                                      |*/
/*  If you need some more information*/
/*  cedric.notredame@europe.com*/
/*...............................................                                                                                                                                     |*/
/**/
/**/
/*    */
/*********************************COPYRIGHT NOTICE**********************************/

Generated by  Doxygen 1.6.0   Back to index