protein.c

药物开发中的基于结构的从头设计代码

# include "pocket.h"

Protein::Protein()
{
	num_res=0;

	res=new Residue[50];
	if(res==NULL) Memory_Allocation_Error();

	num_atom=0;

	atom=new Patom[10000];
	if(atom==NULL) Memory_Allocation_Error();

	num_pocket_res=0;
	
	pocket_res=new Pocket_res[50];
	if(pocket_res==NULL) Memory_Allocation_Error();
	
}

Protein::~Protein()
{
	delete [] res;
	delete [] atom;
	delete [] pocket_res;
}

void Protein::Read_RESIDUE(char *dirname)
{
	FILE *fp;
	int i,j,num;
	char buf[160],head[6],filename[160];

	strcpy(filename,dirname); strcat(filename,"RESIDUE");
	if((fp=fopen(filename,"r"))==NULL) Openning_File_Error(filename);

	i=0;
	for(;;)
		{
		 if(fgets(buf,160,fp)==NULL) break;
		 else if(Blank_Line_Check(buf)==TRUE) continue;
		 else if(buf[0]=='#') continue;
		 
		 sscanf(buf,"%s",head);

		 if(!strcmp(head,"END")) break;
		 else if(!strcmp(head,"RESI"))
			{
			 num=0;
			 sscanf(buf,"%*s%s%d",res[i].name,&num);
			 res[i].num_atom=num;

			 for(j=0;j<num;j++)
				{
				 fgets(buf,160,fp);
				 sscanf(buf,"%s",head);
				 if(strcmp(head,"ATOM")) continue;
				 sscanf(buf,"%*s%s%*s%f%s%f%f%s",
					res[i].atom[j].name,
				       &res[i].atom[j].logp,
					res[i].atom[j].hb,
                                       &res[i].atom[j].q,
				       &res[i].atom[j].r,
					res[i].atom[j].type);
				}
			 i++;
			}
		}

	fclose(fp);

	num_res=i;

	return;
}

void Protein::Show_RESIDUE() const
{
	int i,j;

	printf("The number of residue templates is: %d\n",num_res);

	for(i=0;i<num_res;i++)
		{
		 printf("Residue\t%3s\t%d\n",res[i].name,res[i].num_atom);
		 for(j=0;j<res[i].num_atom;j++)
			{
			 printf("ATOM %8s %6.3f %4s %6.3f %6.3f %8s\n",
				res[i].atom[j].name,
				res[i].atom[j].logp,
				res[i].atom[j].hb,
				res[i].atom[j].q,
				res[i].atom[j].r,
				res[i].atom[j].type);
			}
		 printf("\n");
		}

	return;
}

void Protein::Read_PDB(char *filename)
{
	int i,j,num;
	FILE *fp;
	char buf[160],head[10],name[10];

	// Now read the information of receptor atoms from the PDB file 

	if((fp=fopen(filename,"r"))==NULL) Openning_File_Error(filename);

	num=0;

	for(;;)
		{
		 if(fgets(buf,160,fp)==NULL) break;

		 sscanf(buf,"%s",head);

		 if(!strcmp(head,"END")) break;
		 else if(!strcmp(head,"ATOM"))
			{
			 // Get the atom name and erase that section 
			 i=0;
			 for(j=12;j<17;j++)
				{
				 if(buf[j]==' ') continue;
				 else {name[i]=buf[j]; buf[j]=' ';i++;}
				}
			 name[i]='\0'; 

			 // Skip lone pairs and hydrogens  
			 if(name[0]=='L'&&name[1]=='P') continue;
			 else if(name[0]=='H') continue;
			 else if((name[0]=='1'||name[0]=='2'||name[0]=='3')
				 &&name[1]=='H') continue;

			 strcpy(atom[num].name,name);

			 // Get the chain label and erase that section 
			 if(buf[21]!=' ') 
                                {
                                 atom[num].chain=buf[21]; 
				 buf[21]=' ';
                                }
			 else atom[num].chain='A';	// default

			 sscanf(buf,"%*s%d%s%s%f%f%f",
				&atom[num].atom_id,
				 atom[num].residue,
				 atom[num].res_id,
				&atom[num].coor[0],
				&atom[num].coor[1],
				&atom[num].coor[2]);
			 atom[num].valid=1;
			 num++;
			}
		 else if(!strcmp(head,"HETATM"))
			{
			 sscanf(buf,"%*s%*d%s",head);

			 // skip hydrogens
			 if(head[0]=='H') continue; 
			 else if(head[0]=='1') continue;
			 else if(head[0]=='2') continue;
			 else if(head[0]=='3') continue;

			 strcpy(atom[num].residue,"HET");

			 buf[21]=' ';	// erase the chain label

			 sscanf(buf,"%*s%d%s%*s%s%f%f%f",
				&atom[num].atom_id,
				 atom[num].name,
				 atom[num].res_id,
				&atom[num].coor[0],
				&atom[num].coor[1],
				&atom[num].coor[2]);
			 atom[num].valid=1;
			 num++;
			}
		 else continue;
		}

	fclose(fp);

	num_atom=num;

	if(num_atom==0) PDB_Format_Error(filename);

	return;
}

