+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <stdarg.h>
-#include <string.h>
-#include <stddef.h>
-#include <ctype.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 * process_repeat (char *aln, char *seq, char *pdb)
-{
- char *tf, *file, *name;
- Alignment *A, *A2;
- Sequence *S, *P;
- int r1, r2, is_r1, is_r2, l1=0, l2=0, a;
- int *pos;
- FILE *fp;
-
- A=main_read_aln (aln, NULL);
- for (a=0; a<A->nseq; a++)ungap(A->seq_al[a]);
-
- S=main_read_seq (seq);
- P=get_pdb_sequence (pdb);
-
-
- A2=align_two_sequences (S->seq[0], P->seq[0], "pam250mt", -10, -1, "myers_miller_pair_wise");
-
- pos=vcalloc ( A2->len_aln+1, sizeof (int));
-
- for (l1=0, l2=0,a=0; a< A2->len_aln; a++)
- {
-
- r1=A2->seq_al[0][a];
- r2=A2->seq_al[1][a];
-
- is_r1=1-is_gap(r1);
- is_r2=1-is_gap(r2);
-
- l1+=is_r1;
- l2+=is_r2;
-
- if (!is_r1);
- else
- {
- pos[l1]=l2;
- }
- }
- tf=vtmpnam(NULL);
- fp=vfopen (tf, "w");
- for (a=0; a<A->nseq; a++)
- {
- int *coor, b, c;
-
- name=A->name[a];
- file=vtmpnam (NULL);
- coor=string2num_list2 (name, "-");
-
- //Check the compatibility between the guide sequence and the coordinates
- for ( c=0,b=coor[1]-1; b<coor[2]; b++, c++)
- {
-
- if (tolower(A->seq_al[a][c])!=tolower(S->seq[0][b]))
- 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]);
- }
-
- printf_system ( "extract_from_pdb %s -coor %d %d -seq_field SEQRES> %s", pdb,pos[coor[1]-1],pos[coor[2]-1], file);
- fprintf (fp, ">%s _P_ %s\n", name, file);
- vfree (coor);
- }
- vfclose (fp);
- return tf;
-}
-
-
-
-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(stdout, "\nERROR[normalize_pdb_file]: %s is not a pdb file[FATAL:%s]\n", name, PROGRAM);
- myexit (EXIT_FAILURE);
- }
-
- 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);
- if ( getenv4debug ("DEBUG_EXTRACT_FROM_PDB"))fprintf ( stderr, "\n[DEBUG_EXTRACT_FROM_PDB:read_ca_trace] %s\n", command);
- 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=0, last_res_num=0;
-
-
- 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;
-}
-
-
-char **struclist2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
-{
- char *list, **pdb;
- int a;
- char **R, *tmpf;
-
- tmpf=vtmpnam (NULL);
-
- list=strstr ( comment, "_P_")+3;
- if ( !strstr (comment, "_P_"))return NULL;
- else
- {
- pdb=string2list ( strstr ( comment, "_P_")+3);
- }
-
- for (a=1; a<atoi(pdb[0]); a++)
- {
- char com[1000];
- sprintf ( com, "_P_ %s", pdb[a]);
- R=struc2nb (name, seq, com, Threshold, atom, tmpf);
- }
- if (output==NULL|| strstr (output, "fasta"))
- {
- for (a=0; a< strlen (seq); a++)
- if (isupper (R[a][a]))fprintf (stdout, ">%s R=%d T=%.2f %s\n%s\n", name, a+1, Threshold, comment, R[a]);
- }
- else
- {
- FILE *fp;
- char c;
-
- fp=vfopen (tmpf, "r");
- while ( (c=fgetc(fp))!=EOF)fprintf (stdout, "%c", c);
- vfclose (fp);
- }
- return NULL;
-}
-
-char **struc2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
-{
- char *struc_file;
- char *struc_name;
- Ca_trace *T;
- Atom *A1, *A2;
- int a, b;
- short **hasch;
- FILE *fp;
- float d;
- char command[10000];
- static char **R;
-
- struc_file=vtmpnam (NULL);
- declare_name (struc_name);
-
- sscanf ( strstr(comment, "_P_"), "_P_ %s", struc_name);
- //struc_name=strstr(name, "_S_");
-
- if (!R)
- {
- int l;
- l=strlen (seq);
- R=declare_char (l+1, l+1);
- for ( a=0; a<l; a++)
- {
- sprintf ( R[a], "%s", seq);
- lower_string (R[a]);
- }
- }
- 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",struc_name,(atom==NULL)?"ALL":atom, struc_file);
- }
-
- my_system (command);
- }
- T=read_ca_trace (struc_file, "ATOM");
- hasch=declare_short (T->len, T->len);
-
- if (!R)
- {
- int l;
- l=strlen (seq);
- R=declare_char (l+1, l+1);
- for ( a=0; a<l; a++)
- {
- sprintf ( R[a], "%s", seq);
- lower_string (R[a]);
- }
- }
- 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;
- }
- }
- fp=vfopen (output, "a");
- fprintf ( fp, "#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 && seq)
- {
- char r1, r2;
- r1=(T->seq)[a];
- r2=seq[c];
- r1=tolower(r1);r2=tolower(r2);
- if ( r1!=r2) continue;
- R[c][c]=toupper (R[c][c]);
- fprintf (fp, "\n%s Residue %d ",struc_name,c+1);
- for ( b=0; b<T->len; b++)
- {
- int d;
- char r3, r4;
- r3=(T->seq)[a];r4=seq[c];
- r3=tolower(r3);r4=tolower(r4);
- if ( r3!=r4) continue;
-
- d=change_residue_coordinate (T->seq,seq,b);
-
- if ( hasch[a][b] && d!=-1)
- {
- fprintf (fp, "%d ",d+1);
- R[c][d]=toupper (R[c][d]);
- }
- }
- fprintf (fp, ";");
- }
- }
- free_short (hasch, -1);
- fprintf (fp, "\n");
- return R;
-}
-
-short **seq2nb (char *seq, char *pdb, float Threshold, char *atom)
-{
- char *struc_file;
- char *struc_name;
- Ca_trace *T;
- Atom *A1, *A2;
- int a, b;
- short **hasch;
- short **result; //Contains the result for every residue of seq
- float d;
- char command[10000];
-
- // get a clean pdb file
- struc_file=vtmpnam (NULL);
- if ( check_file_exists (struc_file) && is_simple_pdb_file(struc_name))
- {
- sprintf (command, "cp %s %s", pdb, struc_file);
- }
- else
- {
- sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",pdb,(atom==NULL)?"ALL":atom, struc_file);
- }
- my_system (command);
-
- //Read and hasch the PDB file
- T=read_ca_trace (struc_file, "ATOM");
- hasch=declare_short (T->len, T->len);
- result=declare_short (strlen (seq)+1, strlen (seq)+1);
- 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;
- }
- }
-
- for ( a=0; a<T->len; a++)
- {
- int c;
- c=change_residue_coordinate ( T->seq,seq,a);
- if ( 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)
- {
- result[c][result[c][0]++]=d;
- }
- }
- }
- }
- free_short (hasch, -1);
- return result;
-}
-/*********************************COPYRIGHT NOTICE**********************************/
-/*© Centro de Regulacio Genomica */
-/*and */
-/*Cedric Notredame */
-/*Tue Oct 27 10:12:26 WEST 2009. */
-/*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**********************************/
git://source.jalview.org / jabaws.git / blobdiff