-
Soeren Peters authoredSoeren Peters authored
pdb.c 14.98 KiB
/************************************************************************/
/*! \file pdb.c
\brief Functions for parsing pdb files.
Pdb reader (parser). Loads arrays of pointers for easy backbone access.
\date Started 10/20/06
\author Kevin
\version $Id: pdb.c 10711 2011-08-31 22:23:04Z karypis $
*/
/************************************************************************/
#include <GKlib.h>
/************************************************************************/
/*! \brief Converts three-letter amino acid codes to one-leter codes.
This function takes a three letter \c * and converts it to a single \c
\param res is the three-letter code to be converted.
\returns A \c representing the amino acid.
*/
/************************************************************************/
char gk_threetoone(char *res) { /* {{{ */
/* make sure the matching works */
res[0] = toupper(res[0]);
res[1] = toupper(res[1]);
res[2] = toupper(res[2]);
if(strcmp(res,"ALA") == 0) {
return 'A';
}
else if(strcmp(res,"CYS") == 0) {
return 'C';
}
else if(strcmp(res,"ASP") == 0) {
return 'D';
}
else if(strcmp(res,"GLU") == 0) {
return 'E';
}
else if(strcmp(res,"PHE") == 0) {
return 'F';
}
else if(strcmp(res,"GLY") == 0) {
return 'G';
}
else if(strcmp(res,"HIS") == 0) {
return 'H';
}
else if(strcmp(res,"ILE") == 0) {
return 'I';
}
else if(strcmp(res,"LYS") == 0) {
return 'K';
}
else if(strcmp(res,"LEU") == 0) {
return 'L';
}
else if(strcmp(res,"MET") == 0) {
return 'M';
}
else if(strcmp(res,"ASN") == 0) {
return 'N';
}
else if(strcmp(res,"PRO") == 0) {
return 'P';
}
else if(strcmp(res,"GLN") == 0) {
return 'Q';
} else if(strcmp(res,"ARG") == 0) {
return 'R';
}
else if(strcmp(res,"SER") == 0) {
return 'S';
}
else if(strcmp(res,"THR") == 0) {
return 'T';
}
else if(strcmp(res,"SCY") == 0) {
return 'U';
}
else if(strcmp(res,"VAL") == 0) {
return 'V';
}
else if(strcmp(res,"TRP") == 0) {
return 'W';
}
else if(strcmp(res,"TYR") == 0) {
return 'Y';
}
else {
return 'X';
}
} /* }}} */
/************************************************************************/
/*! \brief Frees the memory of a pdbf structure.
This function takes a pdbf pointer and frees all the memory below it.
\param p is the pdbf structure to be freed.
*/
/************************************************************************/
void gk_freepdbf(pdbf *p) { /* {{{ */
int i;
if(p != NULL) {
gk_free((void **)&p->resSeq, LTERM);
for(i=0; i<p->natoms; i++) {
gk_free((void **)&p->atoms[i].name, &p->atoms[i].resname, LTERM);
}
for(i=0; i<p->nresidues; i++) {
gk_free((void *)&p->threeresSeq[i], LTERM);
}
/* this may look like it's wrong, but it's just a 1-d array of pointers, and
the pointers themselves are freed above */
gk_free((void **)&p->bbs, &p->cas, &p->atoms, &p->cm, &p->threeresSeq, LTERM);
}
gk_free((void **)&p, LTERM);
} /* }}} */
/************************************************************************/
/*! \brief Reads a pdb file into a pdbf structure
This function allocates a pdbf structure and reads the file fname into
that structure.
\param fname is the file name to be read
\returns A filled pdbf structure.
*/
/************************************************************************/
pdbf *gk_readpdbfile(char *fname) { /* {{{ */
int i=0, res=0;
char linetype[6];
int aserial;
char aname[5] = " \0";
char altLoc = ' ';
char rname[4] = " \0";
char chainid = ' ';
char oldchainid = ' '; int rserial;
int oldRserial = -37;
char icode = ' ';
char element = ' ';
double x;
double y;
double z;
double avgx;
double avgy;
double avgz;
double opcy;
double tmpt;
char line[MAXLINELEN];
int corruption=0;
int nresatoms;
int atoms=0, residues=0, cas=0, bbs=0, firstres=1;
pdbf *toFill = gk_malloc(sizeof(pdbf),"fillme");
FILE *FPIN;
FPIN = gk_fopen(fname,"r",fname);
while(fgets(line, 256, FPIN)) {
sscanf(line,"%s ",linetype);
/* It seems the only reliable parts are through temperature, so we only use these parts */
/* if(strstr(linetype, "ATOM") != NULL || strstr(linetype, "HETATM") != NULL) { */
if(strstr(linetype, "ATOM") != NULL) {
sscanf(line, "%6s%5d%*1c%4c%1c%3c%*1c%1c%4d%1c%*3c%8lf%8lf%8lf%6lf%6lf %c\n",
linetype,&aserial,aname,&altLoc,rname,&chainid,&rserial,&icode,&x,&y,&z,&opcy,&tmpt,&element);
sscanf(linetype, " %s ",linetype);
sscanf(aname, " %s ",aname);
sscanf(rname, " %s ",rname);
if(altLoc != ' ') {
corruption = corruption|CRP_ALTLOCS;
}
if(firstres == 1) {
oldRserial = rserial;
oldchainid = chainid;
residues++;
firstres = 0;
}
if(oldRserial != rserial) {
residues++;
oldRserial = rserial;
}
if(oldchainid != chainid) {
corruption = corruption|CRP_MULTICHAIN;
}
oldchainid = chainid;
atoms++;
if(strcmp(aname,"CA") == 0) {
cas++;
}
if(strcmp(aname,"N") == 0 || strcmp(aname,"CA") == 0 ||
strcmp(aname,"C") == 0 || strcmp(aname,"O") == 0) {
bbs++;
}
}
else if(strstr(linetype, "ENDMDL") != NULL || strstr(linetype, "END") != NULL || strstr(linetype, "TER") != NULL) {
break;
}
}
fclose(FPIN);
/* printf("File has coordinates for %d atoms in %d residues\n",atoms,residues); */
toFill->natoms = atoms;
toFill->ncas = cas;
toFill->nbbs = bbs;
toFill->nresidues = residues;
toFill->resSeq = (char *) gk_malloc (residues*sizeof(char),"residue seq"); toFill->threeresSeq = (char **)gk_malloc (residues*sizeof(char *),"residue seq");
toFill->atoms = (atom *) gk_malloc (atoms*sizeof(atom), "atoms");
toFill->bbs = (atom **)gk_malloc ( bbs*sizeof(atom *),"bbs");
toFill->cas = (atom **)gk_malloc ( cas*sizeof(atom *),"cas");
toFill->cm = (center_of_mass *)gk_malloc(residues*sizeof(center_of_mass),"center of mass");
res=0; firstres=1; cas=0; bbs=0; i=0;
avgx = 0.0; avgy = 0.0; avgz = 0.0;
nresatoms = 0;
FPIN = gk_fopen(fname,"r",fname);
while(fgets(line, 256, FPIN)) {
sscanf(line,"%s ",linetype);
/* It seems the only reliable parts are through temperature, so we only use these parts */
/* if(strstr(linetype, "ATOM") != NULL || strstr(linetype, "HETATM") != NULL) { */
if(strstr(linetype, "ATOM") != NULL ) {
/* to ensure our memory doesn't get corrupted by the biologists, we only read this far */
sscanf(line, "%6s%5d%*1c%4c%1c%3c%*1c%1c%4d%1c%*3c%8lf%8lf%8lf%6lf%6lf %c\n",
linetype,&aserial,aname,&altLoc,rname,&chainid,&rserial,&icode,&x,&y,&z,&opcy,&tmpt,&element);
sscanf(aname, "%s",aname);
sscanf(rname, "%s",rname);
if(firstres == 1) {
toFill->resSeq[res] = gk_threetoone(rname);
toFill->threeresSeq[res] = gk_strdup(rname);
oldRserial = rserial;
res++;
firstres = 0;
}
if(oldRserial != rserial) {
/* we're changing residues. store the center of mass from the last one & reset */
toFill->cm[res-1].x = avgx/nresatoms;
toFill->cm[res-1].y = avgy/nresatoms;
toFill->cm[res-1].z = avgz/nresatoms;
avgx = 0.0; avgy = 0.0; avgz = 0.0;
nresatoms = 0;
toFill->cm[res-1].name = toFill->resSeq[res-1];
toFill->threeresSeq[res] = gk_strdup(rname);
toFill->resSeq[res] = gk_threetoone(rname);
res++;
oldRserial = rserial;
}
avgx += x;
avgy += y;
avgz += z;
nresatoms++;
toFill->atoms[i].x = x;
toFill->atoms[i].y = y;
toFill->atoms[i].z = z;
toFill->atoms[i].opcy = opcy;
toFill->atoms[i].tmpt = tmpt;
toFill->atoms[i].element = element;
toFill->atoms[i].serial = aserial;
toFill->atoms[i].chainid = chainid;
toFill->atoms[i].altLoc = altLoc;
toFill->atoms[i].rserial = rserial;
toFill->atoms[i].icode = icode;
toFill->atoms[i].name = gk_strdup(aname);
toFill->atoms[i].resname = gk_strdup(rname);
/* Set up pointers for the backbone and c-alpha shortcuts */
if(strcmp(aname,"CA") == 0) {
toFill->cas[cas] = &(toFill->atoms[i]);
cas++;
}
if(strcmp(aname,"N") == 0 || strcmp(aname,"CA") == 0 || strcmp(aname,"C") == 0 || strcmp(aname,"O") == 0) {
toFill->bbs[bbs] = &(toFill->atoms[i]);
bbs++;
} i++;
}
else if(strstr(linetype, "ENDMDL") != NULL || strstr(linetype, "END") != NULL || strstr(linetype, "TER") != NULL) {
break;
}
}
/* get that last average */
toFill->cm[res-1].x = avgx/nresatoms;
toFill->cm[res-1].y = avgy/nresatoms;
toFill->cm[res-1].z = avgz/nresatoms;
/* Begin test code */
if(cas != residues) {
printf("Number of residues and CA coordinates differs by %d (!)\n",residues-cas);
if(cas < residues) {
corruption = corruption|CRP_MISSINGCA;
}
else if(cas > residues) {
corruption = corruption|CRP_MULTICA;
}
}
if(bbs < residues*4) {
corruption = corruption|CRP_MISSINGBB;
}
else if(bbs > residues*4) {
corruption = corruption|CRP_MULTIBB;
}
fclose(FPIN);
toFill->corruption = corruption;
/* if(corruption == 0)
printf("File was clean!\n"); */
return(toFill);
} /* }}} */
/************************************************************************/
/*! \brief Writes the sequence of residues from a pdb file.
This function takes a pdbf structure and a filename, and writes out
the amino acid sequence according to the atomic coordinates. The output
is in fasta format.
\param p is the pdbf structure with the sequence of interest
\param fname is the file name to be written
*/
/************************************************************************/
void gk_writefastafrompdb(pdbf *pb, char *fname) {
int i;
FILE *FPOUT;
FPOUT = gk_fopen(fname,"w",fname);
fprintf(FPOUT,"> %s\n",fname);
for(i=0; i<pb->nresidues; i++)
fprintf(FPOUT,"%c",pb->resSeq[i]);
fprintf(FPOUT,"\n");
fclose(FPOUT);
}
/************************************************************************/
/*! \brief Writes all centers of mass in pdb-format to file fname.
This function takes a pdbf structure and writes out the calculated
mass center information to file fname as though each one was a c-alpha.
\param p is the pdbf structure to write out
\param fname is the file name to be written
*/
/************************************************************************/
void gk_writecentersofmass(pdbf *p, char *fname) { int i;
FILE *FPIN;
FPIN = gk_fopen(fname,"w",fname);
for(i=0; i<p->nresidues; i++) {
fprintf(FPIN,"%-6s%5d %4s%1c%3s %1c%4d%1c %8.3lf%8.3lf%8.3lf%6.2f%6.2f\n",
"ATOM ",i,"CA",' ',p->threeresSeq[i],' ',i,' ',p->cm[i].x,p->cm[i].y,p->cm[i].z,1.0,-37.0);
}
fclose(FPIN);
}
/************************************************************************/
/*! \brief Writes all atoms in p in pdb-format to file fname.
This function takes a pdbf structure and writes out all the atom
information to file fname.
\param p is the pdbf structure to write out
\param fname is the file name to be written
*/
/************************************************************************/
void gk_writefullatom(pdbf *p, char *fname) {
int i;
FILE *FPIN;
FPIN = gk_fopen(fname,"w",fname);
for(i=0; i<p->natoms; i++) {
fprintf(FPIN,"%-6s%5d %4s%1c%3s %1c%4d%1c %8.3lf%8.3lf%8.3lf%6.2f%6.2f\n",
"ATOM ",p->atoms[i].serial,p->atoms[i].name,p->atoms[i].altLoc,p->atoms[i].resname,p->atoms[i].chainid,p->atoms[i].rserial,p->atoms[i].icode,p->atoms[i].x,p->atoms[i].y,p->atoms[i].z,p->atoms[i].opcy,p->atoms[i].tmpt);
}
fclose(FPIN);
}
/************************************************************************/
/*! \brief Writes out all the backbone atoms of a structure in pdb format
This function takes a pdbf structure p and writes only the backbone atoms
to a filename fname.
\param p is the pdb structure to write out.
\param fname is the file name to be written.
*/
/************************************************************************/
void gk_writebackbone(pdbf *p, char *fname) {
int i;
FILE *FPIN;
FPIN = gk_fopen(fname,"w",fname);
for(i=0; i<p->nbbs; i++) {
fprintf(FPIN,"%-6s%5d %4s%1c%3s %1c%4d%1c %8.3lf%8.3lf%8.3lf%6.2f%6.2f\n",
"ATOM ",p->bbs[i]->serial,p->bbs[i]->name,p->bbs[i]->altLoc,p->bbs[i]->resname,p->bbs[i]->chainid,p->bbs[i]->rserial,p->bbs[i]->icode,p->bbs[i]->x,p->bbs[i]->y,p->bbs[i]->z,p->bbs[i]->opcy,p->bbs[i]->tmpt);
}
fclose(FPIN);
}
/************************************************************************/
/*! \brief Writes out all the alpha carbon atoms of a structure
This function takes a pdbf structure p and writes only the alpha carbon
atoms to a filename fname.
\param p is the pdb structure to write out.
\param fname is the file name to be written.
*/
/************************************************************************/
void gk_writealphacarbons(pdbf *p, char *fname) {
int i;
FILE *FPIN;
FPIN = gk_fopen(fname,"w",fname);
for(i=0; i<p->ncas; i++) {
fprintf(FPIN,"%-6s%5d %4s%1c%3s %1c%4d%1c %8.3lf%8.3lf%8.3lf%6.2f%6.2f\n",
"ATOM ",p->cas[i]->serial,p->cas[i]->name,p->cas[i]->altLoc,p->cas[i]->resname,p->cas[i]->chainid,p->cas[i]->rserial,p->cas[i]->icode,p->cas[i]->x,p->cas[i]->y,p->cas[i]->z,p->cas[i]->opcy,p->cas[i]->tmpt);
} fclose(FPIN);
}
/************************************************************************/
/*! \brief Decodes the corruption bitswitch and prints any problems
Due to the totally unreliable nature of the pdb format, reading a pdb
file stores a corruption bitswitch, and this function decodes that switch
and prints the result on stdout.
\param p is the pdb structure to write out.
\param fname is the file name to be written.
*/
/************************************************************************/
void gk_showcorruption(pdbf *p) {
int corruption = p->corruption;
if(corruption&CRP_ALTLOCS)
printf("Multiple coordinate sets for at least one atom\n");
if(corruption&CRP_MISSINGCA)
printf("Missing coordiantes for at least one CA atom\n");
if(corruption&CRP_MISSINGBB)
printf("Missing coordiantes for at least one backbone atom (N,CA,C,O)\n");
if(corruption&CRP_MULTICHAIN)
printf("File contains coordinates for multiple chains\n");
if(corruption&CRP_MULTICA)
printf("Multiple CA atoms found for the same residue (could be alternate locators)\n");
if(corruption&CRP_MULTICA)
printf("Multiple copies of backbone atoms found for the same residue (could be alternate locators)\n");
}
/* sscanf(line, "%6s%5d%*1c%4s%1c%3s%*1c%1c%4d%1c%*3c%8lf%8lf%8lf%6lf%6lf%*6c%4s%2s%2s\n",
linetype,&aserial,aname,&altLoc,rname,&chainid,&rserial,&icode,&x,&y,&z,&opcy,&tmpt,segId,element,charge);
printf(".%s.%s.%s.\n",segId,element,charge);
printf("%-6s%5d%-1s%-4s%1c%3s%1s%1c%4d%1c%3s%8.3lf%8.3lf%8.3lf%6.2f%6.2f%6s%4s%2s%2s\n",
linetype,aserial," ",aname,altLoc,rname," ",chainid,rserial,icode," ",x,y,z,opcy,tmpt," ",segId,element,charge); */
/* and we could probably get away with this using astral files, */
/* sscanf(line, "%6s%5d%*1c%4s%1c%3s%*1c%1c%4d%1c%*3c%8lf%8lf%8lf%6lf%6lf%*6c%6s\n",
linetype,&aserial,aname,&altLoc,rname,&chainid,&rserial,&icode,&x,&y,&z,&opcy,&tmpt,element);
printf("%-6s%5d%-1s%-4s%1c%3s%1s%1c%4d%1c%3s%8.3lf%8.3lf%8.3lf%6.2f%6.2f%6s%6s\n",
linetype,aserial," ",aname,altLoc,rname," ",chainid,rserial,icode," ",x,y,z,opcy,tmpt," ",element); */