pocket.c

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

			 (grid[i].acceptor_score>0.001)) type='A';
                 else if((grid[i].hydrophobic_score>grid[i].donor_score)&&
			 (grid[i].hydrophobic_score>grid[i].acceptor_score)&&
			 (grid[i].hydrophobic_score>0.001)) type='H';
                 else type=grid[i].type;	// 'V' & 'S'

		 if(type=='H')
                        {
                         fprintf(fp,"HETATM %5d %-3s ", i+1,"C");
                         fprintf(fp,"POK   %-5d   ", 3);
                         fprintf(fp,"%7.3f %7.3f %7.3f  1.00\n",
                                 grid[i].coor[0], grid[i].coor[1], grid[i].coor[2]);
                        }
		 else if(type=='D')
			{
			 fprintf(fp,"HETATM %5d %-3s ", i+1,"N");
		    	 fprintf(fp,"POK   %-5d   ", 1);
			 fprintf(fp,"%7.3f %7.3f %7.3f  1.00\n",
			 	 grid[i].coor[0], grid[i].coor[1], grid[i].coor[2]);
			}
		 else if(type=='A')
                        {
                         fprintf(fp,"HETATM %5d %-3s ", i+1,"O");
                         fprintf(fp,"POK   %-5d   ", 2);
                         fprintf(fp,"%7.3f %7.3f %7.3f  1.00\n",
                                 grid[i].coor[0], grid[i].coor[1], grid[i].coor[2]);
                        }
		 else continue;
                }

        fclose(fp);

	return;
}

void Pocket::Count_Grid_Neighbor(Grid &point,float r)
{
	int i;
	float d,cutoff;

	point.homo_neib=0; point.hetero_neib=0;

	if(point.type=='E'||point.type=='S'||point.type=='V') return; 

	cutoff=r*r;

	for(i=0;i<num_grid;i++)
		{
		 if(grid[i].type!=point.type) continue;
		 else
			{
			 d=Distance2(grid[i].coor,point.coor);
			 if(d>cutoff) continue;
			 else point.homo_neib++;
			}
		}

	return;
}

void Pocket::Output_Grid(char *filename) const
{
	FILE *fp;
	int i;

	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 grids within the binding pocket\n");
        fprintf(fp,"# Creation time: %s", Get_Time());
        fprintf(fp,"#\n");

	fprintf(fp,"# The first three lines record the grid boundary\n");
	fprintf(fp,"# The 4th line records the total number of the grids\n"); 
	fprintf(fp,"#\n");
	fprintf(fp,"# The 1st column: grid ID\n");
	fprintf(fp,"# The 2nd column: grid X coordinate\n");
	fprintf(fp,"# The 3rd column: grid Y coordinate\n");
	fprintf(fp,"# The 4th column: grid Z coordinate\n");
	fprintf(fp,"# The 5th column: grid type\n");
	fprintf(fp,"#                 E represents excluded grid;\n");
	fprintf(fp,"#                 S represents solvent grid;\n");
	fprintf(fp,"#                 V represents vacant grid;\n");
	fprintf(fp,"#                 D represents h-bond donor grid;\n");
	fprintf(fp,"#                 A represents h-bond acceptor grid;\n");
	fprintf(fp,"#                 H represents hydrophobic grid;\n");
	fprintf(fp,"#\n");

	fprintf(fp,"<Start>\n");
	fprintf(fp,"<X_boundary>\t%d\t%d\n",min_x,max_x);
	fprintf(fp,"<Y_boundary>\t%d\t%d\n",min_y,max_y);
	fprintf(fp,"<Z_boundary>\t%d\t%d\n",min_z,max_z);
	fprintf(fp,"<Grid_number>\t%d\n",num_grid);

	for(i=0;i<num_grid;i++)
		{
		 fprintf(fp,"%-6d  ",i+1);
		 fprintf(fp,"%7.2f  %7.2f  %7.2f  ", 
			 grid[i].coor[0], grid[i].coor[1], grid[i].coor[2]);
		 fprintf(fp,"%c\n", grid[i].type);
		}

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

	fclose(fp);

	return;
}

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

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

        fprintf(fp,"HEADER    User defined\n");
	fprintf(fp,"COMPND    Key interaction sites\n");
        fprintf(fp,"AUTHOR    Generated by LigBuilder/Pocket\n");

	fprintf(fp,"REMARK   1\n");
        fprintf(fp,"REMARK   1 Creation time %s", Get_Time());
	fprintf(fp,"REMARK   1\n");
        fprintf(fp,"REMARK   2 Nitrogen atoms represent h-bond donor sites\n");
	fprintf(fp,"REMARK   2 Oxygen atoms represent h-bond acceptor sites\n");
	fprintf(fp,"REMARK   2 Carbon atoms represent hydrophobic sites\n"); 

	num=1;

        for(i=0;i<num_grid;i++)
                {
                 if(grid[i].type=='E'||grid[i].type=='S') continue;
		 else if(grid[i].type=='H')
                        {
                         fprintf(fp,"HETATM %4d  %-3s ", num,"C");
                         fprintf(fp,"POK     3     ");
                         fprintf(fp,"%7.3f %7.3f %7.3f  1.00\n",
                                 grid[i].coor[0], grid[i].coor[1],
                                 grid[i].coor[2]);
			 num++;
                        }
		 else if(grid[i].type=='D')
			{
			 fprintf(fp,"HETATM %4d  %-3s ", num,"N");
			 fprintf(fp,"POK     1     ");
			 fprintf(fp,"%7.3f %7.3f %7.3f  1.00\n",
			 	 grid[i].coor[0], grid[i].coor[1],
				 grid[i].coor[2]);
			 num++;
			}
		 else if(grid[i].type=='A')
                        {
                         fprintf(fp,"HETATM %4d  %-3s ", num,"O");
			 fprintf(fp,"POK     2     ");
                         fprintf(fp,"%7.3f %7.3f %7.3f  1.00\n",
                                 grid[i].coor[0], grid[i].coor[1],
                                 grid[i].coor[2]);
			 num++;
                        }
		 else continue;
                }

	fprintf(fp,"END\n");

        fclose(fp);

	return;
}

