pocket.c

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

# include "pocket.h"

Pocket::Pocket()
{
	max_x=max_y=max_z=-999;
	min_x=min_y=min_z=999;

	grid=NULL; num_grid=0;
	feature=NULL; num_feature=0;
}

Pocket::~Pocket()
{
	delete [] grid;
	delete [] feature;
}

void Pocket::Define_Box()
{
	extern Protein *prot;
	extern Ligand *lig;
	int num,i,j,k,dx,dy,dz,tmp,neib;
	float d,d0,cutoff;
	struct Tmp_grid	{char type; float coor[3];};

	// first, make a rough box to contain all the pocket atoms

        for(i=0;i<prot->num_atom;i++)
        {
         if(prot->atom[i].valid!=2) continue;
         if(prot->atom[i].coor[0]>max_x) max_x=(int)prot->atom[i].coor[0];
         if(prot->atom[i].coor[0]<min_x) min_x=(int)prot->atom[i].coor[0];
         if(prot->atom[i].coor[1]>max_y) max_y=(int)prot->atom[i].coor[1];
         if(prot->atom[i].coor[1]<min_y) min_y=(int)prot->atom[i].coor[1];
         if(prot->atom[i].coor[2]>max_z) max_z=(int)prot->atom[i].coor[2];
         if(prot->atom[i].coor[2]<min_z) min_z=(int)prot->atom[i].coor[2];
        }

        max_x+=2; max_y+=2; max_z+=2;
        min_x-=2; min_y-=2; min_z-=2;   // add 2.0A margin to the box

/*
	printf("The rough box:\n");
        printf("Max_x = %d\tMin_x = %d\n", max_x, min_x);
        printf("Max_y = %d\tMin_y = %d\n", max_y, min_y);
        printf("Max_z = %d\tMin_z = %d\n", max_z, min_z);
*/

	// then, try to refine the box by cutting the edges

	cutoff=(DIST_CUTOFF)*(DIST_CUTOFF);

	// Now the grid spacing is 1.00A for a time of being

	dx=(max_x-min_x)+1;
	dy=(max_y-min_y)+1;
	dz=(max_z-min_z)+1;

	num_grid=dx*dy*dz;

	Tmp_grid *tmp_grid;
	tmp_grid=new Tmp_grid[num_grid];
	if(tmp_grid==NULL) Memory_Allocation_Error();

	for(num=0;num<num_grid;num++)
		{
		 tmp=num;
		 i=tmp/(dy*dz); tmp-=(i*dy*dz);
		 j=tmp/dz; tmp-=(j*dz);
		 k=tmp;

		 // tmp_grid[num]=grid[i][j][k];

		 tmp_grid[num].coor[0]=i*1.00+min_x;
		 tmp_grid[num].coor[1]=j*1.00+min_y;
		 tmp_grid[num].coor[2]=k*1.00+min_z;

		 // check if the grid is close to the ligand

		 neib=0;

                 for(tmp=0;tmp<lig->num_atom;tmp++)
                        {
                         if(lig->atom[tmp].valid<=0) continue;

                         d=Distance2(lig->atom[tmp].coor,tmp_grid[num].coor);

                         if(d<cutoff) {neib++; break;}
                         else continue;
                        }

		 if(neib>0) tmp_grid[num].type='V';		// vacant 
		 else {tmp_grid[num].type='E'; continue;}	// excluded

		 // then check if the grid bumps to any protein atom

		 neib=0;
		 for(tmp=0;tmp<prot->num_atom;tmp++)
			{
			 if(prot->atom[tmp].valid<=0) continue;

			 d=Distance(prot->atom[tmp].coor,tmp_grid[num].coor);
			 d0=prot->atom[tmp].r+0.50;  // using a hydrogen probe 

			 if(d<d0) {neib++; break;}
			 else continue;  
			}

		 if(neib!=0) {tmp_grid[num].type='E'; continue;}
		}

	// refine the box boundary

	max_x=max_y=max_z=-999;
	min_x=min_y=min_z=999;

	for(i=0;i<num_grid;i++)
	{
	 if(tmp_grid[i].type!='V') continue;
	 else
		{
		 if(tmp_grid[i].coor[0]>max_x) max_x=(int)tmp_grid[i].coor[0];
		 if(tmp_grid[i].coor[0]<min_x) min_x=(int)tmp_grid[i].coor[0];
		 if(tmp_grid[i].coor[1]>max_y) max_y=(int)tmp_grid[i].coor[1];
		 if(tmp_grid[i].coor[1]<min_y) min_y=(int)tmp_grid[i].coor[1];
		 if(tmp_grid[i].coor[2]>max_z) max_z=(int)tmp_grid[i].coor[2];
		 if(tmp_grid[i].coor[2]<min_z) min_z=(int)tmp_grid[i].coor[2];
		}
	}

	max_x++; max_y++; max_z++;	// add 1.0A margin to the box
	min_x--; min_y--; min_z--;

/*
	printf("The refined box:\n");
        printf("Max_x = %d\tMin_x = %d\n", max_x, min_x);
        printf("Max_y = %d\tMin_y = %d\n", max_y, min_y);
        printf("Max_z = %d\tMin_z = %d\n", max_z, min_z);
*/

	delete [] tmp_grid;

	return;
}

