Difference between revisions of "2017 DOCK tutorial 1 with PDB 4QMZ NEW"

From Rizzo_Lab
Jump to: navigation, search
(→‎V. Docking a Single Molecule for Pose Reproduction)
(→‎Rescoring Docked Molecules)
 
(32 intermediate revisions by the same user not shown)
Line 149: Line 149:
 
The image that appears should resemble this:
 
The image that appears should resemble this:
  
 +
[[Image:4qmz_spheres.png|thumb|center|1000px|4qmz receptor with all spheres]]
  
  
−
 
+
Then we need to select the spheres pertinent to our docking experiment. Usually these spheres will be the closest N spheres to the native ligand molecule.
−
Then we need to select the spheres pertinent to our docking experiment. Usually these speheres will be the closest N spheres to the native ligand molecule.
 
  
 
Run  
 
Run  
−
     sphere_selector 4qmz.rec.sph ../01.dockprep/4qmz.lig.mol2 8.0
+
     sphere_selector 4qmz.rec.sph ../01.dockprep/4qmz.lig.mol2 10.0
  
−
This command will select all of the spheres within 8.0 angstroms of the ligand and output them to selected_spheres.sph
+
This command will select all of the spheres within 10.0 angstroms of the ligand and output them to selected_spheres.sph
  
 
To visualize the spheres using Chimera as previously done:
 
To visualize the spheres using Chimera as previously done:
Line 165: Line 165:
 
Actions -> Atoms/Bonds -> sphere
 
Actions -> Atoms/Bonds -> sphere
 
The selected spheres with the receptor surface should look similar to that as seen below:
 
The selected spheres with the receptor surface should look similar to that as seen below:
 +
 +
[[Image:Selectedspheres.png|thumb|center|1000px|4qmz receptor surface with the selected spheres within 10.0 Angstroms]]
  
 
==IV. Generating Box and Grid==
 
==IV. Generating Box and Grid==
Line 180: Line 182:
  
 
This input designates to dock that:
 
This input designates to dock that:
−
we want to create a box
+
  we want to create a box,
−
the box length should be 8.0 Angstroms
+
  the box length should be 8.0 Angstroms,
−
use the selected spheres in the file desginated
+
  use the selected spheres in the file designated,
−
output the box to the file specified
+
  output the box to the file specified
  
  
Line 218: Line 220:
  
 
this script should output verbosely to the terminal and produce 2 output files, grid.bmp and grid.nrg, both binary files.
 
this script should output verbosely to the terminal and produce 2 output files, grid.bmp and grid.nrg, both binary files.
 +
 +
  grid -i grid.in > gridinfo.out
 +
 +
Check the standard output file gridinfo.out to ensure that individual charges on charged residues are integer values (-1 or +1).
 +
 +
  University of California at San Francisco, DOCK 4.0.1
 +
 
 +
  __________________Job_Information_________________
 +
  launch_time                    Wed Feb  8 17:05:50 2017
 +
  host_name                      unknown
 +
  memory_limit                  -1
 +
  working_directory              /home/campus.stonybrook.edu/ronassar/MySBFiles/projects/dock_tutorial/03.box-grid2_savegridoutput
 +
  user_name                      ronassar
 +
 
 +
  ________________General_Parameters________________
 +
  compute_grids                  yes
 +
  grid_spacing                  0.4
 +
  output_molecule                no
 +
 
 +
  ________________Scoring_Parameters________________
 +
  contact_score                  no
 +
  energy_score                  yes
 +
  energy_cutoff_distance        9999
 +
  atom_model                    a
 +
  attractive_exponent            6
 +
  repulsive_exponent            12
 +
  distance_dielectric            yes
 +
  dielectric_factor              4
 +
  bump_filter                    yes
 +
  bump_overlap                  0.75
 +
 
 +
  ____________________File_Input____________________
 +
  receptor_file                  ../01.dockprep/4qmz.rec.mol2
 +
  box_file                      ../03.box-grid/4qmz.box.pdb
 +
  vdw_definition_file            /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
 +
 
 +
  ____________________File_Output___________________
 +
  score_grid_prefix              grid
 +
 
 +
 
 +
 
 +
  Reading in coordinates of receptor.
 +
  CHARGED RESIDUE MET                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE ASP                    :  -1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE ARG                    :    1.000
 +
  CHARGED RESIDUE LYS                    :    1.000
 +
  CHARGED RESIDUE GLU                    :  -1.000
 +
  CHARGED RESIDUE SER                    :  -1.000
 +
 
 +
  Total charge on 4qmz_rec_leap.pdb      :  -5.000
 +
 
 +
  Reading in grid box information.
 +
  Box center of mass                      :  -1.912    0.054  14.111
 +
  Box dimensions                          :  39.856  41.687  38.838
 +
  Number of grid points per side [x y z]  :      101      106      99
 +
  Total number of grid points            :  1059894
 +
 
 +
  Generating scoring grids.
 +
  Percent of protein atoms processed      :        0
 +
  Percent of protein atoms processed      :      10
 +
  Percent of protein atoms processed      :      20
 +
  Percent of protein atoms processed      :      30
 +
  Percent of protein atoms processed      :      40
 +
  Percent of protein atoms processed      :      50
 +
  Percent of protein atoms processed      :      60
 +
  Percent of protein atoms processed      :      70
 +
  Percent of protein atoms processed      :      80
 +
  Percent of protein atoms processed      :      90
 +
  Percent of protein atoms processed      :      100
 +
  Writing general grid info to grid.bmp
 +
  Writing bump grid to grid.bmp
 +
  Writing energy grids to grid.nrg
 +
    Writing attractive VDW energy grid
 +
    Writing repulsive VDW energy grid
 +
    Writing electrostatic energy grid
 +
 
 +
  Finished calculation.
  
 
==V. Docking a Single Molecule for Pose Reproduction==
 
==V. Docking a Single Molecule for Pose Reproduction==
Line 223: Line 381:
 
Create or enter directory 4, 04.dock
 
Create or enter directory 4, 04.dock
  
−
Minimization
+
===Minimization===
  
 
create min.in
 
create min.in
Line 291: Line 449:
  
 
This structure can be visualized using chimera and loading the surface into chimera along with the minimized and unminimized ligand structures:
 
This structure can be visualized using chimera and loading the surface into chimera along with the minimized and unminimized ligand structures:
 +
 +
[[Image:image_min.png|thumb|center|1000px|4qmz with ligand before and after minimization]]
 +
 +
The RMSD of the minimized ligand to crystal structure ligand is 0.297 A.
 +
 +
 +
At this point we are going to calculate the van der waals and electrostatic footprints for the unminimized and minimized ligand structures in relation to the receptor active site.
 +
 +
create the input file: footprint.in
 +
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          no
 +
  ligand_atom_file                                            ./4qmz.lig.min_scored.mol2
 +
  limit_max_ligands                                            no
 +
  skip_molecule                                                no
 +
  read_mol_solvation                                          no
 +
  calculate_rmsd                                              no
 +
  use_database_filter                                          no
 +
  orient_ligand                                                no
 +
  bump_filter                                                  no
 +
  score_molecules                                              yes
 +
  contact_score_primary                                        no
 +
  contact_score_secondary                                      no
 +
  grid_score_primary                                          no
 +
  grid_score_secondary                                        no
 +
  multigrid_score_primary                                      no
 +
  multigrid_score_secondary                                    no
 +
  dock3.5_score_primary                                        no
 +
  dock3.5_score_secondary                                      no
 +
  continuous_score_primary                                    no
 +
  continuous_score_secondary                                  no
 +
  footprint_similarity_score_primary                          yes
 +
  footprint_similarity_score_secondary                        no
 +
  fps_use_footprint_reference_mol2                            yes
 +
  fps_footprint_reference_mol2_filename                        ../01.dockprep/4qmz.lig.mol2
 +
  fps_foot_compare_type                                        Euclidean
 +
  fps_normalize_foot                                          no
 +
  fps_foot_comp_all_residue                                    yes
 +
  fps_receptor_filename                                        ../01.dockprep/4qmz.rec.mol2
 +
  fps_vdw_att_exp                                              6
 +
  fps_vdw_rep_exp                                              12
 +
  fps_vdw_rep_rad_scale                                        1
 +
  fps_use_distance_dependent_dielectric                        yes
 +
  fps_dielectric                                              4.0
 +
  fps_vdw_fp_scale                                            1
 +
  fps_es_fp_scale                                              1
 +
  fps_hb_fp_scale                                              0
 +
  pharmacophore_score_secondary                                no
 +
  descriptor_score_secondary                                  no
 +
  gbsa_zou_score_secondary                                    no
 +
  gbsa_hawkins_score_secondary                                no
 +
  SASA_score_secondary                                        no
 +
  amber_score_secondary                                        no
 +
  minimize_ligand                                              no
 +
  atom_model                                                  all
 +
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /opt/AMS536/dock6/parameters/flex.defn
 +
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        fps.min.output
 +
  write_footprints                                            yes
 +
  write_hbonds                                                yes
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 +
 +
this footprint calculation can be run in the same way as the minimization:
 +
  dock6 -i footprint.in
 +
 +
This will output three files: fps.min.output_scored.mol2
 +
                              fps.min.output_scored_footprint_scored.txt
 +
                              fps.min.output_hbond_scored.txt
 +
 +
With an in house script, plot_footprint_single_magnitude.py the two sets of footprints will be plotted for comparison
 +
 +
this script will be run using:
 +
  python plot_footprint_single_magnitude.py fps.min.output_scored_footprint_scored.txt 50
 +
 +
The output will look something like this:
 +
 +
[[Image:figure_fps_minimized.png|thumb|center|1000px|Footprint comparison between unminimized and minimized ligand structure]]
 +
 +
Notice the vdW clashes are now reduced (vdW positive peaks disappeared).
 +
 +
 +
===Rigid Docking===
 +
 +
To help identify any errors in preparing the ligand and receptor files, a rigid docking, fixed growth docking and flexible docking runs will be performed. This is a check that the program can put the ligand into the crystal structure pose before proceeding.
 +
 +