void Pocket::Deduce_Pharmacophore(float min_d, int max_num)
{
	int i,j,k,mark,max_neib,id;
	float d;
	Feature temp;

	feature=new Feature[max_num];
	if(feature==NULL) Memory_Allocation_Error();

	for(i=0;i<max_num;i++)
		{
		 feature[i].valid=0;
		 feature[i].type=' ';
		 feature[i].coor[0]=feature[i].coor[1]=feature[i].coor[2]=0.000;
		 feature[i].neib=0;
		 feature[i].score=0.000;
		}

	// find the center of prospective features

	num_feature=0; 

	for(i=0;i<max_num;i++)
	{
	 max_neib=0; id=num_grid;

	 for(j=0;j<num_grid;j++)
		{
		 if(grid[j].type=='E'||grid[j].type=='S'||grid[j].type=='V') continue;

		 // check whether this grid is closed to any existing feature 

		 k=0; mark=FALSE;

		 while(k<num_feature)
			{
			 d=Distance(grid[j].coor,feature[k].coor);
			 if(d<min_d) {mark=TRUE;break;}
			 else {k++;continue;}
			}

		 if(mark==TRUE) continue;

		 if(grid[j].homo_neib<20) continue; // eliminate the minor ones
		 else if(grid[j].homo_neib<=max_neib) continue;
		 else {id=j; max_neib=grid[j].homo_neib;}
		}

	 if(id<num_grid)	// add a new feature center
		{	
		 feature[i].valid=1;
		 feature[i].type=grid[id].type;
		 feature[i].coor[0]=grid[id].coor[0];
		 feature[i].coor[1]=grid[id].coor[1];
		 feature[i].coor[2]=grid[id].coor[2];
		 feature[i].neib=grid[id].homo_neib;
		 Count_Feature_Score(feature[i],2.00);
		 num_feature++;
		}
	 else	
		{
		 feature[i].valid=0; 
		}

	}

	// now re-arrange the features according to their neighbor numbers

	for(i=0;i<num_feature-1;i++)
	for(j=i+1;j<num_feature;j++)
		{
		 if(feature[j].valid==0) continue;
		 else if(feature[i].valid==0)
			{
			 temp=feature[i]; feature[i]=feature[j]; feature[j]=temp;
			}
		 else if(feature[i].neib>=feature[j].neib) continue;
		 else
			{
			 temp=feature[i]; feature[i]=feature[j]; feature[j]=temp;
			}
		}

	// printf("Totally %d features are found.\n", num_feature);

	return;
}

void Pocket::Count_Feature_Score(Feature &point, float r)
{
	int i,num;
	float d,cutoff,score;

	cutoff=r*r;

	score=0.000; num=0;

	for(i=0;i<num_grid;i++)
		{
		 if(grid[i].type!=point.type) continue;
		 else
			{
			 d=Distance2(grid[i].coor,point.coor);
			 if(d>cutoff) continue;
			 else {score+=grid[i].score; num++;}
			}
		}

	score/=(float)num;

	point.score=score;

	return;
}