void Pocket::Make_Grid()
{
	extern Protein *prot;
	extern Ligand *lig;
	int num,i,j,k,dx,dy,dz,tmp,neib;
	float d,d0,cutoff;

	// Now the grid spacing is 0.50A 

	dx=(max_x-min_x)*2+1;
	dy=(max_y-min_y)*2+1;
	dz=(max_z-min_z)*2+1;

	num_grid=dx*dy*dz;

	grid=new Grid[num_grid];
	if(grid==NULL) Memory_Allocation_Error();

	// initialize all the grids as 'V'

	for(i=0;i<num_grid;i++)
		{
		 grid[i].valid=1;
		 grid[i].type='V';
		 grid[i].patom_neib=0;
		 grid[i].homo_neib=0;
		 grid[i].hetero_neib=0;
		 grid[i].donor_score=0.000;
		 grid[i].acceptor_score=0.000;
		 grid[i].hydrophobic_score=0.000;
		 grid[i].score=0.000;
		}

	// first, find out vacant grids and excluded grids

	cutoff=(DIST_CUTOFF)*(DIST_CUTOFF);

	for(num=0;num<num_grid;num++)
		{
		 tmp=num;
		 i=tmp/(dy*dz); tmp-=(i*dy*dz);
		 j=tmp/dz; tmp-=(j*dz);
		 k=tmp;

		 // grid[num]=grid[i][j][k];

		 grid[num].coor[0]=i*0.50+min_x;
		 grid[num].coor[1]=j*0.50+min_y;
		 grid[num].coor[2]=k*0.50+min_z;

		 // check if the grid is closed to the ligand

		 neib=0;

                 for(tmp=0;tmp<lig->num_atom;tmp++)
                        {
                         if(lig->atom[tmp].valid<=0) continue;

                         d=Distance2(lig->atom[tmp].coor,grid[num].coor);

                         if(d<cutoff) {neib++;break;}
                         else continue;
                        }

		 if(neib==0) {grid[num].type='E'; continue;}

		 // then check if the grid bumps to any pocket atom

		 neib=0;

		 for(tmp=0;tmp<prot->num_atom;tmp++)
			{
			 if(prot->atom[tmp].valid<=0) continue;

			 d=Distance(prot->atom[tmp].coor,grid[num].coor);
			 d0=prot->atom[tmp].r+0.50; 

			 if(d<d0) {neib++; break;}
			 else continue;
			}

		 if(neib>0) {grid[num].type='E'; continue;}

		 // then count the neighboring pocket atoms for each grid

		 neib=0;

		 for(tmp=0;tmp<prot->num_atom;tmp++)
			{
			 if(prot->atom[tmp].valid!=2) continue;

			 d=Distance2(prot->atom[tmp].coor,grid[num].coor);

			 if(d<cutoff) {neib++; continue;}
			 else continue;
			}

		 if(neib<=2) grid[num].type='S';	// empirical criteria
		 else grid[num].patom_neib=neib; 
		}

	// first find out and label the remaining grids outside the pocket

	// unfortunately, this algorithm does not work well

/*
	float tmp_max,tmp_min,tmp_ave,score;

	tmp=0; tmp_max=-999.0; tmp_min=999.0; tmp_ave=0.000; 

	for(i=0;i<num_grid;i++)
		{
		 if(grid[i].type=='E'||grid[i].type=='S') continue;
		 else
			{
			 score=grid[i].patom_neib;
			 if(score>tmp_max) tmp_max=score;
			 if(score<tmp_min) tmp_min=score;
			 tmp_ave+=score;
			 tmp++;
			}
		}

	tmp_ave/=tmp;

        printf("Grid neighbor:\n");
        printf("\tpoint number %d\n",tmp);
        printf("\tmaximum %6.2f\n",tmp_max);
        printf("\tminimum %6.2f\n",tmp_min);
        printf("\taverage %6.2f\n",tmp_ave);

	for(i=0;i<num_grid;i++)
		{
		 if(grid[i].type=='E'||grid[i].type=='S') continue;
		 else if(grid[i].patom_neib>tmp_ave) continue;
		 else grid[i].type='S';
		}
*/
			
	// second, calculate donor score for each grid 

	float r0=HB_D_R;			// the probe atom is N.4
	float angle;
	float v1[3],v2[3];
	struct
		{
		 int valid;
		 float d;
		 float overlap;
		 int type;	// 1 for ordinary h-bond 
				// 2 for water-involved h-bond
				// 3 for metal-acceptor bond
		 int loss;
		} candidate[20];
	int num_candidate;
		 
	for(num=0;num<num_grid;num++)
	{
	 if(grid[num].type=='E'||grid[num].type=='S') continue;

         for(i=0;i<prot->num_atom;i++)
                {
                 if(prot->atom[i].valid!=2) continue;
                 else
                        {
                         d0=r0+prot->atom[i].r;
                         d=Distance(grid[num].coor,prot->atom[i].coor);
                         if((d-d0)<-0.60) grid[num].donor_score+=(VB);
			 if(d<3.50&&prot->atom[i].q<0.000) 
				grid[num].donor_score+=SB_AWARD; 
                        }
                }

	 num_candidate=0;
	 for(i=0;i<20;i++) candidate[i].valid=0;

	 for(i=0;i<prot->num_atom;i++)
		{
		 if(prot->atom[i].valid!=2) continue;
		 else if(!strcmp(prot->atom[i].hb,"A")||
			 !strcmp(prot->atom[i].hb,"DA")) 
		{
		 d0=r0+prot->atom[i].r;
		 d=Distance(grid[num].coor,prot->atom[i].coor);
		 if(d>d0) continue;	// maximum h-bond length

		 for(j=0;j<3;j++)
			{
			 v1[j]=prot->atom[i].coor[j]-prot->atom[i].root[j];
			 v2[j]=prot->atom[i].coor[j]-grid[num].coor[j];
			}
		 angle=Angle_Of_Two_Vectors(v1,v2);
		 if(angle<HB_ANGLE_CUTOFF) continue;	// h-bond angle cutoff

		 candidate[num_candidate].valid=1;
		 candidate[num_candidate].d=d;
		 candidate[num_candidate].overlap=d-d0;
		 candidate[num_candidate].loss=0;

		 if(!strcmp(prot->atom[i].type,"O.w")) 
			candidate[num_candidate].type=2;
		 else candidate[num_candidate].type=1;

		 num_candidate++;
		}
		 else continue;
		}

	// The probe may form many HB pair with nearby pocket atoms
	// only the shortest HBs are taken into account

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

		 for(j=0;j<num_candidate;j++)
			{
			 if(i==j) continue;
			 else if(candidate[i].d<candidate[j].d) continue;
			 else candidate[i].loss++;
			}

		 if(candidate[i].loss>=4) 
			{
			 candidate[i].valid=0;
			 continue;
			}
		 else if(candidate[i].type==1)
			{
			 if(candidate[i].overlap<-0.60) 
                                grid[num].donor_score+=(SHB);
                         else if(candidate[i].overlap<-0.30)
                                grid[num].donor_score+=(MHB);
                         else if(candidate[i].overlap<0.00)
                                grid[num].donor_score+=(WHB);
                         else continue;
			}
		 else if(candidate[i].type==2)
			{
                         if(candidate[i].overlap<-0.60)
                                grid[num].donor_score+=(SWH);
                         else if(candidate[i].overlap<-0.30)
                                grid[num].donor_score+=(MWH);
                         else if(candidate[i].overlap<0.00)
                                grid[num].donor_score+=(WWH);
                         else continue;
			}
		 else continue;
		}
	
	}

	// third, calculate the acceptor score for each grid 

	r0=HB_A_R;	// probe atom O.co2

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

         for(i=0;i<prot->num_atom;i++)
                {
                 if(prot->atom[i].valid!=2) continue;
                 else
                        {
                         d0=r0+prot->atom[i].r;
                         d=Distance(grid[num].coor,prot->atom[i].coor);
                         if((d-d0)<-0.60) grid[num].acceptor_score+=(VB);
			 if(d<3.50&&prot->atom[i].q>0.000)
				grid[num].acceptor_score+=SB_AWARD;  
                        }
                }

	 num_candidate=0;
	 for(i=0;i<20;i++) candidate[i].valid=0;

	 for(i=0;i<prot->num_atom;i++)
		{
		 if(prot->atom[i].valid!=2) continue;
		 else if(!strcmp(prot->atom[i].hb,"D")||
			 !strcmp(prot->atom[i].hb,"DA")) 
		{