Create an input file 'rigid.in' for the rigid docking run with the following parameters:
 +
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100.0
 +
  ligand_atom_file                                            ../01.dockprep/4qmz.lig.mol2
 +
  limit_max_ligands                                            no
 +
  skip_molecule                                                no
 +
  read_mol_solvation                                          no
 +
  calculate_rmsd                                              yes
 +
  use_rmsd_reference_mol                                      yes
 +
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
 +
  use_database_filter                                          no
 +
  orient_ligand                                                yes
 +
  automated_matching                                          yes
 +
  receptor_site_file                                          ../02.surface-spheres/selected_spheres.sph
 +
  max_orientations                                            1000
 +
  critical_points                                              no
 +
  chemical_matching                                            no
 +
  use_ligand_spheres                                          no
 +
  bump_filter                                                  no
 +
  score_molecules                                              yes
 +
  contact_score_primary                                        no
 +
  contact_score_secondary                                      no
 +
  grid_score_primary                                          yes
 +
  grid_score_secondary                                        no
 +
  grid_score_rep_rad_scale                                    1
 +
  grid_score_vdw_scale                                        1
 +
  grid_score_es_scale                                          1
 +
  grid_score_grid_prefix                                      ../03.box-grid/grid
 +
  multigrid_score_secondary                                    no
 +
  dock3.5_score_secondary                                      no
 +
  continuous_score_secondary                                  no
 +
  footprint_similarity_score_secondary                        no
 +
  pharmacophore_score_secondary                                no
 +
  descriptor_score_secondary                                  no
 +
  gbsa_zou_score_secondary                                    no
 +
  gbsa_hawkins_score_secondary                                no
 +
  SASA_score_secondary                                        no
 +
  amber_score_secondary                                        no
 +
  minimize_ligand                                              yes
 +
  simplex_max_iterations                                      1000
 +
  simplex_tors_premin_iterations                              0
 +
  simplex_max_cycles                                          1
 +
  simplex_score_converge                                      0.1
 +
  simplex_cycle_converge                                      1.0
 +
  simplex_trans_step                                          1.0
 +
  simplex_rot_step                                            0.1
 +
  simplex_tors_step                                            10.0
 +
  simplex_random_seed                                          0
 +
  simplex_restraint_min                                        no
 +
  atom_model                                                  all
 +
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /opt/AMS536/dock6/parameters/flex.defn
 +
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        4qmz.rigid_norestraint
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        50
 +
  write_conformations                                          no
 +
  cluster_conformations                                        yes
 +
  cluster_rmsd_threshold                                      2.0
 +
  rank_ligands                                                no
 +
 +
Outputs the file "4qmz.rigid_norestraint" the no restraint refers to the option of the parameter simplex_restraint_min in the input file. This removes the tether when doing the minimization. The best scoring structure had an RMSD of 1.??? to the native ligand pose.
 +
 +
[[Image:rigiddockingresults.png|thumb|center|800px|Rigid docking results]]
 +
 +
===Fixed Anchor Docking===
 +
 +
In fixed anchor docking, the ligand location is fixed in the active site and the rotatable bonds space is searched (i.e angles are flexible).
 +
 +
Create an input file 'fad.in' for the rigid docking run with the following parameters:
 +
Notice the search algorithm is now rigidset to flex (for flexible) and orient option is set to no (for fixing the anchor in place).
 +
 +
  conformer_search_type                                        flex
 +
  user_specified_anchor                                        no
 +
  limit_max_anchors                                            no
 +
  min_anchor_size                                              5
 +
  pruning_use_clustering                                      yes
 +
  pruning_max_orients                                          1000
 +
  pruning_clustering_cutoff                                    100
 +
  pruning_conformer_score_cutoff                              100.0
 +
  pruning_conformer_score_scaling_factor                      1.0
 +
  use_clash_overlap                                            no
 +
  write_growth_tree                                            no
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100.0
 +
  ligand_atom_file                                            ../01.dockprep/4qmz.lig.mol2
 +
  limit_max_ligands                                            no
 +
  skip_molecule                                                no
 +
  read_mol_solvation                                          no
 +
  calculate_rmsd                                              yes
 +
  use_rmsd_reference_mol                                      yes
 +
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
 +
  use_database_filter                                          no
 +
  orient_ligand                                                no
 +
  automated_matching                                          yes
 +
  receptor_site_file                                          ../02.surface-spheres/selected_spheres.sph
 +
  max_orientations                                            1000
 +
  critical_points                                              no
 +
  chemical_matching                                            no
 +
  use_ligand_spheres                                          no
 +
  bump_filter                                                  no
 +
  score_molecules                                              yes
 +
  contact_score_primary                                        no
 +
  contact_score_secondary                                      no
 +
  grid_score_primary                                          yes
 +
  grid_score_secondary                                        no
 +
  grid_score_rep_rad_scale                                    1
 +
  grid_score_vdw_scale                                        1
 +
  grid_score_es_scale                                          1
 +
  grid_score_grid_prefix                                      ../03.box-grid/grid
 +
  multigrid_score_secondary                                    no
 +
  dock3.5_score_secondary                                      no
 +
  continuous_score_secondary                                  no
 +
  footprint_similarity_score_secondary                        no
 +
  pharmacophore_score_secondary                                no
 +
  descriptor_score_secondary                                  no
 +
  gbsa_zou_score_secondary                                    no
 +
  gbsa_hawkins_score_secondary                                no
 +
  SASA_score_secondary                                        no
 +
  amber_score_secondary                                        no
 +
  minimize_ligand                                              yes
 +
  minimize_anchor                                              yes
 +
  minimize_flexible_growth                                    yes
 +
  use_advanced_simplex_parameters                              no
 +
  simplex_max_cycles                                          1
 +
  simplex_score_converge                                      0.1
 +
  simplex_cycle_converge                                      1.0
 +
  simplex_trans_step                                          1.0
 +
  simplex_rot_step                                            0.1
 +
  simplex_tors_step                                            10.0
 +
  simplex_anchor_max_iterations                                500
 +
  simplex_grow_max_iterations                                  500
 +
  simplex_grow_tors_premin_iterations                          0
 +
  simplex_random_seed                                          0
 +
  simplex_restraint_min                                        no
 +
  atom_model                                                  all
 +
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /opt/AMS536/dock6/parameters/flex.defn
 +
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        4qmz.fad
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        100
 +
  write_conformations                                          no
 +
  cluster_conformations                                        yes
 +
  cluster_rmsd_threshold                                      2.0
 +
  rank_ligands                                                no
 +
 +
[[Image:fixedanchorresults.png|thumb|center|800px|Fixed anchor docking results]]
 +
 +
===Flexible Docking===
 +
 +
In flexible docking, the ligand location and rotatable bonds are search.
 +
 +
Create an input file 'flex.in' for the flexible docking run with the following parameters:
 +
Notice the search algorithm is now set to flex (for flexible) and orient option is set to yes (for fixing the anchor in place).
 +
 +
  conformer_search_type                                        flex
 +
  user_specified_anchor                                        no
 +
  limit_max_anchors                                            no
 +
  min_anchor_size                                              5
 +
  pruning_use_clustering                                      yes
 +
  pruning_max_orients                                          1000
 +
  pruning_clustering_cutoff                                    100
 +
  pruning_conformer_score_cutoff                              100.0
 +
  pruning_conformer_score_scaling_factor                      1.0
 +
  use_clash_overlap                                            no
 +
  write_growth_tree                                            no
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100.0
 +
  ligand_atom_file                                            ../01.dockprep/4qmz.lig.mol2
 +
  limit_max_ligands                                            no
 +
  skip_molecule                                                no
 +
  read_mol_solvation                                          no
 +
  calculate_rmsd                                              yes
 +
  use_rmsd_reference_mol                                      yes
 +
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
 +
  use_database_filter                                          no
 +
  orient_ligand                                                yes
 +
  automated_matching                                          yes
 +
  receptor_site_file                                          ../02.surface-spheres/selected_spheres.sph
 +
  max_orientations                                            1000
 +
  critical_points                                              no
 +
  chemical_matching                                            no
 +
  use_ligand_spheres                                          no
 +
  bump_filter                                                  no
 +
  score_molecules                                              yes
 +
  contact_score_primary                                        no
 +
  contact_score_secondary                                      no
 +
  grid_score_primary                                          yes
 +
  grid_score_secondary                                        no
 +
  grid_score_rep_rad_scale                                    1
 +
  grid_score_vdw_scale                                        1
 +
  grid_score_es_scale                                          1
 +
  grid_score_grid_prefix                                      ../03.box-grid/grid
 +
  multigrid_score_secondary                                    no
 +
  dock3.5_score_secondary                                      no
 +
  continuous_score_secondary                                  no
 +
  footprint_similarity_score_secondary                        no
 +
  pharmacophore_score_secondary                                no
 +
  descriptor_score_secondary                                  no
 +
  gbsa_zou_score_secondary                                    no
 +
  gbsa_hawkins_score_secondary                                no
 +
  SASA_score_secondary                                        no
 +
  amber_score_secondary                                        no
 +
  minimize_ligand                                              yes
 +
  minimize_anchor                                              yes
 +
  minimize_flexible_growth                                    yes
 +
  use_advanced_simplex_parameters                              no
 +
  simplex_max_cycles                                          1
 +
  simplex_score_converge                                      0.1
 +
  simplex_cycle_converge                                      1.0
 +
  simplex_trans_step                                          1.0
 +
  simplex_rot_step                                            0.1
 +
  simplex_tors_step                                            10.0
 +
  simplex_anchor_max_iterations                                500
 +
  simplex_grow_max_iterations                                  500
 +
  simplex_grow_tors_premin_iterations                          0
 +
  simplex_random_seed                                          0
 +
  simplex_restraint_min                                        no
 +
  atom_model                                                  all
 +
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /opt/AMS536/dock6/parameters/flex.defn
 +
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        4qmz.flex
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        100
 +
  write_conformations                                          no
 +
  cluster_conformations                                        yes
 +
  cluster_rmsd_threshold                                      2.0
 +
  rank_ligands                                                no
 +
 +