void Protein::Show_PDB() const
{
	int i;

	printf("Total number of atoms in the receptor is %d\n",num_atom);

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid==0) continue;
		 printf("ATOM %5d %4s %3s %c %5s ",
			atom[i].atom_id, 
			atom[i].name,
			atom[i].residue, 
			atom[i].chain, 
			atom[i].res_id); 
		 printf("%7.2f %7.2f %7.2f\n",
			atom[i].coor[0], 
			atom[i].coor[1], 
			atom[i].coor[2]); 
		}

	return;
}

void Protein::Value_Atom()
{
	int i,j;
	int mark,tmp;

	for(i=0;i<num_atom;i++)
	{
	 // First, check the atom to see whether it has a legal residue

	 mark=FALSE;
	 for(j=0;j<num_res;j++)
		{
		 if(strcmp(atom[i].residue,res[j].name)) continue;
		 mark=TRUE; break;
		}

	 if(mark==FALSE)
		{
		 printf("Warning: cannot find parameter for the protein atom %d %s %s\n",
			atom[i].atom_id, atom[i].name, atom[i].residue);
		 printf("This atom is skipped.\n");
		 atom[i].valid=0; continue;
		}
	 else tmp=j;

	// Second, check the atom in the specific residue

	 mark=FALSE;
	 for(j=0;j<res[tmp].num_atom;j++)
		{
		 if(strcmp(atom[i].name,res[tmp].atom[j].name)) continue;
		 else
			{
			 strcpy(atom[i].type,res[tmp].atom[j].type);
			 atom[i].r=res[tmp].atom[j].r;
			 atom[i].q=res[tmp].atom[j].q;
			 atom[i].logp=res[tmp].atom[j].logp;
			 strcpy(atom[i].hb,res[tmp].atom[j].hb);
			 atom[i].valid=1;
			 mark=TRUE; break; 
			}
		}

	 if(mark==FALSE)
		{
		 printf("Warning: cannot find parameter for the protein atom %d %s %s\n",
                       	atom[i].atom_id,atom[i].name,atom[i].residue);
                 printf("This atom is skipped.\n");
                 atom[i].valid=0;
		}
	}

	// again, skip the hydrogens

	for(i=0;i<num_atom;i++)
		{
		 if(!strcmp(atom[i].type,"H")) atom[i].valid=0;
		}
	
	return;
}

void Protein::Cal_HB_Root()
{
	int i,j;
	int num;
	float tmpx,tmpy,tmpz,tmp,cutoff;
	int neib[50];

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid<=0) continue;
		 else if(!strcmp(atom[i].hb,"N")) continue;
		 else if(!strcmp(atom[i].hb,"M")) cutoff=9.00; // metal
		 else if(!strcmp(atom[i].type,"O.w")) cutoff=12.50; // water
		 else cutoff=4.00; // normal

		 for(j=0;j<50;j++) neib[j]=0;

		 num=0;

		 for(j=0;j<num_atom;j++)
			{
			 if(j==i) continue;
			 else if(atom[j].valid<=0) continue;

			 tmp=Distance2(atom[i].coor,atom[j].coor);

			 if(tmp>cutoff) continue;
			 else {neib[num]=j; num++;}
			}

		 if(num==0) atom[i].valid=0;	// no neighbor!
		 else
			{
			 tmpx=tmpy=tmpz=0.000;
			 for(j=0;j<num;j++)
			{
			 tmpx+=(atom[neib[j]].coor[0]/num);
			 tmpy+=(atom[neib[j]].coor[1]/num);
			 tmpz+=(atom[neib[j]].coor[2]/num);
			}
			 atom[i].root[0]=tmpx;
			 atom[i].root[1]=tmpy;
			 atom[i].root[2]=tmpz;
			}
		}

	return;
}

