home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #3 / NN_1993_3.iso / spool / bionet / software / 2532 < prev    next >
Encoding:
Text File  |  1993-01-27  |  12.2 KB  |  417 lines
  1. Newsgroups: bionet.software
  2. Path: sparky!uunet!stanford.edu!ames!haven.umd.edu!darwin.sura.net!gatech!concert!uvaarpa!murdoch!darwin.clas.Virginia.EDU!rbm8p
  3. From: rbm8p@darwin.clas.Virginia.EDU (Rick MacDonald)
  4. Subject: Phi-Psi angles from Brookhaven Protein Data (C code)
  5. Message-ID: <1993Jan26.174853.12244@murdoch.acc.Virginia.EDU>
  6. Followup-To: bionetnet.software
  7. Originator: rbm8p@darwin.clas.Virginia.EDU
  8. Keywords: protein conformation software
  9. Sender: usenet@murdoch.acc.Virginia.EDU
  10. Organization: University of Virginia
  11. Date: Tue, 26 Jan 1993 17:48:53 GMT
  12. Lines: 403
  13.  
  14.  
  15.  
  16. George Johnson <geojohn@casbah.acns.nwu.edu> recently asked for a
  17. routine to calculate phi and psi angles from Brookhaven (pdb) 
  18. coordinates.   
  19.  
  20. The short C program that follows calculates phi-psi and alpha-carbon
  21. bend and torsion angles for Brookhaven protein data files.  The
  22. people who manage the data bank also have FORTRAN routines to do
  23. this.
  24.  
  25. I have run it on SGI 4D and IBM RS6000's, but you may have to fiddle
  26. with your local math library routines.
  27.  
  28. To Build: cc sifi.c -o sifi -lm  
  29.  
  30. Usage: sifi <pdbfile>
  31.  
  32. /*********************   CUT HERE      ********************************/
  33.  
  34. /***********************************************************************
  35. This program was written by 
  36.  
  37. Rick MacDonald
  38. University of Virginia Biophysics
  39.  
  40. rbm8p@virginia.edu
  41.  
  42. This program will print for each amino acid residue in a pdb file: 
  43.     psi
  44.     phi
  45.     alpha carbon bend
  46.     alpha carbon dihedral
  47.  
  48. To Build: cc sifi.c -o sifi -lm  
  49.  
  50. Usage: sifi <pdbfile>
  51. ***********************************************************************/
  52. #include <stdio.h>
  53. #include <math.h>
  54.  
  55. #define TRUE    1
  56. #define FALSE    0
  57.  
  58. #define LINELEN 100
  59.  
  60. #define MAXATOMS 10000
  61. #define MAXRES   1000
  62.  
  63. typedef struct {
  64.     int    num;
  65.     char    a_type[6];
  66.     char    r_type[4];
  67.     char    subunit;
  68.     int    resnum;
  69.     float    p[3];        /* position x,y,z */
  70. } ATOM;
  71.  
  72. typedef struct {
  73.     int    nitrog;
  74.     int    alphac;
  75.     int    carbon;
  76.     int    oxygen;
  77.     char    *r_type;    /* 3 letter code */
  78.     char    res_type;   /* 1 letter code */
  79.     char    subunit;
  80.     float    phi;
  81.     float    psi;
  82.     float    cab;
  83.     float    cad;
  84.     int     resnum;
  85. } AMINO;
  86.  
  87. ATOM    atom[MAXATOMS];
  88. AMINO    res[MAXRES];
  89.  
  90. int     err;
  91.  
  92.  
  93. main(int argc, char **argv)
  94. {
  95. char    *strstr(), *strncpy(), *strcat();
  96. char    *next_str(char *string);
  97. float    angle(float *v1, float *v2);
  98. float    dot(float *v1, float *v2);
  99. float    dihedral(float *p1, float *p2, float *p3, float *p4);
  100. void    cross(float *v1, float *v2, float *p);
  101. float    len(float *v1);
  102. int     cmpstr(char *s1, char *s2);
  103. void    read_atom_record(char *rec, ATOM *atom1);
  104. char    single_letter_code(char *string); 
  105.  
  106.     register i;
  107.  
  108.     FILE    *Fpdb;
  109.  
  110.     char    tname[80], fname[80], line[LINELEN], *temp, WAS, WUZ, WLB;
  111.  
  112.     float  v1[3], v2[3];
  113.     int        ndx, n_atoms, n_res, rno, is, was;
  114.     
  115. /*-----------------------------------------------------------------------
  116. Open file specified by command line, if it exists 
  117. -----------------------------------------------------------------------*/
  118.     if(argc<2) err = printf("Usage: sifi <pdbfile>\n");
  119.  
  120.     if(argc!=2) exit(1);
  121.  
  122.     temp = strncpy(fname,argv[argc-1],sizeof(tname));
  123.     if ((Fpdb = fopen(fname,"r"))==NULL)
  124.         err = printf("Could not open file %s\n",fname);
  125.  
  126. /*-----------------------------------------------------------------------
  127. Read through the PDB format file to extract ATOM records
  128. -----------------------------------------------------------------------*/
  129.  
  130.     ndx = 0;
  131.     while (fgets(line,LINELEN,Fpdb)!=0) {
  132.         if (strstr(line,"ATOM  ")!=NULL) {
  133.             read_atom_record(line,&atom[ndx]);
  134.             ndx++;
  135.         }
  136.     }
  137.     n_atoms = ndx;
  138.  
  139. /*-----------------------------------------------------------------------
  140. Now we've got the ATOMs. Re-order them into residues. 
  141. -----------------------------------------------------------------------*/
  142.     rno = 0;
  143.     was = atom[0].resnum;
  144.     for (ndx = 0; ndx < n_atoms; ndx++) {
  145.     is = atom[ndx].resnum;
  146.         if (is != was) rno++;
  147.     if (cmpstr(atom[ndx].a_type,"CA")!=NULL) {
  148.         res[rno].resnum = atom[ndx].resnum;
  149.         res[rno].alphac = ndx;
  150.         res[rno].r_type = atom[ndx].r_type;
  151.         res[rno].subunit= atom[ndx].subunit;
  152.     }
  153.     else if (cmpstr(atom[ndx].a_type,"N")!=NULL)
  154.         res[rno].nitrog = ndx;
  155.     else if (cmpstr(atom[ndx].a_type,"C")!=NULL)
  156.         res[rno].carbon = ndx;
  157.     else if (cmpstr(atom[ndx].a_type,"O")!=NULL)
  158.         res[rno].oxygen = ndx;
  159.     was = is;
  160.     }
  161.     
  162.     n_res = rno; 
  163.  
  164. /*-----------------------------------------------------------------------
  165. Find for each residue:
  166.         single letter residue code
  167.     phi angle
  168.     psi angle
  169.     cab angle = bend angle between adjacent alpha carbons 
  170.     cad angle = is the dihedral() or torsion angle
  171. -----------------------------------------------------------------------*/
  172.  
  173.     res[0].res_type = single_letter_code(res[0].r_type);
  174.     res[0].psi        = -(dihedral(atom[res[1].nitrog].p,
  175.                          atom[res[1].alphac].p,
  176.                          atom[res[1].carbon].p,
  177.                          atom[res[2].nitrog].p));
  178.     res[1].res_type = single_letter_code(res[1].r_type);
  179.     res[1].phi        = -(dihedral(atom[res[0].carbon].p,
  180.                          atom[res[1].nitrog].p,
  181.                          atom[res[1].alphac].p,
  182.                          atom[res[1].carbon].p));
  183.     res[1].psi        = -(dihedral(atom[res[1].nitrog].p,
  184.                          atom[res[1].alphac].p,
  185.                          atom[res[1].carbon].p,
  186.                          atom[res[2].nitrog].p));
  187.  
  188.     for (rno = 2; rno <= n_res; rno++) {
  189.     res[rno].res_type   = single_letter_code(res[rno].r_type);
  190.     
  191.     res[rno].phi        = -(dihedral(atom[res[rno-1].carbon].p,
  192.                          atom[res[rno  ].nitrog].p,
  193.                          atom[res[rno  ].alphac].p,
  194.                          atom[res[rno  ].carbon].p));
  195.                     
  196.     res[rno].psi        = -(dihedral(atom[res[rno  ].nitrog].p,
  197.                          atom[res[rno  ].alphac].p,
  198.                          atom[res[rno  ].carbon].p,
  199.                          atom[res[rno+1].nitrog].p));
  200.                     
  201.         res[rno].cad        = dihedral( atom[res[rno-2].alphac].p,
  202.                                         atom[res[rno-1].alphac].p,
  203.                         atom[res[rno  ].alphac].p,
  204.                         atom[res[rno+1].alphac].p);
  205.         for (i=0;i<=2;i++) {
  206.             v1[i] = atom[res[rno-1].alphac].p[i] - atom[res[rno].alphac].p[i];
  207.             v2[i] = atom[res[rno+1].alphac].p[i] - atom[res[rno].alphac].p[i];
  208.         }    
  209.         res[rno].cab        = angle(v1, v2);
  210.     }
  211.             
  212.  
  213. /*-----------------------------------------------------------------------
  214. Print this stuff out as a table.
  215. -----------------------------------------------------------------------*/
  216.         printf("%s\n", fname);
  217.         printf("   #     resid   phi  psi    Cab  Cat\n");
  218.  
  219.         WAS = 'x';
  220.     WUZ = 'x';
  221.     WLB = res[0].subunit;
  222.     res[n_res+1].subunit = 'x';
  223.         for (rno = 0; rno <= n_res; rno++) {
  224.     
  225.         if (res[rno].subunit != WAS) {      /* This gibberish flags the meaningless      */
  226.         res[rno].phi = 999.;            /* values at the ends of each subunit as 999 */
  227.         }
  228.         if (res[rno].subunit != WUZ) {
  229.         res[rno].cab = 999.;
  230.         res[rno].cad = 999.;
  231.         }
  232.         WUZ = WAS;
  233.         WAS = res[rno].subunit;
  234.         WLB = res[rno+1].subunit;
  235.         if (res[rno].subunit != WLB) {
  236.         res[rno].psi = 999.;
  237.         res[rno].cab = 999.;
  238.         res[rno].cad = 999.;
  239.         }
  240.         
  241.     
  242.         printf("%4d %c   %3s %c ", 
  243.             res[rno].resnum, res[rno].subunit, res[rno].r_type, res[rno].res_type);
  244.         printf("%5.0f", res[rno].phi);
  245.         printf("%5.0f", res[rno].psi);
  246.         printf("  %5.0f", res[rno].cab);
  247.         printf("%5.0f", res[rno].cad);
  248.         printf("\n");
  249.         }
  250.  
  251. } /* end main */
  252.  
  253. /***********************************************************************
  254.  
  255. ***********************************************************************/
  256. void read_atom_record(char *rec, ATOM *atom1) 
  257. {
  258.     char    *next_str(char *string);
  259.     char    *ptr;
  260.  
  261.     float    fl;
  262.  
  263.     ptr = rec;
  264.     
  265.     err = sscanf( (ptr = next_str(ptr)), "%d",  &atom1->num    );
  266.     err = sscanf( (ptr = next_str(ptr)), "%3s",  atom1->a_type );
  267.     err = sscanf( (ptr = next_str(ptr)), "%3s",   atom1->r_type );
  268.     ptr = ptr + 4;
  269.     atom1->subunit = *ptr;
  270.     err = sscanf( (ptr = next_str(ptr)), "%d",  &atom1->resnum );
  271.  
  272.  
  273. /*  Remove this comment to print atom records as read from pdb for debug
  274. printf("\n%5d %4s %s %c %d ", atom1->num, atom1->a_type, atom1->r_type, atom1->subunit, atom1->resnum);
  275. */
  276.  
  277.  
  278. /* This is a silly patch --- but I can't make %Lf read a float on the SGI! */
  279.     err = sscanf( (ptr = next_str(ptr)), "%f", &fl   );
  280.         atom1->p[0] = fl;
  281.     err = sscanf( (ptr = next_str(ptr)), "%f", &fl   );
  282.         atom1->p[1] = fl;
  283.     err = sscanf( (ptr = next_str(ptr)), "%f", &fl   );
  284.         atom1->p[2] = fl;
  285.  
  286.         return;
  287. }
  288.  
  289. /***********************************************************************
  290.  
  291. ***********************************************************************/
  292. char single_letter_code(char *string) 
  293. {
  294.              if (cmpstr(string,"CYS")!=NULL) return 'C';
  295.     else if (cmpstr(string,"HIS")!=NULL) return 'H';
  296.     else if (cmpstr(string,"ILE")!=NULL) return 'I';
  297.     else if (cmpstr(string,"MET")!=NULL) return 'M';
  298.     else if (cmpstr(string,"SER")!=NULL) return 'S';
  299.     else if (cmpstr(string,"VAL")!=NULL) return 'V';
  300.     else if (cmpstr(string,"ALA")!=NULL) return 'A';
  301.     else if (cmpstr(string,"GLY")!=NULL) return 'G';
  302.     else if (cmpstr(string,"LEU")!=NULL) return 'L';
  303.     else if (cmpstr(string,"PRO")!=NULL) return 'P';
  304.     else if (cmpstr(string,"THR")!=NULL) return 'T';
  305.     else if (cmpstr(string,"PHE")!=NULL) return 'F';
  306.     else if (cmpstr(string,"ARG")!=NULL) return 'R';
  307.     else if (cmpstr(string,"TYR")!=NULL) return 'Y';
  308.     else if (cmpstr(string,"TRP")!=NULL) return 'W';
  309.     else if (cmpstr(string,"ASN")!=NULL) return 'N';
  310.     else if (cmpstr(string,"GLU")!=NULL) return 'E';
  311.     else if (cmpstr(string,"GLN")!=NULL) return 'Q';
  312.     else if (cmpstr(string,"LYS")!=NULL) return 'K';
  313.     else if (cmpstr(string,"ASP")!=NULL) return 'D';
  314.     return '*';
  315. }
  316.  
  317. /***********************************************************************
  318.  cmpstr does a simple compare of a string s1 with string s2 and returns a  
  319.  pointer to the beginning of s1 if same or NULL if not same
  320. ***********************************************************************/
  321. int cmpstr(char *s1, char *s2)
  322. {
  323.     loop:  if(*s1 != *s2) return(0);
  324.     if(*s1 == '\0') return(1);
  325.     s1++;   s2++;
  326.     goto loop;
  327. }
  328.  
  329. /***********************************************************************
  330.  
  331. ***********************************************************************/
  332. char *next_str(char *string)
  333. {
  334.     char    c;
  335.     while (((c = *string) != ' ') && (c != '\t')) string++;
  336.     while (((c = *string) == ' ') || (c == '\t')) string++;
  337.     return string;
  338. }
  339.  
  340. /***********************************************************************
  341. Calculate the dihedral angle between two vectors given endpoints.
  342. ***********************************************************************/
  343. float  dihedral(float *p1, float *p2, float *p3, float *p4)
  344. {
  345. float  dot(float *v1, float *v2);
  346. void    cross(float *v1, float *v2, float *p);
  347.  
  348.     float    v1[3], v2[3], v3[3], p[3], q[3], r[3], theta;
  349.     register i;
  350.  
  351.         for (i = 0; (i < 3); i++) {
  352.         v1[i] = (p2[i] - p1[i]);
  353.         v2[i] = (p2[i] - p3[i]);
  354.         v3[i] = (p4[i] - p3[i]);
  355.     }
  356.     cross(v1,v2,p);
  357.     cross(v2,v3,q);
  358.     theta = angle(p,q);
  359.     cross(p,q,r);
  360.     if(dot(v2,r) < 0) theta = -theta;
  361.     return theta;
  362. }
  363.  
  364. /***********************************************************************
  365. The vector cross product is returned as the third parameter
  366. ***********************************************************************/
  367. void  cross(float *v1, float *v2, float *p)
  368. {
  369.         p[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
  370.         p[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
  371.         p[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
  372. }
  373.  
  374. /***********************************************************************
  375. Return angle (degrees) between two vectors
  376. ***********************************************************************/
  377. float  angle(float *v1, float *v2)
  378. {
  379.     
  380. float    len(float *v1);
  381. float    dot(float *v1, float *v2);
  382.     float cos_val, sin_val, sinsq;
  383.  
  384.     cos_val = dot(v1,v2)/(len(v1)*len(v2));
  385.     sinsq = 1.0 - cos_val*cos_val;
  386.     if( sinsq < 0.0 ) sinsq = 0.0;
  387.     sin_val = sqrt(sinsq);
  388.     return(atan2(sin_val,cos_val) * 57.29578); /* pi radians */
  389.  
  390. }
  391.  
  392. /***********************************************************************
  393. Return vector dot product
  394. ***********************************************************************/
  395. float  dot(float *v1, float *v2)
  396. {
  397.     int    i;
  398.     float    sum;
  399.  
  400.     sum = 0;
  401.     for (i = 0; (i <= 2); i++)  {
  402.         sum = sum + (v1[i] * v2[i]);
  403.     }
  404.  
  405.     return sum;
  406. }
  407.  
  408. /***********************************************************************
  409. Return the length of a vector
  410. ***********************************************************************/
  411. float len (float *v1)
  412. {
  413.     return ( fexp(0.5*flog ( (v1[0]*v1[0]) + (v1[1]*v1[1]) + (v1[2]*v1[2]) )) );
  414. }
  415.  
  416.  
  417.