[[Image:flexibledockingresults.png|thumb|center|800px|Flexible docking results]]
  
 
==VI. Virtual Screening==
 
==VI. Virtual Screening==
 +
 +
At this point we will begin docking potential inhibitors into the active site of 4qmz. The compound library we will use contains 25,000 molecules and was previously prepared for another docking experiment. The compound library contains molecules with rotatable bonds number between 7 and 15. The library used is titled small_ligand_library.mol2, this library can be found at /gpfs/projects/AMS536/virtual_screen_library on seawulf.
 +
 +
We began by copying all of our local directories outlined in the first section to a project directory on seawulf to allow use to expedite the docking process over the 25,000 compounds in our test library.
 +
 +
We will create a virtual-screen.in file that looks like the following:
 +
 +
  conformer_search_type                                        flex
 +
  user_specified_anchor                                        no
 +
  limit_max_anchors                                            no
 +
  min_anchor_size                                              5
 +
  pruning_use_clustering                                      yes
 +
  pruning_max_orients                                          1000
 +
  pruning_clustering_cutoff                                    100
 +
  pruning_conformer_score_cutoff                              100.0
 +
  pruning_conformer_score_scaling_factor                      1.0
 +
  use_clash_overlap                                            no
 +
  write_growth_tree                                            no
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100.0
 +
  ligand_atom_file                                            small_ligand_library.mol2
 +
  limit_max_ligands                                            no
 +
  skip_molecule                                                no
 +
  read_mol_solvation                                          no
 +
  calculate_rmsd                                              yes
 +
  use_rmsd_reference_mol                                      yes
 +
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
 +
  use_database_filter                                          no
 +
  orient_ligand                                                yes
 +
  automated_matching                                          yes
 +
  receptor_site_file                                          ../02.surface-spheres/selected_spheres.sph
 +
  max_orientations                                            1000
 +
  critical_points                                              no
 +
  chemical_matching                                            no
 +
  use_ligand_spheres                                          no
 +
  bump_filter                                                  no
 +
  score_molecules                                              yes
 +
  contact_score_primary                                        no
 +
  contact_score_secondary                                      no
 +
  grid_score_primary                                          yes
 +
  grid_score_secondary                                        no
 +
  grid_score_rep_rad_scale                                    1
 +
  grid_score_vdw_scale                                        1
 +
  grid_score_es_scale                                          1
 +
  grid_score_grid_prefix                                      ../03.box-grid/grid
 +
  multigrid_score_secondary                                    no
 +
  dock3.5_score_secondary                                      no
 +
  continuous_score_secondary                                  no
 +
  footprint_similarity_score_secondary                        no
 +
  pharmacophore_score_secondary                                no
 +
  descriptor_score_secondary                                  no
 +
  gbsa_zou_score_secondary                                    no
 +
  gbsa_hawkins_score_secondary                                no
 +
  SASA_score_secondary                                        no
 +
  amber_score_secondary                                        no
 +
  minimize_ligand                                              yes
 +
  minimize_anchor                                              yes
 +
  minimize_flexible_growth                                    yes
 +
  use_advanced_simplex_parameters                              no
 +
  simplex_max_cycles                                          1
 +
  simplex_score_converge                                      0.1
 +
  simplex_cycle_converge                                      1.0
 +
  simplex_trans_step                                          1.0
 +
  simplex_rot_step                                            0.1
 +
  simplex_tors_step                                            10.0
 +
  simplex_anchor_max_iterations                                500
 +
  simplex_grow_max_iterations                                  500
 +
  simplex_grow_tors_premin_iterations                          0
 +
  simplex_random_seed                                          0
 +
  simplex_restraint_min                                        no
 +
  atom_model                                                  all
 +
  vdw_defn_file                                            /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn
 +
  flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        4qmz.virtualscreen
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  write_conformations                                          no
 +
  cluster_conformations                                        yes
 +
  cluster_rmsd_threshold                                      2.0
 +
  rank_ligands                                                no
 +
 +
This script is submitted to the seawulf cluster using a script called screen.sh which contains the following:
 +
 +
  #! /bin/tcsh
 +
  #PBS -l nodes=4:ppn=28
 +
  #PBS -l walltime=24:00:00
 +
  #PBS -o jobstatus.out
 +
  #PBS -e joberror.err
 +
  #PBS -N large_vs
 +
  #PBS -V
 +
  #PBS -q long
 +
  cd $PBS_O_WORKDIR
 +
  mpirun -np 112 dock6.mpi -i virtual-screen.in -o virtualscreen1.out
 +
 +
 +
This screen.sh input specifies that we will use 112 processors to dock our virtual library. We have a wall clock for the job of 24 hours and we will output jobstatus.out, joberror.err and virtualscreen1.out. In addition to these output files there is an output for the work performed by each processor individually.
 +
 +
===Cartesian Minimization of the docked molecules===
 +
 +
Now we will do cartesian minimization of the docked molecules. Cartesian here is referring to the receptor being explicitly present during the minimization.
 +
 +
Write the minimization input script minimize-vs.in:
 +
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100.0
 +
  ligand_atom_file                                            4qmz.virtualscreen_scored.mol2
 +
  limit_max_ligands                                            no
 +
  skip_molecule                                                no
 +
  read_mol_solvation                                          no
 +
  calculate_rmsd                                              no
 +
  use_database_filter                                          no
 +
  orient_ligand                                                no
 +
  bump_filter                                                  no
 +
  score_molecules                                              yes
 +
  contact_score_primary                                        no
 +
  contact_score_secondary                                      no
 +
  grid_score_primary                                          no
 +
  grid_score_secondary                                        no
 +
  multigrid_score_primary                                      no
 +
  multigrid_score_secondary                                    no
 +
  dock3.5_score_primary                                        no
 +
  dock3.5_score_secondary                                      no
 +
  continuous_score_primary                                    yes
 +
  continuous_score_secondary                                  no
 +
  cont_score_rec_filename                                      ../01.dockprep/4qmz.rec.mol2
 +
  cont_score_att_exp                                          6
 +
  cont_score_rep_exp                                          12
 +
  cont_score_rep_rad_scale                                    1
 +
  cont_score_use_dist_dep_dielectric                          yes
 +
  cont_score_dielectric                                        4.0
 +
  cont_score_vdw_scale                                        1
 +
  cont_score_es_scale                                          1
 +
  footprint_similarity_score_secondary                        no
 +
  pharmacophore_score_secondary                                no
 +
  descriptor_score_secondary                                  no
 +
  gbsa_zou_score_secondary                                    no
 +
  gbsa_hawkins_score_secondary                                no
 +
  SASA_score_secondary                                        no
 +
  amber_score_secondary                                        no
 +
  minimize_ligand                                              yes
 +
  simplex_max_iterations                                      1000
 +
  simplex_tors_premin_iterations                              0
 +
  simplex_max_cycles                                          1
 +
  simplex_score_converge                                      0.1
 +
  simplex_cycle_converge                                      1.0
 +
  simplex_trans_step                                          1.0
 +
  simplex_rot_step                                            0.1
 +
  simplex_tors_step                                            10.0
 +
  simplex_random_seed                                          0
 +
  simplex_restraint_min                                        no
 +
  atom_model                                                  all
 +
  vdw_defn_file                                                /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn
 +
  flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        4qmz.virtualscreen1.minimized
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 +
Submit the minimization job using with the script "job-min.sh" following lines:
 +
 +
  #!/bin/tcsh
 +
  #PBS -l nodes=16:ppn=28
 +
  #PBS -l walltime=12:00:00
 +
  #PBS -o jobstatusmin.out
 +
  #PBS -e jobstatusmin.err
 +
  #PBS -N chunk0min
 +
  #PBS -V
 +
  #PBS -q medium
 +
  module load intel/mpi/64/2017/0.098
 +
  cd $PBS_O_WORKDIR
 +
  mpirun -np 448 /gpfs/projects/AMS536/zzz.programs/dock6/bin/dock6.mpi -i minimize-vs.in -o minimize-vs.out
 +
 +
 +
At this point it would be good to rank our docked ligands by score and extract the 100 best scoring (100 lowest scoring, most negative) ligands to investigate further.
 +
 +
 +
===Rescoring Docked Molecules===
 +
 +
Dock has several built in scoring functions. In an effort to utilize these alternative scoring functions it is worthwhile to rescore the results from the initial virtual screen that were collected. The virtual screen will be rescored using footprint similarity, pharmacophore score, tanimoto score, the hungarian and the volume overlap score. Each of these scoring functions focuses in on a different chemical aspect of the ligand in question. The rescore.in script for this run can be found here:
 +
 +
conformer_search_type                                        rigid
 +
use_internal_energy                                          yes
 +
internal_energy_rep_exp                                      12
 +
internal_energy_cutoff                                      100.0
 +
ligand_atom_file                            ../05.virtual- screen/4qmz.virtualscreen1.minimized_scored.mol2
 +
limit_max_ligands                                            no
 +
skip_molecule                                                no
 +
read_mol_solvation                                          no
 +
calculate_rmsd                                              no
 +
use_database_filter                                          no
 +
orient_ligand                                                no
 +
bump_filter                                                  no
 +
score_molecules                                              yes
 +
contact_score_primary                                        no
 +
contact_score_secondary                                      no
 +
grid_score_primary                                          no
 +
grid_score_secondary                                        no
 +
multigrid_score_primary                                      no
 +
multigrid_score_secondary                                    no
 +
dock3.5_score_primary                                        no
 +
dock3.5_score_secondary                                      no
 +
continuous_score_primary                                    no
 +
continuous_score_secondary                                  no
 +
footprint_similarity_score_primary                          no
 +
footprint_similarity_score_secondary                        no
 +
pharmacophore_score_primary                                  no
 +
pharmacophore_score_secondary                                no
 +
descriptor_score_primary                                    yes
 +
descriptor_score_secondary                                  no
 +
descriptor_use_grid_score                                    no
 +
descriptor_use_multigrid_score                              no
 +
descriptor_use_continuous_energy                            no
 +
descriptor_use_footprint_similarity                          yes
 +
descriptor_use_pharmacophore_score                          yes
 +