void Protein::Check_Metal()
{
	int i,j;
	float d=0.000;

/* For any protein atom within a certain range around a metal ion,
   its charge and HB capacity are switched off */

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid<=0) continue;
		 else if(strcmp(atom[i].hb,"M")) continue;
/*
		 printf("Find metal %d %s %8.3f %8.3f %8.3f\n",
			 atom[i].atom_id,
			 atom[i].name,
			 atom[i].coor[0],
			 atom[i].coor[1],
			 atom[i].coor[2]);
*/
                 for(j=0;j<num_atom;j++)
                        {
                         if(i==j) continue;
                         else if(atom[j].valid<=0) continue;
			 else if(!strcmp(atom[j].hb,"N")) continue;
			 else if(!strcmp(atom[j].hb,"M")) continue;

                         d=Distance2(atom[i].coor,atom[j].coor);

			 if(d<9.00)	// d<3.00A 
				{
				 strcpy(atom[j].hb,"N");
				 atom[j].q=0.000;
				}
                        }
                }

        return;
}

void Protein::Check_Water()
{
	int i,j;
	int count,contact;
	float d;

	count=0;

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid<=0) continue;
		 else if(strcmp(atom[i].type,"O.w")) continue;

		 contact=0;

		 for(j=0;j<num_atom;j++)
			{
			 if(i==j) continue;
			 else if(atom[j].valid<=0) continue;
			 else if(!strcmp(atom[j].type,"O.w")) continue;
	 
			 d=Distance2(atom[i].coor,atom[j].coor);

			 if(d<16.00) contact++;	// d<4.00A
			}
		
		 if(contact<=1) {atom[i].valid=0;continue;}
		 else {atom[i].valid=1; count++; continue;}
		}

	// printf("Totally %d water are reserved in the protein.\n", count);

	return;
}

void Protein::Show_Atom() const
{
	int i;

	printf("\nTotal number of atoms in the protein is %d\n",num_atom);

        for(i=0;i<num_atom;i++)
                {
                 printf("%5d ", atom[i].atom_id);
		 printf("%1d ", atom[i].valid);
		 printf("%4s ", atom[i].name);
		 printf("%5s ", atom[i].type);
		 printf("%3s ", atom[i].residue);
		 printf("%5s ", atom[i].res_id);
		 printf("%c ", atom[i].chain);
		 printf("%7.2f ", atom[i].coor[0]);
		 printf("%7.2f ", atom[i].coor[1]);
		 printf("%7.2f ", atom[i].coor[2]);
		 printf("%4.2f ", atom[i].r);
		 printf("%5.2f ", atom[i].q);
		 printf("%5.2f ", atom[i].logp);
		 printf("%2s ", atom[i].hb);
		 printf("%7.2f ", atom[i].root[0]);
                 printf("%7.2f ", atom[i].root[1]);
                 printf("%7.2f", atom[i].root[2]);
		 printf("\n");                
		}

	return;
}

void Protein::Find_Pocket_Atom()
{
	int i,j,mark,tmp,count;
	float d,cutoff;
	extern Ligand *lig; 

	// find pocket atoms on the protein and label their valid as '2'

	cutoff=(DIST_CUTOFF)*(DIST_CUTOFF);

	for(i=0;i<num_atom;i++)
	{
	 if(atom[i].valid<=0) continue;

	 // check whether this patom is close to the ligand

	 mark=FALSE;

	 for(j=0;j<lig->num_atom;j++)
		{
		 if(lig->atom[j].valid<=0) continue;
		 else
			{
			 d=Distance2(atom[i].coor,lig->atom[j].coor);
			 if(d>cutoff) continue;
                         else {mark=TRUE; break;}
			}
		}

	 if(mark==TRUE) 
		{
		 atom[i].valid=2;	// pocket atom

		 // check if it has been included in an known pocket residue	

		 tmp=FALSE;

		 for(j=0;j<num_pocket_res;j++)
        		{
         		 if(!strcmp(atom[i].residue,pocket_res[j].name) 
			  &&!strcmp(atom[i].res_id,pocket_res[j].id)
			  &&(atom[i].chain==pocket_res[j].chain))
				{tmp=TRUE;break;}
			 else continue;
			}

		 // if not, add the newly found pocket residue to the list

		 if((tmp==FALSE)&&(strcmp(atom[i].residue,"HET")))
			{
			 j=num_pocket_res;
			 strcpy(pocket_res[j].name,atom[i].residue);
			 strcpy(pocket_res[j].id,atom[i].res_id);
			 pocket_res[j].chain=atom[i].chain;
			 pocket_res[j].valid=1;
			 num_pocket_res++;
			}
		}
	}

	// define all the left atoms in pocket residues as pocket atoms 

	count=0;

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid<=0) continue;
		 else if(atom[i].valid==2) {count++; continue;}

	 	 for(j=0;j<num_pocket_res;j++)
		{
		 if(strcmp(atom[i].residue,pocket_res[j].name)) continue;
	 	 else if(strcmp(atom[i].res_id,pocket_res[j].id)) continue;
	 	 else if(atom[i].chain!=pocket_res[j].chain) continue;
	 	 else {atom[i].valid=2; count++; break;}
		}
		}

	if(count<12) // there should be at least 3 residues as pocket!
		{
		 printf("Error: cannot find enough pocket residues.\n");
		 printf("Probably the ligand is not complexed with the protein.\n");
		 printf("Please check it and try again.\n");
		 exit(1);
		}
		 	
	return;
}

