autosphere.c

最经典的分子对结软件

/****************************************************************************************
**    AutoSphere:  A program to form a docking set of spheres from the 0 cluster 
**                 spheres genenerated by SPHGEN.
**
**    Created: 11-8-99
**    Author: Scott Pegg (Kuntz Group, Dept. of Pharm. Chem., UCSF)
**            Jim Arnold (Kuntz Group, Dept. of Pharm. Chem., UCSF)
**
**   Synopsis:
**
**     This program reads in an input file given as a command line argument.
**     This file specifies a SPHGEN output file, starting coordinates, an
**     enclosed space, and output filenames (see the README for more info about
**     the parameters) The spheres in "cluster 0" of the SPHGEN output file are
**     selected with the following procedure:
**
**     1. Remove spheres not within a box, large sphere, or close to a protein residue.
**     2. Remove spheres distant from the protein surface.
**     3. Remove spheres too close to other spheres, as such:
**           x = sphere closest to the center of the enclosed space
**           loop (until all spheres visited or removed):
**              remove all spheres closer than given filt_cut to x
**              x = next closest remaining sphere
**
**     The spheres are treated as points--ie. the distance between two spheres
**     is calculated as the distance between their centers. (I treat them this
**     way because they're treated as points in DOCK.)
**
**     Additions by Jim Arnold, Nov, 22, 1999.
**        - Remove spheres distant from a protein residue.
**        - Remove spheres distant from the protein surface.
**        - Read molecule files using Dock4.0.1 functions (archived Dock).
**        - Write out sphere sets in both .pdb and .sph formats.
**        - Added debug and guards to the program.
****************************************************************************************/



#include  "AutoSphere.h"
#include  <stdio.h>
#include  <malloc.h>
#include  <stdlib.h>
#include  <string.h>
#include  <math.h>
#include  "logical.h"
#include  "debug.h"

#include  "/idk1/people/arnold/fulcrum/Dock4.0.1/source/dock/define.h"
#include  "/idk1/people/arnold/fulcrum/Dock4.0.1/source/dock/global.h"
#include  "/idk1/people/arnold/fulcrum/Dock4.0.1/source/dock/utility.h"
#include  "/idk1/people/arnold/fulcrum/Dock4.0.1/source/dock/mol.h"
#include  "/idk1/people/arnold/fulcrum/Dock4.0.1/source/dock/label.h"



/**************
** mark defs.
**************/
#define    DEBUG_FILE      "as.dbf"                            /** file with debug flags. **/
#define    PARAM_DIR       "/idk1/people/arnold/fulcrum/Tsp"   /** directory for debug file. **/




/****************
** mark globals
****************/
GLOBAL global = {0};            /* See global.h for global variables */



/*******************
** mark prototypes
*******************/
void         read_args(FILE *file, TS_INFO *info);
int          TestParams(TS_INFO *info);
SPHERE*      read_spheres(TS_INFO *info, SPHERE *spheres);
SPHERE*      AddLigSpheres(MOLECULE *mol, SPHERE *spheres);
int          ReadMol(FILE_NAME molS, MOLECULE *mol_ref, MOLECULE *mol_init);
void         SetLabel(LABEL *label);
void         write_spheres(TS_INFO *info, SPHERE *spheres);
SPHERE*      find_enclosed_spheres(TS_INFO *info, SPHERE *spheres, MOLECULE *mol);
SPHERE*      FilterByResidue(TS_INFO *info, SPHERE *spheres, MOLECULE *mol);
SPHERE*      SurfDistThin(TS_INFO *info, SPHERE *spheres, MOLECULE *mol);
SPHERE*      greedy_thin(TS_INFO *info, SPHERE *spheres);
SPHERE*      add_sphere(SPHERE *sphere, SPHERE *head);
SPHERE*      remove_sphere(SPHERE *sphere, SPHERE *head);
int          SphereGuard(SPHERE *spheres, STRING200 errorS);
int          SphLen(SPHERE *sph);
int          IsHeavy(ATOM *at);
float        dist(SPHERE *a, SPHERE *b);
void         show_info(TS_INFO *info);
void         show_spheres(SPHERE *spheres);