descriptor_use_tanimoto                                      yes
 +
descriptor_use_hungarian                                    yes
 +
descriptor_use_volume_overlap                                yes
 +
descriptor_fps_use_footprint_reference_mol2                  yes
 +
descriptor_fps_footprint_reference_mol2_filename            ../04.dock/4qmz.lig.min_scored.mol2
 +
descriptor_fps_foot_compare_type                            Euclidean
 +
descriptor_fps_normalize_foot                                no
 +
descriptor_fps_foot_comp_all_residue                        yes
 +
descriptor_fps_receptor_filename                            ../01.dockprep/4qmz.rec.mol2
 +
descriptor_fps_vdw_att_exp                                  6
 +
descriptor_fps_vdw_rep_exp                                  12
 +
descriptor_fps_vdw_rep_rad_scale                            1
 +
descriptor_fps_use_distance_dependent_dielectric            yes
 +
descriptor_fps_dielectric                                    4.0
 +
descriptor_fps_vdw_fp_scale                                  1
 +
descriptor_fps_es_fp_scale                                  1
 +
descriptor_fps_hb_fp_scale                                  0
 +
descriptor_fms_score_use_ref_mol2                            yes
 +
descriptor_fms_score_ref_mol2_filename                      ../04.dock/4qmz.lig.min_scored.mol2
 +
descriptor_fms_score_write_reference_pharmacophore_mol2      no
 +
descriptor_fms_score_write_reference_pharmacophore_txt      no
 +
descriptor_fms_score_write_candidate_pharmacophore          no
 +
descriptor_fms_score_write_matched_pharmacophore            no
 +
descriptor_fms_score_compare_type                            overlap
 +
descriptor_fms_score_full_match                              yes
 +
descriptor_fms_score_match_rate_weight                      5.0
 +
descriptor_fms_score_match_dist_cutoff                      1.0
 +
descriptor_fms_score_match_proj_cutoff                      0.7071
 +
descriptor_fms_score_max_score                              20
 +
descriptor_fingerprint_ref_filename                          ../04.dock/4qmz.lig.min_scored.mol2
 +
descriptor_hungarian_ref_filename                            ../04.dock/4qmz.lig.min_scored.mol2
 +
descriptor_hungarian_matching_coeff                          -5
 +
descriptor_hungarian_rmsd_coeff                              1
 +
descriptor_volume_reference_mol2_filename                    ../04.dock/4qmz.lig.min_scored.mol2
 +
descriptor_volume_overlap_compute_method                    analytical
 +
descriptor_weight_fps_score                                  1
 +
descriptor_weight_pharmacophore_score                        1
 +
descriptor_weight_fingerprint_tanimoto                      -1
 +
descriptor_weight_hungarian                                  1
 +
descriptor_weight_volume_overlap_score                      -1
 +
gbsa_zou_score_secondary                                    no
 +
gbsa_hawkins_score_secondary                                no
 +
SASA_score_secondary                                        no
 +
amber_score_secondary                                        no
 +
minimize_ligand                                              no
 +
atom_model                                                  all
 +
vdw_defn_file                    /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn
 +
flex_defn_file                    /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn
 +
flex_drive_file                  /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl
 +
chem_defn_file                    /gpfs/projects/AMS536/zzz.programs/dock6/parameters/chem.defn
 +
pharmacophore_defn_file          /gpfs/projects/AMS536/zzz.programs/dock6/parameters/ph4.defn
 +
ligand_outfile_prefix                                        descriptor.output
 +
write_footprints                                            yes
 +
write_hbonds                                                yes
 +
write_orientations                                          no
 +
num_scored_conformers                                        1
 +
rank_ligands                                                no
 +
 +
 +
From here the results should be sorted according the score based on the scoring function in question. A script for this can be written fairly quickly in whatever language you are most comfortable with. If you are not very experienced with computer programming, a python script developed in house can be found here that will search through the output .mol2 file and select the top N scoring molecules. These top scoring molecule s are taken from the virtual screen "top'N'_molid.txt" and finds and their coordinates saved to an output file "top'N'_molid.mol2".
 +
 +
rank_scores_get_coords.py:
 +
 +
  #!/usr/bin/python
 +
  import sys
 +
  from operator import itemgetter, attrgetter, methodcaller
 +
  # input mol2 file as command line argument
 +
  f0=open(sys.argv[1],'r') ; data0 = f0.readlines() ; f0.close()                  #open mol2 data file with  scored molecules
 +
  # define variables
 +
  nummol = 50                                                                      # number of molecules to  keep after ranking by score
 +
  #scoretype="Grid_Score:"                                                        # scoring function to rank  by
 +
  scoretype="Footprint_Similarity_Score:"
 +
  #--------------------------------------------------------------------------------------------------------------------------
 +
  list_molid=[]                                                              #create a list to store molecule ids
 +
  list_molscores=[]                                                          #create list to store scores for each respective molecule
 +
  for line in data0:                                                          #reads through the input file
 +
        if line != "\n":                                                    #if the line is not empty  continue
 +
                if line.split()[0]=="##########":                          #if the line starts with ############ then continue
 +
                                if line.split()[1]=="Name:":                #if the line has the next entry as Name: continue
 +
                                        molid = line.split()[2]            #then add the molID to list_molid
 +
                                        list_molid.append(molid)
 +
                                elif line.split()[1]==scoretype:
 +
                                        molscore = line.split()[2]
 +
                                        list_molscores.append(molscore)
 +
  #print list_molid
 +
  #print list_molscores
 +
  both = zip(list_molid,list_molscores)                                      #concatenates list_molid and list_molscores into one array
 +
  sorted_both=sorted(both, key=itemgetter(1), reverse=True)                  #sorts the "both" list
 +
  sorted_top=[]                                                            #names a new list sorted_top5
 +
  #for element in sorted_both:                                                #print out first n entries in sorted_both
 +
  #      print element[0], element[1]
 +
  o=open("top"+str(nummol)+"_molid.txt","w")                                #opens output file topx_molid.txt
 +
  for element in sorted_both[:nummol]:                                           
 +
        sorted_top.append(element[0])                                         
 +
        o.write(element[0]+"\n")                                          #stores the top nummol entries in sorted_top
 +
  o.close()
 +
  #print sorted_top                                                         
 +
  ##
 +
  ##--------------------------------------------------------------------------------------------------------------------------
 +
  ##
 +
  ## Get coordinates of molecules of interest from scored mol2
 +
  f1=open("top"+str(nummol)+"_molid.txt","r") ; data1=f1.readlines(); f1.close() # takes in sorted list of top molecules
 +
  f2=open(sys.argv[1],"r") ; data2=f2.readlines(); f2.close()                    # takes in mol2 results
 +
  identifier=0    # this will switch between 0 and 1. When it is 0 no line printing, when it is 1 we print the line.
 +
  molid_top=[]  # list of molecule ids that are of interest
 +
  # go through the first file (has top n molecules) and append those ids to molecule ids list
 +
  for line in data1: 
 +
        molid=line.split()[0]
 +
        molid_top.append(molid)
 +
  output=open("top"+str(nummol)+"_molid.mol2","w")
 +
  # go through the second file (the rescored mol2 file), the next commands print out the coordinates for the molecules of interest
 +
  for line in data2:
 +
  #      if line != "\n":  # skipping empty lines
 +
                if identifier==1:
 +
                        try:
 +
                                if line.split()[0]=="##########" and line.split()[1]=="Name:" and line.split()[2] not in molid_top:  # changing the identifier when hitting a new molecule that is not of interest
 +
                                        identifier=0
 +
                                        output.write('\n')     
 +
                                else:
 +
                                        output.write(line)  # else keep the identifier at 1 and keep printing
 +
                        except IndexError:
 +
                                continue 
 +
                try:   
 +
                        if identifier==0 and line.split()[0]=="##########" and line.split()[1]=="Name:" and line.split()[2] in molid_top:  # change identifier to 1 when we hit molecule of interest
 +
                                identifier=1
 +
                                output.write(line)
 +
                except IndexError:
 +
                                continue
 +
  output.close()
  
 
==VIII. Frequently Encountered Problems==
 
==VIII. Frequently Encountered Problems==
 +
 +
Spending the time to ensure your starting structure is reasonable is of the utmost importance and will save you much headache moving    forward in your study. Reading through the PDB file as downloaded from the protein data bank will give the researcher a lot of  information regarding issues or grey area encountered by the researchers responsible for publishing the structure initially. Threading a structure through some quality control software such as tleap, chimera or another comparable software package will ensure charges are reasonable, hydrogens are appropriately placed, and there are no missing residues in the peptide chain.
 +
 +
A basic understanding of your enzymes active site will also be beneficial while approaching an enzymatic system from a docking perspective. Some questions to ask oneself are: does my protein have an ion (Mg, Ca, etc.) in the active site? Does my ligand make important interactions with conserved water residues? If so these waters could be kept in with the receptor to aid in discovery. The alternative would be to design an inhibitor that displaces the water molecules and supplements favorable interactions mediated by the water between the ligand and protein.

Latest revision as of 17:10, 8 March 2017

For additional Rizzo Lab tutorials see DOCK Tutorials. Use this link Wiki Formatting as a reference for editing the wiki. This tutorial was developed collaboratively by a subsection of the AMS 536 class of 2017, using DOCK v6.8.

I. Introduction

DOCK

DOCK is a molecular docking program used in drug discovery. It was developed by Irwin D. Kuntz, Jr. and colleagues at UCSF (see UCSF DOCK). This program, given a protein binding site and a small molecule, tries to predict the correct binding mode of the small molecule in the binding site, and the associated binding energy. Small molecules with highly favorable binding energies could be new drug leads. This makes DOCK a valuable drug discovery tool. DOCK is typically used to screen massive libraries of millions of compounds against a protein to isolate potential drug leads. These leads are then further studied, and could eventually result in a new, marketable drug. DOCK works well as a screening procedure for generating leads, but is not currently as useful for optimization of those leads.

DOCK 6 uses an incremental construction algorithm called anchor and grow. It is described by a three-step process:

  1. Rigid portion of ligand (anchor) is docked by geometric methods.
  2. Non-rigid segments added in layers; energy minimized.
  3. The resulting configurations are 'pruned' and energy re-minimized, yielding the docked configurations.

