population.c

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

	// then re-arrange the index according to binding scores

	for(i=0;i<num_valid-1;i++)
	for(j=i+1;j<num_valid;j++)
		{
		 if(member[i].binding_score>=member[j].binding_score) continue;
		 else
			{
			 temp=member[i];
			 member[i]=member[j];
			 member[j]=temp;
			}
		}

	// now select top molecules into the output pool

	output=new Ligand[parm->num_output];
	if(output==NULL) Memory_Allocation_Error();

	sim_record=new float[parm->num_output*parm->num_output];
	if(sim_record==NULL) Memory_Allocation_Error();

	for(i=0;i<parm->num_output;i++)
	for(j=0;j<parm->num_output;j++)
		{
		 if(i==j) sim_record[i*parm->num_output+j]=1.000;
		 else sim_record[i*parm->num_output+j]=0.000;
		}

	num_output=0;

	for(i=0;i<num_valid;i++)
		{
         	 if(num_output>=parm->num_output) break;    // already full

		 tmp=member[i].id-1;

	 	 // similarity check

	 	 mark=FALSE;
	 	 for(j=0;j<num_output;j++)
         	{
	 	 sim=Cal_Molecular_Similarity(output[j],pool[tmp]);

		 sim_record[j*parm->num_output+num_output]=sim;
		 sim_record[num_output*parm->num_output+j]=sim;

	 	 if(sim<parm->similarity_cutoff) continue;
               	 else {mark=TRUE; break;}
               	}

	 	 if(mark==TRUE) continue;
		 else {output[num_output]=pool[tmp]; num_output++;}
        	}

	delete [] pool;

	return;
}

void Population::Cluster_Results()
{
	extern Parameter *parm;
	int i,j,tmp,mark;
        int matrix_size,no1,no2,old_id,new_id;
        float *matrix,*tmp_mat;
        float sim,ave_sim,top_sim;

	if(num_output<=1) return;

        // first, make the similarity matrix

        matrix_size=num_output*num_output;

        matrix=new float[matrix_size];
        if(matrix==NULL) Memory_Allocation_Error();

        for(i=0;i<matrix_size;i++) matrix[i]=1.000;

        ave_sim=0.000; tmp=0;

        for(i=0;i<num_output-1;i++)
        for(j=i+1;j<num_output;j++)
                {
		 sim=sim_record[i*parm->num_output+j];
                 matrix[i*num_output+j]=sim;
                 matrix[j*num_output+i]=sim;
                 ave_sim+=sim;
                 tmp++;
                }

        ave_sim/=tmp;

	// second, cluster the result molecules according to the matrix

        tmp_mat=new float[matrix_size];
        if(tmp_mat==NULL) Memory_Allocation_Error();

        for(i=0;i<matrix_size;i++) tmp_mat[i]=matrix[i];

        for(i=0;i<num_output;i++) output[i].valid=i+1;  // initialize clusters

	do
        {
         top_sim=0.000;

         for(i=0;i<num_output-1;i++)
         for(j=i+1;j<num_output;j++)
                {
                 if(matrix[i*num_output+j]<=top_sim) continue;
                 else {top_sim=matrix[i*num_output+j]; no1=i; no2=j;}
                }

         matrix[no1*num_output+no2]=0.000;
         matrix[no2*num_output+no1]=0.000;

         if(output[no1].valid==output[no2].valid) continue;
         else if(output[no1].valid>output[no2].valid)
                {
                 old_id=output[no1].valid;
                 new_id=output[no2].valid;
                }
         else
                {
                 old_id=output[no2].valid;
                 new_id=output[no1].valid;
                }

         mark=TRUE;

	 for(i=0;i<num_output;i++)
                {
                 if(output[i].valid!=new_id) continue;

                 for(j=0;j<num_output;j++)
                        {
                         if(i==j) continue;
                         else if(output[j].valid!=old_id) continue;
                         else if(tmp_mat[i*num_output+j]>=ave_sim) continue;
                         else {mark=FALSE; break;}
                        }

                 if(mark==FALSE) break;
                }

         if(mark==TRUE)         // merge two cluster
                {
                 for(i=0;i<num_output;i++)
                         {
                          if(output[i].valid!=old_id) continue;
                          else output[i].valid=new_id;
                         }
                }
         else continue;

        } while(top_sim>ave_sim);

        delete [] matrix; delete [] tmp_mat;

	// third, re-arrange the family IDs

        int num_family=0;
        int *family_list;

        family_list=new int[num_output];
        if(family_list==NULL) Memory_Allocation_Error();

        for(i=0;i<num_output;i++) family_list[i]=0;

        for(i=0;i<num_output;i++)
                {
                 mark=FALSE;

                 for(j=0;j<num_family;j++)
                        {
                         if(output[i].valid!=family_list[j]) continue;
                         else {mark=TRUE; break;}
                        }

                 if(mark==TRUE) {output[i].valid=j+1;}  // not a new family
                 else                                   // a new family
                        {
                         family_list[num_family]=output[i].valid;
                         num_family++;
                         output[i].valid=num_family;
                        }
                }

        delete [] family_list;

	// finally, arrange all the output molecules according to the family
        // those molecules within a family are ranked by binding score

        Ligand tmp_mol;

        for(i=0;i<num_output-1;i++)
        for(j=i+1;j<num_output;j++)
                {
                 if(output[i].valid==output[j].valid)
                {
                 if(output[i].binding_score>=output[j].binding_score) continue;
                 else
                        {
                         tmp_mol=output[i];
                         output[i]=output[j];
                         output[j]=tmp_mol;
                        }
                }
                 else if(output[i].valid<output[j].valid) continue;
                 else
                        {
                         tmp_mol=output[i];
                         output[i]=output[j];
                         output[j]=tmp_mol;
                        }
                }
	
	return;
}