void Protein::Show_Pocket_Residue() const
{
	int i;

	printf("The binding pocket is formed by:\n\n");

	for(i=0;i<num_pocket_res;i++)
                {
                 printf("Pocket residue %d: %s %s %c\n", 
			 i+1,
                         pocket_res[i].name, 
			 pocket_res[i].id,
                         pocket_res[i].chain); 
                }

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid!=2) continue;
		 else if(strcmp(atom[i].residue,"HET")) continue;
		 else printf("%d %s %s\n", 
			     atom[i].atom_id, 
			     atom[i].residue, 
			     atom[i].type);
		}

	return;
}

void Protein::Output_Pocket_Atom(char *filename) const
{
	FILE *fp;
	int i,num;

	if((fp=fopen(filename,"w"))==NULL) Openning_File_Error(filename);

	fprintf(fp,"#\n");
	fprintf(fp,"# This file is create by LigBuilder/Pocket.\n");
	fprintf(fp,"# It stores the parameters for the pocket atoms.\n");
	fprintf(fp,"# Creation time: %s", Get_Time());
	fprintf(fp,"#\n");

	fprintf(fp,"# 1st column: ID of the atom\n");
	fprintf(fp,"# 2nd column: ID of the atom in the original PDB file\n");
	fprintf(fp,"# 3rd column: atom type\n");
	fprintf(fp,"# 4th column: X coordinate\n");
	fprintf(fp,"# 5th column: Y coordinate\n");
	fprintf(fp,"# 6th column: Z coordinate\n");
	fprintf(fp,"# 7th column: vdw radius\n");
	fprintf(fp,"# 8th column: atomic hydrophobic scale\n");
	fprintf(fp,"# 9th column: H-bond character\n");
	fprintf(fp,"# 10th column: H-bond root X coordinate\n"); 
	fprintf(fp,"# 11th column: H-bond root Y coordinate\n"); 
	fprintf(fp,"# 12th column: H-bond root Z coordinate\n");
	fprintf(fp,"#\n");

	num=0;

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid!=2) continue;
		 else num++;
		}

	fprintf(fp,"<Start>\n");
	fprintf(fp,"<Atom_number> %d\n", num);

	num=1;

	for(i=0;i<num_atom;i++)
		{
		 if(atom[i].valid!=2) continue;
	
		 fprintf(fp,"%3d ", num);
		 fprintf(fp,"%5d ", atom[i].atom_id);
		 fprintf(fp,"%5s ", atom[i].type);
		 fprintf(fp,"%7.2f ", atom[i].coor[0]);
		 fprintf(fp,"%7.2f ", atom[i].coor[1]);
		 fprintf(fp,"%7.2f  ", atom[i].coor[2]);
		 fprintf(fp,"%4.2f  ", atom[i].r);
		 fprintf(fp,"%5.2f  ", atom[i].logp);
		 fprintf(fp,"%2s ", atom[i].hb);
		 fprintf(fp,"%7.2f ", atom[i].root[0]);
                 fprintf(fp,"%7.2f ", atom[i].root[1]);
                 fprintf(fp,"%7.2f", atom[i].root[2]);
		 fprintf(fp,"\n");

		 num++;
		}

	fprintf(fp,"<End>\n");

	fclose(fp);
			
	return;
}