4QMZ

In this tutorial we will use PDB code 4QMZ, the deposited crystal structure of MST3 in complex with SUNITINIB.

Organizing Directories

While performing docking, it is convenient to adopt a standard directory structure / naming scheme, so that files are easy to find / identify.For this tutorial, we will use something similar to the following: 

~username/AMS536-Spring2016/dock-tutorial/00.files/
                                         /01.dockprep/
                                         /02.surface-spheres/
                                         /03.box-grid/
                                         /04.dock/
                                         /05.large-virtual-screen/
                                         /06.virtual-screen/
                                         /07.footprint/
                                         /08.print_fps

In addition, most of the important files that are derived from the original crystal structure will be given a prefix that is the same as the PDB code, '4QMZ'. The following sections in this tutorial will adhere to this directory structure/naming scheme.

II. Preparing the Receptor and Ligand

Download the PDB file (4QMZ)

4QMZ was moved into 00.files

4qmz.pdb was copied to raw_4qmz.pdb

raw_4qmz.pdb was opened with VI terminal editor 

    The header information, connect records, ions (atoms 2333 and 2334) and waters were deleted

    Res 178 = TPO, or phosphonothreonine
    Res 178 (TPO) was renamed to THR (Threonine) and HETATM renamed to ATOM, in addition the acanonical atoms were removed from the pdb leaving a deprotonated threonine (Atoms 1311-1314 in 4qmz.pdb)

    Res B49 was renamed to LIG and made Chain B

raw_4qmz.pdb was copied twice to 4qmz_rec.pdb and 4qmz_lig.pdb

4qmz_rec.pdb was opened with VI terminal editor 

    LIG atoms, or chain B, was deleted and the file saved

4qmz_lig.pdb was opened with VI terminal editor 

    Protein atoms, or chain A, was deleted and the file saved

4qmz_rec.pdb was loaded into tleap as a quality control measure 

    tleap
    source leaprc.protein.ff14SB
    lin = loadpdb /path/to/4qmz_rec.pdb
         2340 Hydrogens added, 1 heavy atom added (CSER RES 299, Chain A, OXT 12) 
    check lin
    saveamberparm lin /path/to/4qmz_rec_leap.parm7 /path/to/4qmz_rec_leap.crd

Running the receptor through leap ensures a reasonable starting structure and can help identify obvious issues sooner rather than later.

At this point the .parm7 and .crd have been created via tleap ambpdb can be used to obtain the clean pdb 4qmz_rec_leap.pdb 

    ambpdb -p 4qmz_rec_leap.parm7 -c 4qmz_rec_leap.crd > 4qmz_rec_leap.pdb

Now add partial charges to the receptor and save file in .mol2 format: 

    open chimera
    load 4qmz_rec_leap.pdb
    Tools --> Structure editing --> Add charge --> AMBER ff99SB with AM1-BCC charges 
    File --> Save Mol2... --> 4qmz.rec.mol2

A no-hydrogen receptor pdb file will now be created: 

    Chimera, load 4qmz.rec.mol2
    Select --> Chemistry --> element --> H
    Actions --> Atoms --> Delete
    File --> Save PDB... --> 4qmz.rec.noH.pdb

Ligand (4qmz_lig.pdb) will now be charged and saved in mol2 format 

    open chimera
    load 4qmz_lig.pdb
    Tools --> structure editing --> AddH
    Tools --> Structure editing --> Add charge --> AMBER ff99SB with AM1-BCC charges 
    File --> Save mol2 --> 4qmz.lig.mol2

Placement of partial charges can be verified by examining the saved files 4qmz.lig.mol2 and 4qmz.rec.mol2.

III. Generating Receptor Surface and Spheres

Open 4qmz.rec.noH.pdb in Chimera

To generate the molecular surface: 

  Action --> Surface --> Show


Save the .dms file 

  Tools --> Structure editing --> write DMS

.dms save to 4qmz.rec.noH.dms



Create surface spheres

create input file INSPH 

    4qmz.rec.noH.dms     
    R
    X
    0.0
    4.0
    1.4
    4qmz.rec.sph

line 1 designates input file line 2 designates the generated spheres will be outside the receptor surface line 3 designates that all points on the receptor will be used line 4 designates the maximum surface radius of the spheres line 5 designates the minimum surface radius of the spheres line 6 designates the output file name

run sphgen 

    sphgen -i INSPH -o OUTSPH

sphgen is the sphere generation program from dock -i desginates the input file: INSPH -o designates the output file

At this point it is beneficial to visualize the spheres that were created. This can be done with chimera: Open chimera from the terminal choose file --> open 4qmz.rec.mol2 choose file --> open 4qmz.rec.sph

The image that appears should resemble this:

4qmz receptor with all spheres


Then we need to select the spheres pertinent to our docking experiment. Usually these spheres will be the closest N spheres to the native ligand molecule.

Run 

    sphere_selector 4qmz.rec.sph ../01.dockprep/4qmz.lig.mol2 10.0

This command will select all of the spheres within 10.0 angstroms of the ligand and output them to selected_spheres.sph

To visualize the spheres using Chimera as previously done: Launch Chimera, choose File -> Open, choose 4qmz.rec.noH.pdb File -> Open, choose output_spheres_selected.pdb Select -> Residue -> SPH Actions -> Atoms/Bonds -> sphere The selected spheres with the receptor surface should look similar to that as seen below:

4qmz receptor surface with the selected spheres within 10.0 Angstroms

IV. Generating Box and Grid

enter directory 03.box-grid

create input showbox.in 

    Y
    8.0
    ../02.surface-spheres/selected_spheres.sph
    1
    4qmz.box.pdb

This input designates to dock that: 

  we want to create a box,
  the box length should be 8.0 Angstroms,
  use the selected spheres in the file designated,
  output the box to the file specified


to use this input type: 

    showbox < showbox.in

This box can be visualized in chimera using a similar approach as visualizing the spheres


Compute the energy grid

create grid.in file 

    vi grid.in

grid.in needs to contain: 

  compute_grids yes
  grid_spacing 0.4
  output_molecule no
  contact_score no
  energy_score yes
  energy_cutoff_distance 9999
  atom_model a
  attractive_exponent 6
  repulsive_exponent 12
  distance_dielectric yes
  dielectric_factor              4
  bump_filter yes
  bump_overlap 0.75
  receptor_file ../01.dockprep/4qmz.rec.mol2
  box_file ../03.box-grid/4qmz.box.pdb
  vdw_definition_file /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
  score_grid_prefix grid


this script should output verbosely to the terminal and produce 2 output files, grid.bmp and grid.nrg, both binary files. 

  grid -i grid.in > gridinfo.out

Check the standard output file gridinfo.out to ensure that individual charges on charged residues are integer values (-1 or +1). 

  University of California at San Francisco, DOCK 4.0.1
  
  __________________Job_Information_________________
  launch_time                    Wed Feb  8 17:05:50 2017
  host_name                      unknown
  memory_limit                   -1
  working_directory              /home/campus.stonybrook.edu/ronassar/MySBFiles/projects/dock_tutorial/03.box-grid2_savegridoutput
  user_name                      ronassar
  
  ________________General_Parameters________________
  compute_grids                  yes
  grid_spacing                   0.4
  output_molecule                no
  
  ________________Scoring_Parameters________________
  contact_score                  no
  energy_score                   yes
  energy_cutoff_distance         9999
  atom_model                     a
  attractive_exponent            6
  repulsive_exponent             12
  distance_dielectric            yes
  dielectric_factor              4
  bump_filter                    yes
  bump_overlap                   0.75
  
  ____________________File_Input____________________
  receptor_file                  ../01.dockprep/4qmz.rec.mol2
  box_file                       ../03.box-grid/4qmz.box.pdb
  vdw_definition_file            /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
  
  ____________________File_Output___________________
  score_grid_prefix              grid
  
  
  
  Reading in coordinates of receptor.
  CHARGED RESIDUE MET                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE ASP                     :   -1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE ARG                     :    1.000
  CHARGED RESIDUE LYS                     :    1.000
  CHARGED RESIDUE GLU                     :   -1.000
  CHARGED RESIDUE SER                     :   -1.000
  
  Total charge on 4qmz_rec_leap.pdb       :   -5.000
  
  Reading in grid box information.
  Box center of mass                      :   -1.912    0.054   14.111
  Box dimensions                          :   39.856   41.687   38.838
  Number of grid points per side [x y z]  :      101      106       99
  Total number of grid points             :  1059894
  
  Generating scoring grids.
  Percent of protein atoms processed      :        0
  Percent of protein atoms processed      :       10
  Percent of protein atoms processed      :       20
  Percent of protein atoms processed      :       30
  Percent of protein atoms processed      :       40
  Percent of protein atoms processed      :       50
  Percent of protein atoms processed      :       60
  Percent of protein atoms processed      :       70
  Percent of protein atoms processed      :       80
  Percent of protein atoms processed      :       90
  Percent of protein atoms processed      :      100
  Writing general grid info to grid.bmp
  Writing bump grid to grid.bmp
  Writing energy grids to grid.nrg
    Writing attractive VDW energy grid
    Writing repulsive VDW energy grid
    Writing electrostatic energy grid
  
  Finished calculation.

V. Docking a Single Molecule for Pose Reproduction

Create or enter directory 4, 04.dock

Minimization

create min.in 

 vi min.in
     or
 touch min.in