void Pocket::Output_Pharmacophore_PDB(char *filename) const
{
	FILE *fp;
        int i;

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

        fprintf(fp,"HEADER    User defined\n");
	fprintf(fp,"COMPND    Pharmacophore features\n");
        fprintf(fp,"AUTHOR    Generated by LigBuilder/Pocket\n");

	fprintf(fp,"REMARK   1\n");
        fprintf(fp,"REMARK   1 Creation time %s", Get_Time());
	fprintf(fp,"REMARK   1\n");
        fprintf(fp,"REMARK   2 Nitrogen atoms represent h-bond donor\n");
	fprintf(fp,"REMARK   2 Oxygen atoms represent h-bond acceptor\n");
	fprintf(fp,"REMARK   2 Carbon atoms represent hydrophobic center\n"); 

        for(i=0;i<num_feature;i++)
                {
		 if(feature[i].valid==0) continue;
		 else if(feature[i].type=='H')
                        {
                         fprintf(fp,"HETATM %4d  %-3s ",i+1,"C");
                         fprintf(fp,"POK     1     ");
                         fprintf(fp,"%7.3f %7.3f %7.3f %5.2f\n",
                                 feature[i].coor[0], feature[i].coor[1],
                                 feature[i].coor[2], feature[i].score);
                        }
		 else if(feature[i].type=='D')
			{
			 fprintf(fp,"HETATM %4d  %-3s ",i+1, "N");
			 fprintf(fp,"POK     1     ");
			 fprintf(fp,"%7.3f %7.3f %7.3f %5.2f\n",
			 	 feature[i].coor[0], feature[i].coor[1],
				 feature[i].coor[2], feature[i].score);
			}
		 else if(feature[i].type=='A')
                        {
                         fprintf(fp,"HETATM %4d  %-3s ",i+1,"O");
			 fprintf(fp,"POK     1     ");
                         fprintf(fp,"%7.3f %7.3f %7.3f %5.2f\n",
                                 feature[i].coor[0], feature[i].coor[1],
                                 feature[i].coor[2], feature[i].score);
                        }
		 else continue;
                }

	fprintf(fp,"END\n");

        fclose(fp);

	return;
}

void Pocket::Output_Pharmacophore_TXT(char *filename) const
{
        FILE *fp;
        int i,j,num;
	float d;

        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 lists the pharmacophore features detected\n");
        fprintf(fp,"# Creation time: %s", Get_Time());
        fprintf(fp,"#\n");

        fprintf(fp,"# The 1st column: feature ID number\n");
        fprintf(fp,"# The 2nd column: feature property\n");
	fprintf(fp,"# The 3rd column: feature symbol\n");
	fprintf(fp,"# The 4th column: X coordinate\n");
	fprintf(fp,"# The 5th column: Y coordinate\n");
	fprintf(fp,"# The 6th column: Z coordinate\n");
	fprintf(fp,"# The 7th column: average score\n"); 
	fprintf(fp,"# The 8th column: neighboring grids around this feature\n");
        fprintf(fp,"#\n");

	fprintf(fp,"<Start>\n");
	fprintf(fp,"<Feature_description>\n");

	num=1;

        for(i=0;i<num_feature;i++)
                {
                 if(feature[i].valid==0) continue;
                 else 
		{
		 fprintf(fp,"%2d  ",num);

		 if(feature[i].type=='D') fprintf(fp,"H-bond_donor_site      ");
		 else if(feature[i].type=='A') fprintf(fp,"H-bond_acceptor_site   ");
		 else fprintf(fp,"Hydrophobic_site       ");

		 fprintf(fp,"%c   ",feature[i].type);
		 fprintf(fp,"%7.2f  %7.2f  %7.2f  ", feature[i].coor[0],
			 feature[i].coor[1], feature[i].coor[2]); 
		 fprintf(fp,"%7.2f  ", feature[i].score);
		 fprintf(fp,"%d\n", feature[i].neib);
		
		 num++;
		}
		}

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

	fprintf(fp,"    ");
	for(i=1;i<num;i++) fprintf(fp,"%5d  ",i);
	fprintf(fp,"\n");

	num=1;

	for(i=0;i<num_feature;i++)
	{
	 if(feature[i].valid==0) continue;
	 fprintf(fp,"%2d  ",num);

	 for(j=0;j<num_feature;j++)
		{
		 if(feature[j].valid==0) continue;
		 else if(i==j) fprintf(fp,"-----  ");
		 else 
			{
			 d=Distance(feature[i].coor,feature[j].coor);
			 fprintf(fp,"%5.2f  ",d); 
			}
		}

	 fprintf(fp,"\n"); 
	
	 num++;
	} 

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

	fclose(fp);
			
	return;
}