population.c

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

# include "process.h"

Population::Population()
{
	num_member=0;
	member=NULL;

	num_valid=0;
	pool=NULL;

	num_output=0;
	output=NULL;

	max_weight=0;
	min_weight=1000;
	ave_weight=0;

	max_binding_score=-100.00;
	min_binding_score= 100.00;
	ave_binding_score=   0.00;

	max_chemical_score=-1000.00;
	min_chemical_score= 1000.00;
	ave_chemical_score=    0.00;

	max_logp=-10.00;
	min_logp= 10.00;
	ave_logp=  0.00;
}

Population::~Population()
{
	delete [] member;
	delete [] output;
	delete [] sim_record;
}

void Population::Analyze_Population(char *filename)
{
	extern Parameter *parm;
	FILE *fp;
	int i;
	char buf[160],head[80];

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

	// first check the number of members

	num_member=0;

	for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else sscanf(buf,"%s",head);

                 if(strcmp(head,"<MOLECULE>")) continue;
                 else num_member++;
                }

	if(num_member==0) Lig_Format_Error(filename); 

	// then, read in the necessary information

	member=new Record[num_member];
	if(member==NULL) Memory_Allocation_Error();

	rewind(fp);

	for(i=0;i<num_member;i++)
		{
		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<MOLECULE>"));

                 sscanf(buf,"%*s%s",member[i].name);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<FORMULA>"));

                 sscanf(buf,"%*s%s",member[i].formula);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<WEIGHT>"));

                 sscanf(buf,"%*s%d",&member[i].weight);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<LOGP>"));

                 sscanf(buf,"%*s%f",&member[i].logp);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<BINDING_SCORE>"));

                 sscanf(buf,"%*s%f",&member[i].binding_score);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<CHEMICAL_SCORE>"));

                 sscanf(buf,"%*s%f",&member[i].chemical_score);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<INDEX>"));

                 sscanf(buf,"%*s%s",member[i].index);

		 member[i].valid=1; member[i].id=i+1;
		}

	fclose(fp);

	// give out the general information

	for(i=0;i<num_member;i++)
		{
		 if(member[i].weight>max_weight) max_weight=member[i].weight;
		 if(member[i].weight<min_weight) min_weight=member[i].weight;
		 ave_weight+=(member[i].weight);

		 if(member[i].binding_score>max_binding_score) 
			max_binding_score=member[i].binding_score;
		 if(member[i].binding_score<min_binding_score) 
			min_binding_score=member[i].binding_score;
		 ave_binding_score+=(member[i].binding_score);

		 if(member[i].chemical_score>max_chemical_score) 
			max_chemical_score=member[i].chemical_score;
		 if(member[i].chemical_score<min_chemical_score) 
			min_chemical_score=member[i].chemical_score;
		 ave_chemical_score+=(member[i].chemical_score);

		 if(member[i].logp>max_logp) max_logp=member[i].logp;
		 if(member[i].logp<min_logp) min_logp=member[i].logp;
		 ave_logp+=(member[i].logp);
		}

	ave_weight/=(num_member);
	ave_binding_score/=(num_member);
	ave_chemical_score/=(num_member);
	ave_logp/=(num_member);

	printf("\n****************************************************************\n");

	printf("There are %d molecules in '%s':\n", num_member, filename);
	printf("Molecular weight range: %d --- %d,  average = %d;\n",
		min_weight, max_weight, ave_weight);
	printf("LogP range: %6.2f --- %6.2f,  average = %6.2f;\n",
		min_logp, max_logp, ave_logp);
	printf("Binding score range: %6.2f --- %6.2f,  average = %6.2f;\n",
		min_binding_score, max_binding_score, ave_binding_score);
	// printf("Chemical score range: %7.2f --- %7.2f,  average = %7.2f;\n",
        //        min_chemical_score, max_chemical_score, ave_chemical_score); 

	printf("****************************************************************\n");

	// filter out unwanted members

	num_valid=0;

	for(i=0;i<num_member;i++)
		{
		 if((member[i].weight>parm->max_mol_weight)|| 
		    (member[i].weight<parm->min_mol_weight)||
		    (member[i].logp>parm->max_logp)||
		    (member[i].logp<parm->min_logp)||
		    (member[i].binding_score>parm->max_pkd)||
		    (member[i].binding_score<parm->min_pkd))
			{
			 member[i].valid=0;
			}
		 else 
			{
			 num_valid++;
			}
		}

	printf("\n%d molecules meet the user's criteria.\n", num_valid);

	return;
}