min.in should contain: 

  conformer_search_type                                        rigid
  use_internal_energy                                          yes
  internal_energy_rep_exp                                      12
  internal_energy_cutoff                                       100.0
  ligand_atom_file                                             ../01.dockprep/4qmz.lig.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               yes
  use_rmsd_reference_mol                                       yes
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
  use_database_filter                                          no
  orient_ligand                                                no
  bump_filter                                                  no
  score_molecules                                              yes
  contact_score_primary                                        no
  contact_score_secondary                                      no
  grid_score_primary                                           yes
  grid_score_secondary                                         no
  grid_score_rep_rad_scale                                     1
  grid_score_vdw_scale                                         1
  grid_score_es_scale                                          1
  grid_score_grid_prefix                                       ../03.box-grid/grid
  multigrid_score_secondary                                    no
  dock3.5_score_secondary                                      no
  continuous_score_secondary                                   no
  footprint_similarity_score_secondary                         no
  pharmacophore_score_secondary                                no
  descriptor_score_secondary                                   no
  gbsa_zou_score_secondary                                     no
  gbsa_hawkins_score_secondary                                 no
  SASA_score_secondary                                         no
  amber_score_secondary                                        no
  minimize_ligand                                              yes
  simplex_max_iterations                                       1000
  simplex_tors_premin_iterations                               0
  simplex_max_cycles                                           1
  simplex_score_converge                                       0.1
  simplex_cycle_converge                                       1.0
  simplex_trans_step                                           1.0
  simplex_rot_step                                             0.1
  simplex_tors_step                                            10.0
  simplex_random_seed                                          0
  simplex_restraint_min                                        yes
  simplex_coefficient_restraint                                10.0
  atom_model                                                   all
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
  flex_defn_file                                               /opt/AMS536/dock6/parameters/flex.defn
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
  ligand_outfile_prefix                                        4qmz.lig.min
  write_orientations                                           no
  num_scored_conformers                                        1
  rank_ligands                                                 no

run: 

  dock6 -i min.in

this command will output 4qmz.lig.min_scored.mol2

This structure can be visualized using chimera and loading the surface into chimera along with the minimized and unminimized ligand structures:

4qmz with ligand before and after minimization

The RMSD of the minimized ligand to crystal structure ligand is 0.297 A.


At this point we are going to calculate the van der waals and electrostatic footprints for the unminimized and minimized ligand structures in relation to the receptor active site.

create the input file: footprint.in 

  conformer_search_type                                        rigid
  use_internal_energy                                          no
  ligand_atom_file                                             ./4qmz.lig.min_scored.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               no
  use_database_filter                                          no
  orient_ligand                                                no
  bump_filter                                                  no
  score_molecules                                              yes
  contact_score_primary                                        no
  contact_score_secondary                                      no
  grid_score_primary                                           no
  grid_score_secondary                                         no
  multigrid_score_primary                                      no
  multigrid_score_secondary                                    no
  dock3.5_score_primary                                        no
  dock3.5_score_secondary                                      no
  continuous_score_primary                                     no
  continuous_score_secondary                                   no
  footprint_similarity_score_primary                           yes
  footprint_similarity_score_secondary                         no
  fps_use_footprint_reference_mol2                             yes
  fps_footprint_reference_mol2_filename                        ../01.dockprep/4qmz.lig.mol2
  fps_foot_compare_type                                        Euclidean
  fps_normalize_foot                                           no
  fps_foot_comp_all_residue                                    yes
  fps_receptor_filename                                        ../01.dockprep/4qmz.rec.mol2
  fps_vdw_att_exp                                              6
  fps_vdw_rep_exp                                              12
  fps_vdw_rep_rad_scale                                        1
  fps_use_distance_dependent_dielectric                        yes
  fps_dielectric                                               4.0
  fps_vdw_fp_scale                                             1
  fps_es_fp_scale                                              1
  fps_hb_fp_scale                                              0
  pharmacophore_score_secondary                                no
  descriptor_score_secondary                                   no
  gbsa_zou_score_secondary                                     no
  gbsa_hawkins_score_secondary                                 no
  SASA_score_secondary                                         no
  amber_score_secondary                                        no
  minimize_ligand                                              no
  atom_model                                                   all
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
  flex_defn_file                                               /opt/AMS536/dock6/parameters/flex.defn
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
  ligand_outfile_prefix                                        fps.min.output
  write_footprints                                             yes
  write_hbonds                                                 yes
  write_orientations                                           no
  num_scored_conformers                                        1
  rank_ligands                                                 no


this footprint calculation can be run in the same way as the minimization: 

  dock6 -i footprint.in

This will output three files: fps.min.output_scored.mol2 

                             fps.min.output_scored_footprint_scored.txt
                             fps.min.output_hbond_scored.txt

With an in house script, plot_footprint_single_magnitude.py the two sets of footprints will be plotted for comparison

this script will be run using: 

  python plot_footprint_single_magnitude.py fps.min.output_scored_footprint_scored.txt 50

The output will look something like this:

Footprint comparison between unminimized and minimized ligand structure

Notice the vdW clashes are now reduced (vdW positive peaks disappeared).


Rigid Docking

To help identify any errors in preparing the ligand and receptor files, a rigid docking, fixed growth docking and flexible docking runs will be performed. This is a check that the program can put the ligand into the crystal structure pose before proceeding.

Create an input file 'rigid.in' for the rigid docking run with the following parameters: 

  conformer_search_type                                        rigid
  use_internal_energy                                          yes
  internal_energy_rep_exp                                      12
  internal_energy_cutoff                                       100.0
  ligand_atom_file                                             ../01.dockprep/4qmz.lig.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               yes
  use_rmsd_reference_mol                                       yes
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
  use_database_filter                                          no
  orient_ligand                                                yes
  automated_matching                                           yes
  receptor_site_file                                           ../02.surface-spheres/selected_spheres.sph
  max_orientations                                             1000
  critical_points                                              no
  chemical_matching                                            no
  use_ligand_spheres                                           no
  bump_filter                                                  no
  score_molecules                                              yes
  contact_score_primary                                        no
  contact_score_secondary                                      no
  grid_score_primary                                           yes
  grid_score_secondary                                         no
  grid_score_rep_rad_scale                                     1
  grid_score_vdw_scale                                         1
  grid_score_es_scale                                          1
  grid_score_grid_prefix                                       ../03.box-grid/grid
  multigrid_score_secondary                                    no
  dock3.5_score_secondary                                      no
  continuous_score_secondary                                   no
  footprint_similarity_score_secondary                         no
  pharmacophore_score_secondary                                no
  descriptor_score_secondary                                   no
  gbsa_zou_score_secondary                                     no
  gbsa_hawkins_score_secondary                                 no
  SASA_score_secondary                                         no
  amber_score_secondary                                        no
  minimize_ligand                                              yes
  simplex_max_iterations                                       1000
  simplex_tors_premin_iterations                               0
  simplex_max_cycles                                           1
  simplex_score_converge                                       0.1
  simplex_cycle_converge                                       1.0
  simplex_trans_step                                           1.0
  simplex_rot_step                                             0.1
  simplex_tors_step                                            10.0
  simplex_random_seed                                          0
  simplex_restraint_min                                        no
  atom_model                                                   all
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
  flex_defn_file                                               /opt/AMS536/dock6/parameters/flex.defn
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
  ligand_outfile_prefix                                        4qmz.rigid_norestraint
  write_orientations                                           no
  num_scored_conformers                                        50
  write_conformations                                          no
  cluster_conformations                                        yes
  cluster_rmsd_threshold                                       2.0
  rank_ligands                                                 no

Outputs the file "4qmz.rigid_norestraint" the no restraint refers to the option of the parameter simplex_restraint_min in the input file. This removes the tether when doing the minimization. The best scoring structure had an RMSD of 1.??? to the native ligand pose.

Rigid docking results

Fixed Anchor Docking

In fixed anchor docking, the ligand location is fixed in the active site and the rotatable bonds space is searched (i.e angles are flexible).

Create an input file 'fad.in' for the rigid docking run with the following parameters: Notice the search algorithm is now rigidset to flex (for flexible) and orient option is set to no (for fixing the anchor in place). 

  conformer_search_type                                        flex
  user_specified_anchor                                        no
  limit_max_anchors                                            no
  min_anchor_size                                              5
  pruning_use_clustering                                       yes
  pruning_max_orients                                          1000
  pruning_clustering_cutoff                                    100
  pruning_conformer_score_cutoff                               100.0
  pruning_conformer_score_scaling_factor                       1.0
  use_clash_overlap                                            no
  write_growth_tree                                            no
  use_internal_energy                                          yes
  internal_energy_rep_exp                                      12
  internal_energy_cutoff                                       100.0
  ligand_atom_file                                             ../01.dockprep/4qmz.lig.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               yes
  use_rmsd_reference_mol                                       yes
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
  use_database_filter                                          no
  orient_ligand                                                no
  automated_matching                                           yes
  receptor_site_file                                           ../02.surface-spheres/selected_spheres.sph
  max_orientations                                             1000
  critical_points                                              no
  chemical_matching                                            no
  use_ligand_spheres                                           no
  bump_filter                                                  no
  score_molecules                                              yes
  contact_score_primary                                        no
  contact_score_secondary                                      no
  grid_score_primary                                           yes
  grid_score_secondary                                         no
  grid_score_rep_rad_scale                                     1
  grid_score_vdw_scale                                         1
  grid_score_es_scale                                          1
  grid_score_grid_prefix                                       ../03.box-grid/grid
  multigrid_score_secondary                                    no
  dock3.5_score_secondary                                      no
  continuous_score_secondary                                   no
  footprint_similarity_score_secondary                         no
  pharmacophore_score_secondary                                no
  descriptor_score_secondary                                   no
  gbsa_zou_score_secondary                                     no
  gbsa_hawkins_score_secondary                                 no
  SASA_score_secondary                                         no
  amber_score_secondary                                        no
  minimize_ligand                                              yes
  minimize_anchor                                              yes
  minimize_flexible_growth                                     yes
  use_advanced_simplex_parameters                              no
  simplex_max_cycles                                           1
  simplex_score_converge                                       0.1
  simplex_cycle_converge                                       1.0
  simplex_trans_step                                           1.0
  simplex_rot_step                                             0.1
  simplex_tors_step                                            10.0
  simplex_anchor_max_iterations                                500
  simplex_grow_max_iterations                                  500
  simplex_grow_tors_premin_iterations                          0
  simplex_random_seed                                          0
  simplex_restraint_min                                        no
  atom_model                                                   all
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
  flex_defn_file                                               /opt/AMS536/dock6/parameters/flex.defn
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
  ligand_outfile_prefix                                        4qmz.fad
  write_orientations                                           no
  num_scored_conformers                                        100
  write_conformations                                          no
  cluster_conformations                                        yes
  cluster_rmsd_threshold                                       2.0
  rank_ligands                                                 no