/** mark main **/
/**************************************************************************
**  Author:   Scott Pegg, Jim Arnold.
**  Purpose:  Main routine for TS.
**  Notes:    
**************************************************************************/
main(int argc, char *argv[]) 
{
  const char   c[101];
  FILE         *param_file;
  TS_INFO      *info;
  SPHERE       *spheres = NULL;  /* head of the list of spheres */
  SPHERE       *temp, *temp2;
  MOLECULE     mol_ref = {0}, mol_init = {0}, lig_ref = {0}, lig_init = {0};

  if (argc == 1) {
    printf("Usage: TS <parameter file>\n");
    exit(0);
  }

  DebugInitI(PARAM_DIR, DEBUG_FILE);     /** read debug flags from file. **/
  if (MONITOR) {printf("\nRunning AutoSphere.\n");}

  param_file = fopen(argv[1], "r");
  if (param_file == NULL) {
    printf("Error: Could not open input file: %s\n", argv[1]);
    exit(0);
  }

  info = malloc(sizeof(struct ts_info));
  read_args(param_file, info);
  if (TestParams(info) equal true)
  {
     spheres = read_spheres(info, spheres);
     if (SphereGuard(spheres, "reading spheres, check the sph_file value in the .in file"))
     {
        if (MONITOR) {printf("  Read |%d| spheres from cluster 0 of file %s.\n", SphLen(spheres),info->sph_input);}

        if (info->lig_sph_file[0] != NULL)
	{
           allocate_molecule(&lig_ref);     lig_ref.info.name = NULL;     
           allocate_molecule(&lig_init);    lig_init.info.name = NULL;   
           ReadMol(info->lig_sph_file, &lig_ref, &lig_init);       
           if (MONITOR) {printf("  Read lig |%s| with %d atoms.\n",info->lig_sph_file, lig_ref.total.atoms);}
           spheres = AddLigSpheres(&lig_ref, spheres);
	}

        allocate_molecule(&mol_ref);     mol_ref.info.name = NULL;     
        allocate_molecule(&mol_init);    mol_init.info.name = NULL;   
        ReadMol(info->receptor_filename, &mol_ref, &mol_init);       
        if (MONITOR) {printf("  Read receptor |%s| with %d atoms.\n",info->receptor_filename,mol_ref.total.atoms);}

        spheres = find_enclosed_spheres(info, spheres, &mol_ref);
        if (SphereGuard(spheres, 
            "finding spheres in volume, check residue or res_cut values in .in file - or remove TERs from .pdb"))
        {
           if (MONITOR) {printf("  Removed spheres outside defined volume:  %d remain.\n", SphLen(spheres));}

           spheres = SurfDistThin(info, spheres, &mol_ref);
           if (SphereGuard(spheres, 
               "surf dist filter, check surface_cutoff value in the .in file, or mol/sph coords"))
  	   {
              if (MONITOR){printf("  Removed spheres distant from receptor atoms:  %d remain.\n", SphLen(spheres));}
 
              spheres = greedy_thin(info, spheres);
              if (SphereGuard(spheres, "greedy filter, check the filter_cutoff value in the .in file"))
   	      {
                 if (MONITOR) {printf("  Filtered spheres using greedy algorithm:  %d remain.\n", SphLen(spheres));}

                 if (MONITOR) 
                    {printf("  Writing %d spheres to .pdb and .sph: |%s|.\n",SphLen(spheres),info->sph_output);}
                 write_spheres(info, spheres);
  	      }
	   }
        }
     }
  }

  free(info);
  temp = spheres;
  while (temp != NULL) 
  {
    temp2 = temp->next;
    free(temp);
    temp = temp2;
  }

  if (MONITOR) {printf("  Exiting AutoSphere.\n\n");}
}

/** mark AddLigSpheres **/
/**************************************************************************
**  Author:   Jim Arnold
**  Purpose:  Adds spheres from the ligand to the 0 cluster spheres.
**  Notes:    
**************************************************************************/
SPHERE* AddLigSpheres(MOLECULE *mol, SPHERE *spheres)
{
   int      i, ctI;
   SPHERE   *new_sphere;

   if (DB_LIG) 
      {printf("**AddLigSpheres, mol: %s, num sph: %d.\n", mol->info.name,SphLen(spheres));}

   for (i = 0; i < mol->total.atoms; i++)
   {
      if (IsHeavy(&(mol->atom[i])))
      {
         new_sphere = malloc(sizeof(struct sphere));

         new_sphere->x = mol->coord[i][0]; 
         new_sphere->y = mol->coord[i][1]; 
         new_sphere->z =mol->coord[i][2];

         new_sphere->radius = 1.0;
         new_sphere->visited = 0;

         if (DB_LIG) {printf("  add lig sphere %d to the list.\n", ctI);}
         spheres = add_sphere(new_sphere, spheres);
         ctI++;
     }
   }

   if (DB_LIG) {printf("  ret: %d sph from AddLigSpheres, added: %d.\n", SphLen(spheres), ctI);}
   return (spheres);
}

/** mark SphereGuard **/
/**************************************************************************
**  Author:   Jim Arnold
**  Purpose:  Generates an error message if the number of spheres is 0.
**  Notes:    
**************************************************************************/
int SphereGuard(SPHERE *spheres, STRING200 errorS)
{
   int retI;

   if (DB_SPH_G) {printf("**SphereGuard, num sph: %d, error: |%s|.\n", SphLen(spheres), errorS);}

   retI = SphLen(spheres);

   if (retI equal 0)
   {
      printf("\n\nERROR:  %s.\n\n\n", errorS);
   }

   if (DB_SPH_G) {printf("  returning %d from SphereGuard.\n", retI);}
   return (retI);
}

