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