Fixed anchor docking results

Flexible Docking

In flexible docking, the ligand location and rotatable bonds are search.

Create an input file 'flex.in' for the flexible docking run with the following parameters: Notice the search algorithm is now set to flex (for flexible) and orient option is set to yes (for fixing the anchor in place). 

  conformer_search_type                                        flex
  user_specified_anchor                                        no
  limit_max_anchors                                            no
  min_anchor_size                                              5
  pruning_use_clustering                                       yes
  pruning_max_orients                                          1000
  pruning_clustering_cutoff                                    100
  pruning_conformer_score_cutoff                               100.0
  pruning_conformer_score_scaling_factor                       1.0
  use_clash_overlap                                            no
  write_growth_tree                                            no
  use_internal_energy                                          yes
  internal_energy_rep_exp                                      12
  internal_energy_cutoff                                       100.0
  ligand_atom_file                                             ../01.dockprep/4qmz.lig.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               yes
  use_rmsd_reference_mol                                       yes
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
  use_database_filter                                          no
  orient_ligand                                                yes
  automated_matching                                           yes
  receptor_site_file                                           ../02.surface-spheres/selected_spheres.sph
  max_orientations                                             1000
  critical_points                                              no
  chemical_matching                                            no
  use_ligand_spheres                                           no
  bump_filter                                                  no
  score_molecules                                              yes
  contact_score_primary                                        no
  contact_score_secondary                                      no
  grid_score_primary                                           yes
  grid_score_secondary                                         no
  grid_score_rep_rad_scale                                     1
  grid_score_vdw_scale                                         1
  grid_score_es_scale                                          1
  grid_score_grid_prefix                                       ../03.box-grid/grid
  multigrid_score_secondary                                    no
  dock3.5_score_secondary                                      no
  continuous_score_secondary                                   no
  footprint_similarity_score_secondary                         no
  pharmacophore_score_secondary                                no
  descriptor_score_secondary                                   no
  gbsa_zou_score_secondary                                     no
  gbsa_hawkins_score_secondary                                 no
  SASA_score_secondary                                         no
  amber_score_secondary                                        no
  minimize_ligand                                              yes
  minimize_anchor                                              yes
  minimize_flexible_growth                                     yes
  use_advanced_simplex_parameters                              no
  simplex_max_cycles                                           1
  simplex_score_converge                                       0.1
  simplex_cycle_converge                                       1.0
  simplex_trans_step                                           1.0
  simplex_rot_step                                             0.1
  simplex_tors_step                                            10.0
  simplex_anchor_max_iterations                                500
  simplex_grow_max_iterations                                  500
  simplex_grow_tors_premin_iterations                          0
  simplex_random_seed                                          0
  simplex_restraint_min                                        no
  atom_model                                                   all
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
  flex_defn_file                                               /opt/AMS536/dock6/parameters/flex.defn
  flex_drive_file                                              /opt/AMS536/dock6/parameters/flex_drive.tbl
  ligand_outfile_prefix                                        4qmz.flex
  write_orientations                                           no
  num_scored_conformers                                        100
  write_conformations                                          no
  cluster_conformations                                        yes
  cluster_rmsd_threshold                                       2.0
  rank_ligands                                                 no
Flexible docking results

VI. Virtual Screening

At this point we will begin docking potential inhibitors into the active site of 4qmz. The compound library we will use contains 25,000 molecules and was previously prepared for another docking experiment. The compound library contains molecules with rotatable bonds number between 7 and 15. The library used is titled small_ligand_library.mol2, this library can be found at /gpfs/projects/AMS536/virtual_screen_library on seawulf.

We began by copying all of our local directories outlined in the first section to a project directory on seawulf to allow use to expedite the docking process over the 25,000 compounds in our test library.

We will create a virtual-screen.in file that looks like the following: 

  conformer_search_type                                        flex
  user_specified_anchor                                        no
  limit_max_anchors                                            no
  min_anchor_size                                              5
  pruning_use_clustering                                       yes
  pruning_max_orients                                          1000
  pruning_clustering_cutoff                                    100
  pruning_conformer_score_cutoff                               100.0
  pruning_conformer_score_scaling_factor                       1.0
  use_clash_overlap                                            no
  write_growth_tree                                            no
  use_internal_energy                                          yes
  internal_energy_rep_exp                                      12
  internal_energy_cutoff                                       100.0
  ligand_atom_file                                             small_ligand_library.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               yes
  use_rmsd_reference_mol                                       yes
  rmsd_reference_filename                                      ../01.dockprep/4qmz.lig.mol2
  use_database_filter                                          no
  orient_ligand                                                yes
  automated_matching                                           yes
  receptor_site_file                                           ../02.surface-spheres/selected_spheres.sph
  max_orientations                                             1000
  critical_points                                              no
  chemical_matching                                            no
  use_ligand_spheres                                           no
  bump_filter                                                  no
  score_molecules                                              yes
  contact_score_primary                                        no
  contact_score_secondary                                      no
  grid_score_primary                                           yes
  grid_score_secondary                                         no
  grid_score_rep_rad_scale                                     1
  grid_score_vdw_scale                                         1
  grid_score_es_scale                                          1
  grid_score_grid_prefix                                       ../03.box-grid/grid
  multigrid_score_secondary                                    no
  dock3.5_score_secondary                                      no
  continuous_score_secondary                                   no
  footprint_similarity_score_secondary                         no
  pharmacophore_score_secondary                                no
  descriptor_score_secondary                                   no
  gbsa_zou_score_secondary                                     no
  gbsa_hawkins_score_secondary                                 no
  SASA_score_secondary                                         no
  amber_score_secondary                                        no
  minimize_ligand                                              yes
  minimize_anchor                                              yes
  minimize_flexible_growth                                     yes
  use_advanced_simplex_parameters                              no
  simplex_max_cycles                                           1
  simplex_score_converge                                       0.1
  simplex_cycle_converge                                       1.0
  simplex_trans_step                                           1.0
  simplex_rot_step                                             0.1
  simplex_tors_step                                            10.0
  simplex_anchor_max_iterations                                500
  simplex_grow_max_iterations                                  500
  simplex_grow_tors_premin_iterations                          0
  simplex_random_seed                                          0
  simplex_restraint_min                                        no
  atom_model                                                   all
  vdw_defn_file                                            /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn
  flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn
  flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl
  ligand_outfile_prefix                                        4qmz.virtualscreen
  write_orientations                                           no
  num_scored_conformers                                        1
  write_conformations                                          no
  cluster_conformations                                        yes
  cluster_rmsd_threshold                                       2.0
  rank_ligands                                                 no

This script is submitted to the seawulf cluster using a script called screen.sh which contains the following: 

  #! /bin/tcsh
  #PBS -l nodes=4:ppn=28
  #PBS -l walltime=24:00:00
  #PBS -o jobstatus.out
  #PBS -e joberror.err
  #PBS -N large_vs
  #PBS -V
  #PBS -q long
  cd $PBS_O_WORKDIR
  mpirun -np 112 dock6.mpi -i virtual-screen.in -o virtualscreen1.out


This screen.sh input specifies that we will use 112 processors to dock our virtual library. We have a wall clock for the job of 24 hours and we will output jobstatus.out, joberror.err and virtualscreen1.out. In addition to these output files there is an output for the work performed by each processor individually.

Cartesian Minimization of the docked molecules

Now we will do cartesian minimization of the docked molecules. Cartesian here is referring to the receptor being explicitly present during the minimization.

Write the minimization input script minimize-vs.in: 

  conformer_search_type                                        rigid
  use_internal_energy                                          yes
  internal_energy_rep_exp                                      12
  internal_energy_cutoff                                       100.0
  ligand_atom_file                                             4qmz.virtualscreen_scored.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               no
  use_database_filter                                          no
  orient_ligand                                                no
  bump_filter                                                  no
  score_molecules                                              yes
  contact_score_primary                                        no
  contact_score_secondary                                      no
  grid_score_primary                                           no
  grid_score_secondary                                         no
  multigrid_score_primary                                      no
  multigrid_score_secondary                                    no
  dock3.5_score_primary                                        no
  dock3.5_score_secondary                                      no
  continuous_score_primary                                     yes
  continuous_score_secondary                                   no
  cont_score_rec_filename                                      ../01.dockprep/4qmz.rec.mol2
  cont_score_att_exp                                           6
  cont_score_rep_exp                                           12
  cont_score_rep_rad_scale                                     1
  cont_score_use_dist_dep_dielectric                           yes
  cont_score_dielectric                                        4.0
  cont_score_vdw_scale                                         1
  cont_score_es_scale                                          1
  footprint_similarity_score_secondary                         no
  pharmacophore_score_secondary                                no
  descriptor_score_secondary                                   no
  gbsa_zou_score_secondary                                     no
  gbsa_hawkins_score_secondary                                 no
  SASA_score_secondary                                         no
  amber_score_secondary                                        no
  minimize_ligand                                              yes
  simplex_max_iterations                                       1000
  simplex_tors_premin_iterations                               0
  simplex_max_cycles                                           1
  simplex_score_converge                                       0.1
  simplex_cycle_converge                                       1.0
  simplex_trans_step                                           1.0
  simplex_rot_step                                             0.1
  simplex_tors_step                                            10.0
  simplex_random_seed                                          0
  simplex_restraint_min                                        no
  atom_model                                                   all
  vdw_defn_file                                                /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn
  flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn
  flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl
  ligand_outfile_prefix                                        4qmz.virtualscreen1.minimized
  write_orientations                                           no
  num_scored_conformers                                        1
  rank_ligands                                                 no

Submit the minimization job using with the script "job-min.sh" following lines: 

  #!/bin/tcsh
  #PBS -l nodes=16:ppn=28
  #PBS -l walltime=12:00:00
  #PBS -o jobstatusmin.out
  #PBS -e jobstatusmin.err
  #PBS -N chunk0min
  #PBS -V
  #PBS -q medium
  module load intel/mpi/64/2017/0.098 
  cd $PBS_O_WORKDIR
  mpirun -np 448 /gpfs/projects/AMS536/zzz.programs/dock6/bin/dock6.mpi -i minimize-vs.in -o minimize-vs.out