void Population::Extract_Members(char *filename)
{
	FILE *fp;
	int i,count,num;
	char buf[160],head[80];
	Ligand temp;

	if(num_valid<=0) return;

	pool=new Ligand[num_valid];
	if(pool==NULL) Memory_Allocation_Error();

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

	count=0; num=0;

	for(;;)
                {
                 if(fgets(buf,160,fp)==NULL) break;
                 else sscanf(buf,"%s",head);

                 if(strcmp(head,"<MOLECULE>")) continue;
                 else count++;

		 if(member[count-1].valid==0) continue; 

		 sscanf(buf,"%*s%s",temp.name);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<FORMULA>"));

                 sscanf(buf,"%*s%s",temp.formula);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<WEIGHT>"));

                 sscanf(buf,"%*s%d",&temp.weight);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<LOGP>"));

                 sscanf(buf,"%*s%f",&temp.logp);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<BINDING_SCORE>"));

                 sscanf(buf,"%*s%f",&temp.binding_score);

                 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<CHEMICAL_SCORE>"));

                 sscanf(buf,"%*s%f",&temp.chemical_score);

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<ATOM>"));

		 sscanf(buf,"%*s%d",&temp.num_atom);

		 delete [] temp.atom;

		 temp.atom=new Atom[temp.num_atom];
		 if(temp.atom==NULL) Memory_Allocation_Error();

        	 for(i=0;i<temp.num_atom;i++)
                	{
                 	 fgets(buf,160,fp);
                 	 sscanf(buf,"%d%s%f%f%f%s%s%*s%d%f",
                        	&temp.atom[i].id,
                         	 temp.atom[i].name,
                        	&temp.atom[i].coor[0],
                        	&temp.atom[i].coor[1],
                        	&temp.atom[i].coor[2],
                         	 temp.atom[i].type,
                         	 temp.atom[i].xtype,
                        	&temp.atom[i].valid,
                        	&temp.atom[i].score);
                	}

        	 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
        	 while(strcmp(head,"<BOND>"));

        	 sscanf(buf,"%*s%d",&temp.num_bond);

        	 delete [] temp.bond;

        	 temp.bond=new Bond[temp.num_bond];
        	 if(temp.bond==NULL) Memory_Allocation_Error();

        	 for(i=0;i<temp.num_bond;i++)
                	{
                 	 fgets(buf,160,fp);
                 	 sscanf(buf,"%d%d%d%s%*s%d",
                        	&temp.bond[i].id,
                        	&temp.bond[i].atom_1,
                        	&temp.bond[i].atom_2,
                         	 temp.bond[i].type,
                        	&temp.bond[i].valid);
                	}
		

		 do {fgets(buf,160,fp); sscanf(buf,"%s",head);}
                 while(strcmp(head,"<INDEX>"));

                 sscanf(buf,"%*s%s",temp.index);

		 temp.valid=1;

		 pool[num]=temp; num++;	
                }

	if(num!=num_valid)
		{
		 printf("Warning: something wrong with %s\n", filename);
		}

	fclose(fp);

	return;
}

void Population::Select_Results()
{
	extern Parameter *parm;
	int i,j,mark,tmp;
	float sim;
	Record temp;

	if(num_valid<=0) return;

	// first, build up the index for all the valid molecules

	delete [] member;

	member=new Record[num_valid];
	if(member==NULL) Memory_Allocation_Error();

	for(i=0;i<num_valid;i++)
		{
		 member[i].valid=1;
		 member[i].id=i+1;
		 member[i].binding_score=pool[i].binding_score;
		 member[i].chemical_score=pool[i].chemical_score;
		}