/** mark find_enclosed_spheres **/
/**************************************************************************
**  Author:   Scott Pegg
**  Purpose:  A function to remove all spheres not within the specified 
**            enclosed area (a sphere or box).
**  Notes:    
**************************************************************************/
SPHERE* find_enclosed_spheres(TS_INFO *info, SPHERE *spheres, MOLECULE *mol) 
{
  SPHERE *temp, *temp2;
  SPHERE center;
  float x_lim1, x_lim2, y_lim1, y_lim2, z_lim1, z_lim2;

  if (DB_ENCLOSED) {printf("**find_enclosed_spheres, num sph: %d.\n", SphLen(spheres));}

  if (info->residue)
  {
     if (DB_ENCLOSED) {printf("  filter by residue num: %d.\n", info->res_num);}
     spheres = FilterByResidue(info, spheres, mol);
  }
  else if (info->sphere) 
  {
    if (DB_ENCLOSED) {printf("  filter by sphere.\n");}
    /* make a dummy sphere for the center of the enclosure */
    center.x = info->sph_x;
    center.y = info->sph_y;
    center.z = info->sph_z;

    /* for each sphere, check if it's inside the enclosure radius, and 
       remove it if it's not */ 
    temp = spheres;
    while (temp != NULL) 
    {
      temp2 = temp->next;
      if (dist(&center, temp) > info->radius) 
      {
	spheres = remove_sphere(temp, spheres);
      }
      temp = temp2;
    }
  }
  else if (info->box) 
  {
    if (DB_ENCLOSED) {printf("  filter with box.\n");}
    /* set the x,y,z limits on the coordinates (Yes, I know I don't really
       need to cast them explicity, but I like really clean compiles...) */
    x_lim1 = info->sph_x + (float)(info->size_x / 2.0);
    x_lim2 = info->sph_x - (float)(info->size_x / 2.0);
    y_lim1 = info->sph_y + (float)(info->size_y / 2.0);
    y_lim2 = info->sph_y - (float)(info->size_y / 2.0);
    z_lim1 = info->sph_z + (float)(info->size_z / 2.0);
    z_lim2 = info->sph_z - (float)(info->size_z / 2.0);

    /* remove any spheres with a coordinate outside the limits */
    temp = spheres;
    while (temp != NULL) 
    {
      temp2 = temp->next; 
      if (temp->x > x_lim1 || temp->x < x_lim2 || temp->y > y_lim1 ||
	  temp->y < y_lim2 || temp->z > z_lim1 || temp->z < z_lim2)
	spheres = remove_sphere(temp, spheres);
      temp = temp2;
    }
  }
  else
  {
     printf("ERROR:  No match for type of enclosed volume filtering.\n");
  }

  return spheres;
}

/** mark FilterByResidue **/
/**************************************************************************
**  Author:   Jim Arnold
**  Purpose:  Removes spheres not within a given distance of a specified
**            protein residue.
**  Notes:    
**************************************************************************/
SPHERE* FilterByResidue(TS_INFO *info, SPHERE *spheres, MOLECULE *mol)
{
   int             i, keep;
   float           d;
   SPHERE          *temp, *temp2;
   static SPHERE   dsp;

   if (DB_RES) 
   {
      printf("**FilterByRes, res: %d, sph: %d nat: %d.\n",info->res_num,SphLen(spheres),mol->total.atoms);
   }

   if (SphLen(spheres) <= 0)
   {
      if (DB_RES) {printf("  GUARD:  %d spheres passed to FilterByResidue.\n", SphLen(spheres));}
   }
   else
   {
      temp = spheres;
      while (temp != NULL)
      {
         keep = false;
         for (i = 0; i < mol->total.atoms; i++)
         {
            if (DB_RES) {printf("  atom: %d, res num: %d.\n", i, mol->atom[i].res_num);}
            if (mol->atom[i].res_num equal info->res_num)
	    {
               if (IsHeavy(&(mol->atom[i])))
      	       {
                  dsp.x = mol->coord[i][0]; dsp.y = mol->coord[i][1]; dsp.z = mol->coord[i][2];
                  d = dist(temp, &dsp);
                  if (DB_RES_L) {printf("    dist is %4.1f.\n", d);}
                  if (d < info->res_cut)
     	          {
                     if (DB_RES_L) {printf("      dist < %4.1f, keep sph.\n", d);}
                     keep = true;
      	          }
	       }
            }
	 }

         if (keep equal false)
         {
            if (DB_RES_L) {printf("    keep is false, remove sphere.\n");}
            temp2 = temp->next; /* temp get's freed if removed */
            spheres = remove_sphere(temp, spheres);
            temp = temp2;
         }