At this point it would be good to rank our docked ligands by score and extract the 100 best scoring (100 lowest scoring, most negative) ligands to investigate further.


Rescoring Docked Molecules

Dock has several built in scoring functions. In an effort to utilize these alternative scoring functions it is worthwhile to rescore the results from the initial virtual screen that were collected. The virtual screen will be rescored using footprint similarity, pharmacophore score, tanimoto score, the hungarian and the volume overlap score. Each of these scoring functions focuses in on a different chemical aspect of the ligand in question. The rescore.in script for this run can be found here: 

conformer_search_type                                        rigid
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
internal_energy_cutoff                                       100.0
ligand_atom_file                            ../05.virtual- screen/4qmz.virtualscreen1.minimized_scored.mol2 
limit_max_ligands                                            no
skip_molecule                                                no
read_mol_solvation                                           no
calculate_rmsd                                               no
use_database_filter                                          no
orient_ligand                                                no
bump_filter                                                  no
score_molecules                                              yes
contact_score_primary                                        no
contact_score_secondary                                      no
grid_score_primary                                           no
grid_score_secondary                                         no
multigrid_score_primary                                      no
multigrid_score_secondary                                    no
dock3.5_score_primary                                        no
dock3.5_score_secondary                                      no
continuous_score_primary                                     no
continuous_score_secondary                                   no
footprint_similarity_score_primary                           no
footprint_similarity_score_secondary                         no
pharmacophore_score_primary                                  no
pharmacophore_score_secondary                                no
descriptor_score_primary                                     yes
descriptor_score_secondary                                   no
descriptor_use_grid_score                                    no
descriptor_use_multigrid_score                               no
descriptor_use_continuous_energy                             no
descriptor_use_footprint_similarity                          yes
descriptor_use_pharmacophore_score                           yes
descriptor_use_tanimoto                                      yes
descriptor_use_hungarian                                     yes
descriptor_use_volume_overlap                                yes
descriptor_fps_use_footprint_reference_mol2                  yes
descriptor_fps_footprint_reference_mol2_filename             ../04.dock/4qmz.lig.min_scored.mol2
descriptor_fps_foot_compare_type                             Euclidean
descriptor_fps_normalize_foot                                no
descriptor_fps_foot_comp_all_residue                         yes
descriptor_fps_receptor_filename                             ../01.dockprep/4qmz.rec.mol2
descriptor_fps_vdw_att_exp                                   6
descriptor_fps_vdw_rep_exp                                   12
descriptor_fps_vdw_rep_rad_scale                             1
descriptor_fps_use_distance_dependent_dielectric             yes
descriptor_fps_dielectric                                    4.0
descriptor_fps_vdw_fp_scale                                  1
descriptor_fps_es_fp_scale                                   1
descriptor_fps_hb_fp_scale                                   0
descriptor_fms_score_use_ref_mol2                            yes
descriptor_fms_score_ref_mol2_filename                       ../04.dock/4qmz.lig.min_scored.mol2
descriptor_fms_score_write_reference_pharmacophore_mol2      no
descriptor_fms_score_write_reference_pharmacophore_txt       no
descriptor_fms_score_write_candidate_pharmacophore           no
descriptor_fms_score_write_matched_pharmacophore             no
descriptor_fms_score_compare_type                            overlap
descriptor_fms_score_full_match                              yes
descriptor_fms_score_match_rate_weight                       5.0
descriptor_fms_score_match_dist_cutoff                       1.0
descriptor_fms_score_match_proj_cutoff                       0.7071
descriptor_fms_score_max_score                               20
descriptor_fingerprint_ref_filename                          ../04.dock/4qmz.lig.min_scored.mol2
descriptor_hungarian_ref_filename                            ../04.dock/4qmz.lig.min_scored.mol2
descriptor_hungarian_matching_coeff                          -5
descriptor_hungarian_rmsd_coeff                              1
descriptor_volume_reference_mol2_filename                    ../04.dock/4qmz.lig.min_scored.mol2
descriptor_volume_overlap_compute_method                     analytical
descriptor_weight_fps_score                                  1
descriptor_weight_pharmacophore_score                        1
descriptor_weight_fingerprint_tanimoto                       -1
descriptor_weight_hungarian                                  1
descriptor_weight_volume_overlap_score                       -1
gbsa_zou_score_secondary                                     no
gbsa_hawkins_score_secondary                                 no
SASA_score_secondary                                         no
amber_score_secondary                                        no
minimize_ligand                                              no
atom_model                                                   all
vdw_defn_file                     /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn
flex_defn_file                    /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn
flex_drive_file                   /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl
chem_defn_file                    /gpfs/projects/AMS536/zzz.programs/dock6/parameters/chem.defn
pharmacophore_defn_file           /gpfs/projects/AMS536/zzz.programs/dock6/parameters/ph4.defn
ligand_outfile_prefix                                        descriptor.output
write_footprints                                             yes
write_hbonds                                                 yes
write_orientations                                           no
num_scored_conformers                                        1
rank_ligands                                                 no


From here the results should be sorted according the score based on the scoring function in question. A script for this can be written fairly quickly in whatever language you are most comfortable with. If you are not very experienced with computer programming, a python script developed in house can be found here that will search through the output .mol2 file and select the top N scoring molecules. These top scoring molecule s are taken from the virtual screen "top'N'_molid.txt" and finds and their coordinates saved to an output file "top'N'_molid.mol2".

rank_scores_get_coords.py: 

 #!/usr/bin/python
 import sys
 from operator import itemgetter, attrgetter, methodcaller
 # input mol2 file as command line argument
 f0=open(sys.argv[1],'r') ; data0 = f0.readlines() ; f0.close()                  #open mol2 data file with   scored molecules
 # define variables
 nummol = 50                                                                      # number of molecules to   keep after ranking by score
 #scoretype="Grid_Score:"                                                        # scoring function to rank  by
 scoretype="Footprint_Similarity_Score:"
 #--------------------------------------------------------------------------------------------------------------------------
 list_molid=[]                                                               #create a list to store molecule ids
 list_molscores=[]                                                           #create list to store scores for each respective molecule
 for line in data0:                                                           #reads through the input file
       if line != "\n":                                                    #if the line is not empty  continue
               if line.split()[0]=="##########":                           #if the line starts with ############ then continue
                               if line.split()[1]=="Name:":                #if the line has the next entry as Name: continue
                                       molid = line.split()[2]             #then add the molID to list_molid
                                       list_molid.append(molid)
                               elif line.split()[1]==scoretype:
                                       molscore = line.split()[2]
                                       list_molscores.append(molscore)
 #print list_molid
 #print list_molscores
 both = zip(list_molid,list_molscores)                                      #concatenates list_molid and list_molscores into one array
 sorted_both=sorted(both, key=itemgetter(1), reverse=True)                  #sorts the "both" list
 sorted_top=[]                                                             #names a new list sorted_top5
 #for element in sorted_both:                                                #print out first n entries in sorted_both
 #       print element[0], element[1]
 o=open("top"+str(nummol)+"_molid.txt","w")                                 #opens output file topx_molid.txt
 for element in sorted_both[:nummol]:                                            
       sorted_top.append(element[0])                                           
       o.write(element[0]+"\n")                                           #stores the top nummol entries in sorted_top
 o.close()
 #print sorted_top                                                          
 ##
 ##--------------------------------------------------------------------------------------------------------------------------
 ##
 ## Get coordinates of molecules of interest from scored mol2 
 f1=open("top"+str(nummol)+"_molid.txt","r") ; data1=f1.readlines(); f1.close() # takes in sorted list of top molecules
 f2=open(sys.argv[1],"r") ; data2=f2.readlines(); f2.close()                    # takes in mol2 results 
 identifier=0    # this will switch between 0 and 1. When it is 0 no line printing, when it is 1 we print the line.
 molid_top=[]  # list of molecule ids that are of interest 
 # go through the first file (has top n molecules) and append those ids to molecule ids list
 for line in data1:   
       molid=line.split()[0]
       molid_top.append(molid)
 output=open("top"+str(nummol)+"_molid.mol2","w")
 # go through the second file (the rescored mol2 file), the next commands print out the coordinates for the molecules of interest
 for line in data2:
 #       if line != "\n":  # skipping empty lines
               if identifier==1:
                       try: 
                               if line.split()[0]=="##########" and line.split()[1]=="Name:" and line.split()[2] not in molid_top:   # changing the identifier when hitting a new molecule that is not of interest
                                       identifier=0
                                       output.write('\n')      
                               else:
                                       output.write(line)   # else keep the identifier at 1 and keep printing
                       except IndexError:
                               continue  
               try:    
                       if identifier==0 and line.split()[0]=="##########" and line.split()[1]=="Name:" and line.split()[2] in molid_top:  # change identifier to 1 when we hit molecule of interest
                               identifier=1
                               output.write(line)
               except IndexError:
                               continue
 output.close()

VIII. Frequently Encountered Problems

Spending the time to ensure your starting structure is reasonable is of the utmost importance and will save you much headache moving forward in your study. Reading through the PDB file as downloaded from the protein data bank will give the researcher a lot of information regarding issues or grey area encountered by the researchers responsible for publishing the structure initially. Threading a structure through some quality control software such as tleap, chimera or another comparable software package will ensure charges are reasonable, hydrogens are appropriately placed, and there are no missing residues in the peptide chain.

A basic understanding of your enzymes active site will also be beneficial while approaching an enzymatic system from a docking perspective. Some questions to ask oneself are: does my protein have an ion (Mg, Ca, etc.) in the active site? Does my ligand make important interactions with conserved water residues? If so these waters could be kept in with the receptor to aid in discovery. The alternative would be to design an inhibitor that displaces the water molecules and supplements favorable interactions mediated by the water between the ligand and protein.

Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2017_DOCK_tutorial_1_with_PDB_4QMZ_NEW&oldid=10281"