void Population::Output_Results()
{
	extern Parameter *parm;
        FILE *fp;
        int i,num1,num2,num3,tmp;
        char name[80],filename[160],str1[2],str2[2],str3[2];

        if(num_output<=0) return;

        strcpy(filename,parm->output_dir); strcat(filename,"INDEX");
        if((fp=fopen(filename,"w"))==NULL) Openning_File_Error(filename);

        fprintf(fp,"#\n");
        fprintf(fp,"# This file is create by LigBuilder/Process\n");
        fprintf(fp,"# It lists the properties of the selected molecules.\n");
        fprintf(fp,"# Creation time: %s", Get_Time());
        fprintf(fp,"#\n");
        fprintf(fp,"# 1st column: ID of the molecule\n");
        fprintf(fp,"# 2nd column: name of the mol2 file\n");
        fprintf(fp,"# 3rd column: family of the molecule\n");
        fprintf(fp,"# 4th column: formula\n");
        fprintf(fp,"# 5th column: molecular weight\n");
        fprintf(fp,"# 6th column: calculated LogP\n");
        fprintf(fp,"# 7th column: binding score (pKd)\n");
        fprintf(fp,"# 8th column: chemical score\n");
        fprintf(fp,"#\n");

	// then, output the top ones to the assigned directory

        for(i=0;i<num_output;i++)
                {
                 tmp=i+1;

                 num1=tmp/100; tmp-=(num1*100);
                 num2=tmp/10;  tmp-=(num2*10);
                 num3=tmp;

                 sprintf(str1,"%d",num1);
                 sprintf(str2,"%d",num2);
                 sprintf(str3,"%d",num3);

                 strcpy(name,"result_");
                 strcat(name,str1);strcat(name,str2);strcat(name,str3);

                 strcpy(output[i].name,name);
                 strcat(name,".mol2");

                 fprintf(fp,"%-6d", i+1);
                 fprintf(fp,"%s\t", name);
                 fprintf(fp,"<%d>\t", output[i].valid);
                 fprintf(fp,"%-16s", output[i].formula);
                 fprintf(fp,"%d\t", output[i].weight);
                 fprintf(fp,"%6.2f\t", output[i].logp);
                 fprintf(fp,"%5.2f\t", output[i].binding_score);
                 fprintf(fp,"%d\n", (int)output[i].chemical_score);

                 strcpy(filename,parm->output_dir);
                 strcat(filename,name);
                 output[i].Write_Out_Mol2(filename);
                }

        fclose(fp);

	printf("%d molecules are dumped into '%s'.\n", 
		num_output, parm->output_dir);

	return;
}

void Population::Target_Similarity(char *filename)
{
        int i,j;
        Ligand target,temp;

        target.Read_From_Mol2(filename);

        for(i=0;i<num_output;i++)
        output[i].binding_score=Cal_Molecular_Similarity(target,output[i]);

        for(i=0;i<num_output-1;i++)
        for(j=i+1;j<num_output;j++)
                {
                 if(output[i].binding_score>=output[j].binding_score) continue;
                 else
                        {
                         temp=output[i];
                         output[i]=output[j];
                         output[j]=temp;
                        }
                }

        printf("\nSimilarity with the target molecule:\n");

        for(i=0;i<num_output;i++)
                {
                 printf("%s\t%5.3f\n",output[i].name,output[i].binding_score);
                }

        printf("\n");

        return;
}