Difference between revisions of "2018 DOCK tutorial 2 with PDBID 1C87"

From Rizzo_Lab
Jump to: navigation, search
m (→‎Flexible Docking)
 
(28 intermediate revisions by one other user not shown)
Line 15: Line 15:
  
 
===Organizing Directories===
 
===Organizing Directories===
−
We are going to create and organize directories so it would be easier for us to find or identify files in each directory.  
+
We are going to create and organize directories so it would be easier for us to find or identify files in each directory. The structure of your directories should look something like this:
  
 
{|
 
{|
 
  |<pre>~/gpfs/projects/AMS536/2018/
 
  |<pre>~/gpfs/projects/AMS536/2018/
 
                           /001.files/
 
                           /001.files/
−
                           /002.spheres/
+
                           /002.surface-spheres/
 
                           /003.box/
 
                           /003.box/
 
                           /004.grid/
 
                           /004.grid/
Line 29: Line 29:
 
  |
 
  |
 
  |}
 
  |}
 +
 +
 +
To locate which directory you are currently in, type '''pwd'''. From your personal directory (...../AMS536/2018/yourName/), you can use the '''mkdir''' command to make directories. One way to make a number of directories is:
 +
 +
 +
  mkdir 001.files 002.spheres 003.box 004.grid 005.dock 006.virtual-screen
 +
 +
To see if you were successful, type
 +
 +
  ls -l
 +
 +
to list directory contents. If you make a mistake, you can change the name of the directory or remove it and try again
 +
 +
to rename a directory..
 +
  mv oldname/ newname/
 +
 +
to remove a directory (must be empty)
 +
  rmdir mydirectory/
  
 
==II. Preparing the Receptor and Ligand==   
 
==II. Preparing the Receptor and Ligand==   
Line 36: Line 54:
 
We are going to the PDB website (https://www.rcsb.org/) to download 1C87.pdb file and transfer this pdb file to your directory.  
 
We are going to the PDB website (https://www.rcsb.org/) to download 1C87.pdb file and transfer this pdb file to your directory.  
 
First, open Chimera and load 1C87.pdb file.   
 
First, open Chimera and load 1C87.pdb file.   
−
Remove the ligand and save the receptor in mol2.format. "mol2" format shows types of bonds whether it is single or double bond.
+
For dock, it is important that the receptor is separate entity and saved as its own structure file. We must also think about what parts of the structure should be included in the receptor and what parts should be ignored. (We may want some crystallographic water which interact with the binding site/ligand but may choose to remove other surrounding waters further away.
−
After saving the receptor with no H (''noh_receptor_1c87.mol2'') and the receptor with H (''rec_withH_1c87.mol2''), load the original 1c87.pdb again. Remove the receptor and water, and save the ligand with H (''lig_withH_1c87.mol2'') and ligand with H and charge (''lig_withH_charge_1c87.mol2'').  
+
Remove the ligand and save the receptor in mol2.format. "mol2" format shows the bond types (whether it is single or double bond), and the connectivity of the atoms (what is bonded to what).
−
In order to add charge on ligand, go to '''Tools''' and '''Structure Editing''' and save the files. '''(NOTE: when you add H on ligand, make the charge -1. Check the article that is related to pdb file to decide whether it should be pronated or depronated.)'''  
+
After saving the receptor with no Hydrogen (''noh_receptor_1c87.mol2'') and the receptor with built Hydrogens and partial charges (''rec_withH_charge_1c87.mol2''), load the original 1c87.pdb again. Remove the receptor and water, and save the ligand with no Hydrogen (''noh_lig_1c87.mol2'') and ligand with Hydrogen and partial charges (''lig_withH_charge_1c87.mol2'').  
 +
In order to build Hydrogens and add charge on these structures, go to '''Tools''' and '''Structure Editing''' and use '''Dockprep'''. Along with building Hydrogens and adding partial charges, dockprep will try to fill in missing info from the pdb and also save the structures in a mol2 format at the end . '''(NOTE: For this system, when you add H on ligand, make the charge -1. In other cases, check the article that is related to pdb file to decide how it should be pronated or depronated.)'''  
 +
 
 +
Go to
 +
  Tools -> Structure Editing -> Dockprep
 +
Check the appropriate boxes:
  
 +
(### Need photo of dockprep setup ###)
 +
 +
Or more explicitly:
 
  Tools -> Structure Editing -> Add H
 
  Tools -> Structure Editing -> Add H
 
  Tools -> Structure Editing-> Add Charge
 
  Tools -> Structure Editing-> Add Charge
Line 45: Line 71:
 
  File -> Save Mol2... -> lig_withH_charge.mol2
 
  File -> Save Mol2... -> lig_withH_charge.mol2
  
−
Lastly, open the receptor file and click '''Surface Editing''' and write DMS file. Make sure to save all files.   
+
Lastly, open the receptor file and tell chimera to show the surface of your receptor. After this is calculated and you can see the surface of the protein, go to '''Structure Editing''' and '''write DMS''' file. Make sure to save all files.   
 
So far, we will have these files ready.
 
So far, we will have these files ready.
 
{|
 
{|
 
  |<pre>
 
  |<pre>
−
1c87.pdb           lig_withH_charge_1c87.mol2    noh_receptor_1c87.mol2
+
  1c87.pdb             lig_withH_charge_1c87.mol2    noh_receptor_1c87.mol2
−
lig_withH_1c87.mol2    noh_lig_1c87.mol2          rec_withH_1c87.mol2
+
  noh_lig_1c87.mol2          rec_withH_1c87.mol2     rec_1c87.dms
 
  |
 
  |
 
  |}
 
  |}
Line 61: Line 87:
 
=== Generating the Receptor Surface ===
 
=== Generating the Receptor Surface ===
 
Make sure sphere directory is created and open the directory:
 
Make sure sphere directory is created and open the directory:
−
  mkdir 002.surface
+
  mkdir 002.surface-spheres
−
  cd 002.surface
+
  cd 002.surface-spheres
  
 
Open Chimera and load the receptor surface:
 
Open Chimera and load the receptor surface:
Line 95: Line 121:
 
Load ''1c87.recep.sph'' file.  
 
Load ''1c87.recep.sph'' file.  
  
−
[[Image:Generated1c87spheres.png|thumb|center|1000px| The generated 1C87 surface spheres]]
+
[[Image:Generated speres and receptor.png|thumb|center|1000px| The generated 1C87 surface spheres with a receptor]]
  
 
Now, we are going to select the spheres that is close to the native ligand molecule, which would be within 8.0 angstroms of the ligand.  
 
Now, we are going to select the spheres that is close to the native ligand molecule, which would be within 8.0 angstroms of the ligand.  
Line 101: Line 127:
  
 
''1c87.recep.pdb''will be created and load this file in Chimera.
 
''1c87.recep.pdb''will be created and load this file in Chimera.
 +
[[Image:Selected spheres and receptor.png|thumb|center|1000px| The generated 1C87 surface spheres with a receptor]]
  
 
==IV. Generating Box and Grid==
 
==IV. Generating Box and Grid==
  
 +
===Creating the box===
 
Enter ''003.box-grid'' directory.
 
Enter ''003.box-grid'' directory.
 
Create input'' showbox.in'':
 
Create input'' showbox.in'':
Line 118: Line 146:
 
  showbox < showbox.in
 
  showbox < showbox.in
 
Now we can visualize the box in Chimera. Load ''1c87.box.pdb''.
 
Now we can visualize the box in Chimera. Load ''1c87.box.pdb''.
 +
 
[[Image:spheres_in_grid.png|thumb|center|1000px| The selected spheres with box]]
 
[[Image:spheres_in_grid.png|thumb|center|1000px| The selected spheres with box]]
 +
[[File:Box-sphere.png|thumb|center|1000px| Receptor w/ spheres in box]]
  
−
 
+
===Calculating Energy Grid===
−
Then, we will compute the energy grid. Create ''grid.in'' file
+
Now, we will compute the energy grid. Create ''grid.in'' file
 +
This saves the interaction energy of the receptor at each grid ''point''
  
 
  vi grid.in
 
  vi grid.in
Line 142: Line 173:
 
   bump_filter              yes
 
   bump_filter              yes
 
   bump_overlap              0.75
 
   bump_overlap              0.75
−
   receptor_file            ../001.files/lig_withH_charge_1c87.mol2
+
   receptor_file            ../001.files/rec_withH_charge_1c87.mol2
 
   box_file                  ../003.box/1c87.box.pdb
 
   box_file                  ../003.box/1c87.box.pdb
 
   vdw_definition_file      /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
 
   vdw_definition_file      /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
Line 148: Line 179:
  
 
Once completed to answer the questions, it will create two output files: grid.bmp and grid.nrg
 
Once completed to answer the questions, it will create two output files: grid.bmp and grid.nrg
−
 
  
 
==V. Docking a Single Molecule for Pose Reproduction==
 
==V. Docking a Single Molecule for Pose Reproduction==
 +
We would like to re-dock 1C87 into a receptor by using these methods: rigid docking, flex docking and fixed anchor docking.
  
 
===Minimization===
 
===Minimization===
  
 +
First, we are going to run minimization. Minimization is important because the structure of the ligand in the crystal may not be suited for the forcefield used in our calculations. Minimization relaxes the structure to the forcefield so that contact geometries are at an energy minimum.
 
Create ''min.in'' file and run dock with this input file.
 
Create ''min.in'' file and run dock with this input file.
  
 
  touch min.in
 
  touch min.in
−
  dock6 -i min.in
+
  ../../../zzz.programs/dock6/bin/dock6 -i min.in
  
 
Answer the questions like we did above:
 
Answer the questions like we did above:
Line 186: Line 218:
 
   continuous_score_primary                                    yes
 
   continuous_score_primary                                    yes
 
   continuous_score_secondary                                  no
 
   continuous_score_secondary                                  no
−
   cont_score_rec_filename                                      ../001.files/lig_withH_charge_1c87.mol2
+
   cont_score_rec_filename                                      ../001.files/rec_withH_charge_1c87.mol2
 
   cont_score_att_exp                                          6
 
   cont_score_att_exp                                          6
 
   cont_score_rep_exp                                          12
 
   cont_score_rep_exp                                          12
Line 213: Line 245:
 
   simplex_restraint_min                                        yes
 
   simplex_restraint_min                                        yes
 
   atom_model                                                  all
 
   atom_model                                                  all
−
  vdw_defn_file                                                /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
+
vdw_defn_file                                                /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn
−
   flex_defn_file                                              /opt/AMS536/dock6/parameters/flex.defn
+
   flex_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn
−
   flex_drive_file                                             /opt/AMS536/dock6/parameters/flex_drive.tbl
+
   flex_drive_file                                             /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl
−
   ligand_outfile_prefix                                        3pgl.min
+
  ligand_outfile_prefix 
 +
   ligand_outfile_prefix                                        1c87.min
 
   write_orientations                                          no
 
   write_orientations                                          no
 
   num_scored_conformers                                        1
 
   num_scored_conformers                                        1
 
   rank_ligands                                                no
 
   rank_ligands                                                no
  
−
Now, ''1c87.min_scored.mol2'' will be produced. We can compare minimized ligand structure with unminimized ligand structure.  
+
Now, ''1c87.min_scored.mol2'' will be produced. We can compare minimized ligand structure with unminimized ligand structure.
  
 
===Rigid Docking===
 
===Rigid Docking===
  
−
We are going to run rigid docking. In order to run dock, we will have to create ''rigid.in'' file.  
+
We are going to run rigid docking. With rigid docking the internal degrees of freedom of the ligand are constrained, and only translation and rotation of the whole structure is considered for docking.
 +
In order to run dock, we will have to create ''rigid.in'' file.  
  
 
  touch rigid.in
 
  touch rigid.in
Line 310: Line 344:
  
 
In order to run flexible docking, we will also create ''flex.in'' file.
 
In order to run flexible docking, we will also create ''flex.in'' file.
−
In this case, ligand and angles are flexible ('''NOTE: ''orient_ligand'' is set to yes''').  
+
In this case, ligand and internal angles are flexible ('''NOTE: ''orient_ligand'' is set to yes''').  
  
 
  touch flex.in
 
  touch flex.in
Line 401: Line 435:
  
  
−
We are going to run fixed anchor docking. Create ''fad.in.''The ligand is fixed in the active site and angles are flexible.  
+
We are going to run fixed anchor docking. Create ''fad.in.''The central ligand structure is fixed in the active site and internal angles are flexible and sampled.  
 
('''NOTE: ''orient_ligand'' is set to no''').  
 
('''NOTE: ''orient_ligand'' is set to no''').  
  
Line 483: Line 517:
 
  cluster_rmsd_threshold                                      2.0
 
  cluster_rmsd_threshold                                      2.0
 
  rank_ligands                                                no
 
  rank_ligands                                                no
 +
 +
Again, we can view fixed anchor docking in Chimera and find the lowest grid score value.
 +
[[Image:1c87 fad.png|thumb|center|1000px| Anchor fixed docking]]
  
 
Now we have created output files for rigid, flexible and fixed anchor docking.  
 
Now we have created output files for rigid, flexible and fixed anchor docking.  
  
−
To submit the job for rigid docking, create ''submit.sh'' file.
+
===Molecular Footprint===
 +
When it comes to examining how a ligand interacts with its receptor, we can break it down to the interaction contribution of side chains in the binding pocket. Some interactions will be much stronger than other, with some being strongly favorable or strongly prohibitive to binding. Typically we want to look at those residues which contribute most to the interaction between ligand and receptor, favorable or not, so that we can deduce ideal physiochemical structures for future lead refinement. Dock can show this in what is called a '''Footprint'''.
 +
The plot shows electrostatic and van der Waals interaction contribution for the top '''X''' highest residue interactions vs the contribution of another ligand (In this case, we want to see the footprint of a given ligand with respect to the known binder ligand 1c87 from either its crystal structure, or its minimized structure). Here we will compare the footprint of minimized 1c87 vs its crystal structure.
 +
 
 +
Create a new directory for the footprint calculations, and create the input file for dock inside that directory:
 +
 
 +
  mkdir 007.footprint
 +
  cd !$
 +
  touch footprint.in
 +
 
 +
Answer the questions to fill in the input parameters like before:
 +
 
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          no
 +
  ligand_atom_file                                            ../001.files/1c87.min.mol2 ### NOTE: This can be changes based on what you want to compare with what ###
 +
  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                        ../001.files/lig_withH_charge_1c87.mol2 ### Reference ligand for comparison of footprints ###
 +
  fps_foot_compare_type                                        Euclidean
 +
  fps_normalize_foot                                          no
 +
  fps_foot_comp_all_residue                                    yes
 +
  fps_receptor_filename                                        ../001.files/rec_withH_1c87.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                                                /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                                        fps.min.output
 +
  write_footprints                                            yes
 +
  write_hbonds                                                yes
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 
 +
This script produces a *.txt file to which you can plot the footprint using the provided python script:
 +
plot_footprint_single_magnitude.py
 +
 
 +
to run this enter the following into your command line:
 +
  python plot_footprint_single_magnitude.py fps.min.output_footprint_scored.txt 50 ### NOTE: Here, we want the 50 highest contribution interactions. The order of these arguments matter! ###
 +
 
 +
The result should look something like this:
 +
 
 +
[[File:Footprint_ex.png|thumb|center|1000px]]
 +
 
 +
==VI. Virtual Screening==
 +
 
 +
 
 +
We will create a ''virtual-screen.in'' file in ''006.virtual-screen'' 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                                      ../001.files/lig_withH_charge_1c87.mol2
 +
  use_database_filter                                          no
 +
  orient_ligand                                                yes
 +
  automated_matching                                          yes
 +
  receptor_site_file                                          ../002.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                                      ../004.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                                        1c87_virtualscreen
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  write_conformations                                          no
 +
  cluster_conformations                                        yes
 +
  cluster_rmsd_threshold                                      2.0
 +
  rank_ligands                                                no
 +
 
 +
 
 +
To submit the job for docking, create ''submit.sh'' file.
 
  vi submit.sh
 
  vi submit.sh
  
Line 500: Line 696:
 
  dock6 -i rigid.in -o 1c87_rigid_docking.out
 
  dock6 -i rigid.in -o 1c87_rigid_docking.out
  
−
After creating the script, submit the job.  
+
After creating this script, submit the job.  
 
  qsub submit.sh
 
  qsub submit.sh
 +
 +
We can view whether the job runs or does not run (if there is an error):
 +
qstat -u ''username''
 +
 +
===Large Library Virtual Screening (with MPI)===
 +
Using 24 processors for two days will not be sufficient to dock all ~25,000 compounds. We can utilize more processors using the Message Passing Interface (MPI) protocol. (This changes the contents of the submission script a little bit)
 +
 +
#! /bin/bash
 +
#PBS -l walltime=48:00:00
 +
#PBS -l '''nodes=4:ppn=28'''
 +
#PBS -q long
 +
#PBS -N 1c87_rigid
 +
#PBS -V
 +
cd $PBS_O_WORKDIR
 +
'''mpirun -np 112 dock6.mpi''' -i rigid.in -o 1c87_rigid_docking.out
 +
 +
make sure you are requesting the same number of processors on both cwulf and the MPI command (the -np flag). With these resources, we should finish all 25,000 compounds within one day
 +
 +
=== Cartesian Minimization of Docked Molecules ===
 +
Here, cartesian just means we will be using an explicit receptor environment as opposed to the energy grid calculated earlier.
 +
These calculations should be calculated in a different sub-directory so first we can create a directory called 007.cartesian-minimization, and create a new min.in file inside that directory.
 +
 +
  mkdir 008.cartesian-minimization
 +
  cd !$
 +
  touch min.in
 +
 +
Now we can call dock and answer the questions like before:
 +
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100.0
 +
  ligand_atom_file                                            1c87_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                                      ../001.files/rec_withH_1c87.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                                        1c87.virtualscreen1.minimized
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 +
 +
Since this is another job which needs lots of resources, we submit this to the computing nodes with the submission script. First, create submit.sh.
 +
 +
vi submit.sh
 +
 +
Now, copy the following into the file:
 +
 +
#! /bin/bash
 +
#PBS -l walltime=48:00:00
 +
#PBS -l nodes=4:ppn=28
 +
#PBS -q long
 +
#PBS -N 1c87_min
 +
#PBS -V
 +
cd $PBS_O_WORKDIR
 +
mpirun -np 112 dock6.mpi -i min.in -o 1c87_vir-screen_cartesian-min.out
 +
 +
 +
===Rescoring Docked Molecules===
 +
 +
Now, we are going to rescore virtual screen using foorprint similarity, pharmacophore score, tantimoto score and hungarian and the volume overlap score.
 +
touch descriptor.in
 +
  ../../../zzz.programs/dock6/bin/dock6 -i descriptor.in
 +
 +
Here is the descriptor.in script:
 +
 +
conformer_search_type                                        rigid
 +
use_internal_energy                                          yes
 +
internal_energy_rep_exp                                      12
 +
internal_energy_cutoff                                      100.0
 +
ligand_atom_file                                            /gpfs/projects/AMS536/2018/seung/006.virtual-screen.mpi/1c87.flex_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_score                              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_score_use_footprint_reference_mol2            yes
 +
descriptor_fps_score_footprint_reference_mol2_filename      /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
 +
descriptor_fps_score_foot_compare_type                      Euclidean
 +
descriptor_fps_score_normalize_foot                          no
 +
descriptor_fps_score_foot_comp_all_residue                  yes
 +
descriptor_fps_score_receptor_filename                      /gpfs/projects/AMS536/2018/seung/001.files/rec_withH_1c87.mol2
 +
descriptor_fps_score_vdw_att_exp                            6
 +
descriptor_fps_score_vdw_rep_exp                            12
 +
descriptor_fps_score_vdw_rep_rad_scale                      1
 +
descriptor_fps_score_use_distance_dependent_dielectric      yes
 +
descriptor_fps_score_dielectric                              4.0
 +
descriptor_fps_score_vdw_fp_scale                            1
 +
descriptor_fps_score_es_fp_scale                            1
 +
descriptor_fps_score_hb_fp_scale                            0
 +
descriptor_fms_score_use_ref_mol2                            yes
 +
descriptor_fms_score_ref_mol2_filename                      /gpfs/projects/AMS536/2018/seung/005.dock/1c87.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                          /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
 +
descriptor_hms_score_ref_filename                            /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
 +
descriptor_hms_score_matching_coeff                          -5
 +
descriptor_hms_score_rmsd_coeff                              1
 +
descriptor_volume_score_reference_mol2_filename              /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
 +
descriptor_volume_score_overlap_compute_method              analytical
 +
descriptor_weight_fps_score                                  1
 +
descriptor_weight_pharmacophore_score                        1
 +
descriptor_weight_fingerprint_tanimoto                      -1
 +
descriptor_weight_hms_score                                  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.defn
 +
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.out
 +
write_footprints                                            yes
 +
write_hbonds                                                yes
 +
write_orientations                                          no
 +
num_scored_conformers                                        1
 +
rank_ligands                                                no
 +
 +
Now, copy the following into the file:
 +
 +
#! /bin/bash
 +
#PBS -l walltime=48:00:00
 +
#PBS -l nodes=1:ppn=24
 +
#PBS -q long
 +
#PBS -N 1c87_rescored
 +
#PBS -V
 +
cd $PBS_O_WORKDIR
 +
dock6 -i descriptor.in -o rescored.out
 +
 +
 +
We are going to submit the job again.
 +
qsub submit.sh
 +
 +
.

Latest revision as of 14:12, 12 November 2018

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 the AMS 536 class of 2018, using DOCK v6.8 and it shows how to dock a ligand into a receptor.

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.

1C87

In this tutorial we will use PDB code 1C87, which is the crystal structure of protein tyrosine phosphatase 1B complexed with 2-(oxalyl-amino-4,7-dihydro-5H-thieno[2,3-C]pyran-3-carboxylic acid.

Organizing Directories

We are going to create and organize directories so it would be easier for us to find or identify files in each directory. The structure of your directories should look something like this: 

~/gpfs/projects/AMS536/2018/
                          /001.files/
                          /002.surface-spheres/
                          /003.box/
                          /004.grid/
                          /005.dock/
                          /006.virtual-screen/


To locate which directory you are currently in, type pwd. From your personal directory (...../AMS536/2018/yourName/), you can use the mkdir command to make directories. One way to make a number of directories is: 


 mkdir 001.files 002.spheres 003.box 004.grid 005.dock 006.virtual-screen

To see if you were successful, type 

 ls -l

to list directory contents. If you make a mistake, you can change the name of the directory or remove it and try again

to rename a directory.. 

 mv oldname/ newname/

to remove a directory (must be empty) 

 rmdir mydirectory/

II. Preparing the Receptor and Ligand

Download the PDB File (1C87)

We are going to the PDB website (https://www.rcsb.org/) to download 1C87.pdb file and transfer this pdb file to your directory. First, open Chimera and load 1C87.pdb file. For dock, it is important that the receptor is separate entity and saved as its own structure file. We must also think about what parts of the structure should be included in the receptor and what parts should be ignored. (We may want some crystallographic water which interact with the binding site/ligand but may choose to remove other surrounding waters further away. Remove the ligand and save the receptor in mol2.format. "mol2" format shows the bond types (whether it is single or double bond), and the connectivity of the atoms (what is bonded to what). After saving the receptor with no Hydrogen (noh_receptor_1c87.mol2) and the receptor with built Hydrogens and partial charges (rec_withH_charge_1c87.mol2), load the original 1c87.pdb again. Remove the receptor and water, and save the ligand with no Hydrogen (noh_lig_1c87.mol2) and ligand with Hydrogen and partial charges (lig_withH_charge_1c87.mol2). In order to build Hydrogens and add charge on these structures, go to Tools and Structure Editing and use Dockprep. Along with building Hydrogens and adding partial charges, dockprep will try to fill in missing info from the pdb and also save the structures in a mol2 format at the end . (NOTE: For this system, when you add H on ligand, make the charge -1. In other cases, check the article that is related to pdb file to decide how it should be pronated or depronated.)

Go to 

 Tools -> Structure Editing -> Dockprep

Check the appropriate boxes:

(### Need photo of dockprep setup ###)

Or more explicitly: 

Tools -> Structure Editing -> Add H
Tools -> Structure Editing-> Add Charge
Select AMBER ff99SB for Standard Residues, and add a charge of -1.
File -> Save Mol2... -> lig_withH_charge.mol2

Lastly, open the receptor file and tell chimera to show the surface of your receptor. After this is calculated and you can see the surface of the protein, go to Structure Editing and write DMS file. Make sure to save all files. So far, we will have these files ready.

<pre> 
 1c87.pdb             lig_withH_charge_1c87.mol2     noh_receptor_1c87.mol2
 noh_lig_1c87.mol2           rec_withH_1c87.mol2     rec_1c87.dms

After saving these two files with H, transfer files into 001.files directory: 

 scp -r ./*1c87* username@login.seawulf.stonybrook.edu:/gpfs/projects/AMS536/2018/your_directory_name/001.files

III. Generating Receptor Surface and Spheres

Generating the Receptor Surface

Make sure sphere directory is created and open the directory: 

mkdir 002.surface-spheres
cd 002.surface-spheres

Open Chimera and load the receptor surface:

Load the receptor (noh_receptor_1c87.mol2) file, which is located in the 001.files directory. Go to Action -> Surface -> Show and it will show the surface receptor. Next, save a DMS file as noh_surface_rec.dms by clicking Tools -> Structure editing -> Write DMS.

Creating Spheres We are going to create spheres by using the Sphgen program, which is a sphere generation program in Dock.

1. Create an input file 

vi INSPH

Copy this following script in INSPH input file. 

./noh_surface_rec.dms
R
X
0.0
4.0
1.4
1c87.spheres.sph

The first line ./noh_surface_rec.dms specifies the input file. R means that spheres generated will be outside of the receptor surface. X specifies all the points will be used. 0.0 is the distance in angstrom and it will avoid steric clashes. 4.0 is the maximum surface radius of the spheres and 1.4 is the minimum radius in angstroms.The last line 1c87.spheres.sph creates the file that has clustered spheres.

2. Run the Sphgen program 

sphgen -i INSPH -o OUTSPH

-i chooses input file and -o writes toutput file.

3. Open Chimera to visualize the generated spheres Load 1c87.recep.sph file.

The generated 1C87 surface spheres with a receptor

Now, we are going to select the spheres that is close to the native ligand molecule, which would be within 8.0 angstroms of the ligand. 

sphere_selector 1c87.recep.sph ../001.files/lig_withH_charge_1c87.mol2 8.0

1c87.recep.pdbwill be created and load this file in Chimera.

The generated 1C87 surface spheres with a receptor

IV. Generating Box and Grid

Creating the box

Enter 003.box-grid directory. Create input showbox.in: 

Y
8.0
../002.spheres/selected_spheres.sph
1
1c87.box.pdb

The first line Y creates a box (grid). 8.0 is the box length in angstroms. 1 means spheres. The last line, 1c87.box.pdb, outputs the box to the file.

Use this input type: 

showbox < showbox.in

Now we can visualize the box in Chimera. Load 1c87.box.pdb.

The selected spheres with box
Receptor w/ spheres in box

Calculating Energy Grid

Now, we will compute the energy grid. Create grid.in file This saves the interaction energy of the receptor at each grid point 

vi grid.in

Save the output in the file grid.out. 

grid -i grid.in > grid.out

Answer the questions: 

 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        9
 distance_dielectric       yes
 dielectric_factor         4
 bump_filter               yes
 bump_overlap              0.75
 receptor_file             ../001.files/rec_withH_charge_1c87.mol2
 box_file                  ../003.box/1c87.box.pdb
 vdw_definition_file       /opt/AMS536/dock6/parameters/vdw_AMBER_parm99.defn
 score_grid_prefix         grid

Once completed to answer the questions, it will create two output files: grid.bmp and grid.nrg

V. Docking a Single Molecule for Pose Reproduction

We would like to re-dock 1C87 into a receptor by using these methods: rigid docking, flex docking and fixed anchor docking.

Minimization

First, we are going to run minimization. Minimization is important because the structure of the ligand in the crystal may not be suited for the forcefield used in our calculations. Minimization relaxes the structure to the forcefield so that contact geometries are at an energy minimum. Create min.in file and run dock with this input file. 

touch min.in
../../../zzz.programs/dock6/bin/dock6 -i min.in

Answer the questions like we did above: 

 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100.0
 ligand_atom_file                                             ../001.files/lig_withH_charge_1c87.mol2
 limit_max_ligands                                            no
 skip_molecule                                                no
 read_mol_solvation                                           no
 calculate_rmsd                                               yes
 use_rmsd_reference_mol                                       yes
 rmsd_reference_filename                                     ../001.files/lig_withH_charge_1c87.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                                           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                                      ../001.files/rec_withH_charge_1c87.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                                        yes
 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   
 ligand_outfile_prefix                                        1c87.min
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no

Now, 1c87.min_scored.mol2 will be produced. We can compare minimized ligand structure with unminimized ligand structure.

Rigid Docking

We are going to run rigid docking. With rigid docking the internal degrees of freedom of the ligand are constrained, and only translation and rotation of the whole structure is considered for docking. In order to run dock, we will have to create rigid.in file. 

touch rigid.in
../../../zzz.programs/dock6/bin/dock6 -i rigid.in

Answer the questions: 

  conformer_search_type                                        rigid
  use_internal_energy                                          yes
  internal_energy_rep_exp                                      12
  internal_energy_cutoff                                       100.0
  ligand_atom_file                                             ../001.files/lig_withH_charge_1c87.mol2
  limit_max_ligands                                            no
  skip_molecule                                                no
  read_mol_solvation                                           no
  calculate_rmsd                                               yes
  use_rmsd_reference_mol                                       yes
  rmsd_reference_filename                                   ../001.files/lig_withH_charge_1c87.mol2
  use_database_filter                                          no
  orient_ligand                                                yes
  automated_matching                                           yes
  receptor_site_file                                           ../002.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                                       ../004.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                                                /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                                        1c87.rigid_restraint
  write_orientations                                           no
  num_scored_conformers                                        40
  write_conformations                                          no
  cluster_conformations                                        yes
  cluster_rmsd_threshold                                       2.0
  rank_ligands                                                 no

This will create 1c87.rigid_restraint_scored.mol2 in 005.dock directory.

We can view rigid docking output file in CHIMERA. 

Tools --> Surface/Binding Analysis -> ViewDock
Select the output file -> select Dock 4,5, or 6
Grid Score of 1C87


We can view the lowest grid score (the best ligand that fits in the pocket) , which is what we want, by clicking Column -> Show -> Grid_Score.

Grid Score of 1C87

Flexible Docking

In order to run flexible docking, we will also create flex.in file. In this case, ligand and internal angles are flexible (NOTE: orient_ligand is set to yes). 

touch flex.in
 ../../../zzz.programs/dock6/bin/dock6 -i flex.in

When prompted, answer the questions as follows 

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                                            ../001.files/lig_withH_charge_1c87.mol2
limit_max_ligands                                            no
skip_molecule                                                no
read_mol_solvation                                           no
calculate_rmsd                                               yes
use_rmsd_reference_mol                                       yes
rmsd_reference_filename                                     ../001.files/lig_withH_charge_1c87.mol2
use_database_filter                                          no
orient_ligand                                                yes
automated_matching                                           yes
receptor_site_file                                           ../002.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                                        ../004.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                                        yes
simplex_coefficient_restraint                                10.0
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                                        1c87_flex
write_orientations                                           no
num_scored_conformers                                        40
write_conformations                                          no
cluster_conformations                                        yes
cluster_rmsd_threshold                                       2.0
rank_ligands                                                 no

The output file is 3pgl.flex_scored.mol2, and is saved in the directory from where you ran dock. Then we can visualize the flexible docking in Chimera:

Flexible docking

Fixed Anchor Docking

We are going to run fixed anchor docking. Create fad.in.The central ligand structure is fixed in the active site and internal angles are flexible and sampled. (NOTE: orient_ligand is set to no). 

touch fad.in
 ../../../zzz.programs/dock6/bin/dock6 -i fad.in

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                                             ../001.files/lig_withH_charge_1c87.mol2
limit_max_ligands                                            no
skip_molecule                                                no
read_mol_solvation                                           no
calculate_rmsd                                               yes
use_rmsd_reference_mol                                       yes
rmsd_reference_filename                                      ../001.files/lig_withH_charge_1c87.mol2
use_database_filter                                          no
orient_ligand                                                yes
automated_matching                                           yes
receptor_site_file                                           ../002.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                                       ../004.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                                        yes
simplex_coefficient_restraint                                10.0
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                                        1c87_fad
write_orientations                                           no
num_scored_conformers                                        40
write_conformations                                          no
cluster_conformations                                        yes
cluster_rmsd_threshold                                       2.0
rank_ligands                                                 no

Again, we can view fixed anchor docking in Chimera and find the lowest grid score value.

Anchor fixed docking

Now we have created output files for rigid, flexible and fixed anchor docking.

Molecular Footprint

When it comes to examining how a ligand interacts with its receptor, we can break it down to the interaction contribution of side chains in the binding pocket. Some interactions will be much stronger than other, with some being strongly favorable or strongly prohibitive to binding. Typically we want to look at those residues which contribute most to the interaction between ligand and receptor, favorable or not, so that we can deduce ideal physiochemical structures for future lead refinement. Dock can show this in what is called a Footprint. The plot shows electrostatic and van der Waals interaction contribution for the top X highest residue interactions vs the contribution of another ligand (In this case, we want to see the footprint of a given ligand with respect to the known binder ligand 1c87 from either its crystal structure, or its minimized structure). Here we will compare the footprint of minimized 1c87 vs its crystal structure.

Create a new directory for the footprint calculations, and create the input file for dock inside that directory: 

 mkdir 007.footprint
 cd !$
 touch footprint.in

Answer the questions to fill in the input parameters like before: 

 conformer_search_type                                        rigid
 use_internal_energy                                          no
 ligand_atom_file                                             ../001.files/1c87.min.mol2 ### NOTE: This can be changes based on what you want to compare with what ###
 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                        ../001.files/lig_withH_charge_1c87.mol2 ### Reference ligand for comparison of footprints ###
 fps_foot_compare_type                                        Euclidean
 fps_normalize_foot                                           no
 fps_foot_comp_all_residue                                    yes
 fps_receptor_filename                                        ../001.files/rec_withH_1c87.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                                                /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                                        fps.min.output
 write_footprints                                             yes
 write_hbonds                                                 yes
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no

This script produces a *.txt file to which you can plot the footprint using the provided python script: plot_footprint_single_magnitude.py

to run this enter the following into your command line: 

 python plot_footprint_single_magnitude.py fps.min.output_footprint_scored.txt 50 ### NOTE: Here, we want the 50 highest contribution interactions. The order of these arguments matter! ###

The result should look something like this:

Footprint ex.png

VI. Virtual Screening

We will create a virtual-screen.in file in 006.virtual-screen 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                                      ../001.files/lig_withH_charge_1c87.mol2
  use_database_filter                                          no
  orient_ligand                                                yes
  automated_matching                                           yes
  receptor_site_file                                           ../002.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                                       ../004.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                                        1c87_virtualscreen
  write_orientations                                           no
  num_scored_conformers                                        1
  write_conformations                                          no
  cluster_conformations                                        yes
  cluster_rmsd_threshold                                       2.0
  rank_ligands                                                 no


To submit the job for docking, create submit.sh file. 

vi submit.sh

Write down the following script in submit.sh file. 

#! /bin/bash
#PBS -l walltime=48:00:00
#PBS -l nodes=1:ppn=24
#PBS -q long
#PBS -N 1c87_rigid 
#PBS -V
cd $PBS_O_WORKDIR 
dock6 -i rigid.in -o 1c87_rigid_docking.out

After creating this script, submit the job. 

qsub submit.sh

We can view whether the job runs or does not run (if there is an error): 

qstat -u username

Large Library Virtual Screening (with MPI)

Using 24 processors for two days will not be sufficient to dock all ~25,000 compounds. We can utilize more processors using the Message Passing Interface (MPI) protocol. (This changes the contents of the submission script a little bit) 

#! /bin/bash
#PBS -l walltime=48:00:00
#PBS -l nodes=4:ppn=28
#PBS -q long
#PBS -N 1c87_rigid 
#PBS -V
cd $PBS_O_WORKDIR 
mpirun -np 112 dock6.mpi -i rigid.in -o 1c87_rigid_docking.out

make sure you are requesting the same number of processors on both cwulf and the MPI command (the -np flag). With these resources, we should finish all 25,000 compounds within one day

Cartesian Minimization of Docked Molecules

Here, cartesian just means we will be using an explicit receptor environment as opposed to the energy grid calculated earlier. These calculations should be calculated in a different sub-directory so first we can create a directory called 007.cartesian-minimization, and create a new min.in file inside that directory. 

 mkdir 008.cartesian-minimization
 cd !$
 touch min.in

Now we can call dock and answer the questions like before: 

 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100.0
 ligand_atom_file                                             1c87_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                                      ../001.files/rec_withH_1c87.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                                        1c87.virtualscreen1.minimized
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no


Since this is another job which needs lots of resources, we submit this to the computing nodes with the submission script. First, create submit.sh. 

vi submit.sh

Now, copy the following into the file: 

#! /bin/bash
#PBS -l walltime=48:00:00
#PBS -l nodes=4:ppn=28
#PBS -q long
#PBS -N 1c87_min 
#PBS -V
cd $PBS_O_WORKDIR 
mpirun -np 112 dock6.mpi -i min.in -o 1c87_vir-screen_cartesian-min.out


Rescoring Docked Molecules

Now, we are going to rescore virtual screen using foorprint similarity, pharmacophore score, tantimoto score and hungarian and the volume overlap score. 

touch descriptor.in
 ../../../zzz.programs/dock6/bin/dock6 -i descriptor.in

Here is the descriptor.in script: 

conformer_search_type                                        rigid 
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
internal_energy_cutoff                                       100.0
ligand_atom_file                                             /gpfs/projects/AMS536/2018/seung/006.virtual-screen.mpi/1c87.flex_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_score                              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_score_use_footprint_reference_mol2            yes
descriptor_fps_score_footprint_reference_mol2_filename       /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
descriptor_fps_score_foot_compare_type                       Euclidean
descriptor_fps_score_normalize_foot                          no
descriptor_fps_score_foot_comp_all_residue                   yes
descriptor_fps_score_receptor_filename                       /gpfs/projects/AMS536/2018/seung/001.files/rec_withH_1c87.mol2
descriptor_fps_score_vdw_att_exp                             6
descriptor_fps_score_vdw_rep_exp                             12
descriptor_fps_score_vdw_rep_rad_scale                       1
descriptor_fps_score_use_distance_dependent_dielectric       yes
descriptor_fps_score_dielectric                              4.0
descriptor_fps_score_vdw_fp_scale                            1
descriptor_fps_score_es_fp_scale                             1 
descriptor_fps_score_hb_fp_scale                             0
descriptor_fms_score_use_ref_mol2                            yes
descriptor_fms_score_ref_mol2_filename                       /gpfs/projects/AMS536/2018/seung/005.dock/1c87.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                          /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
descriptor_hms_score_ref_filename                            /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
descriptor_hms_score_matching_coeff                          -5
descriptor_hms_score_rmsd_coeff                              1
descriptor_volume_score_reference_mol2_filename              /gpfs/projects/AMS536/2018/seung/005.dock/1c87.min_scored.mol2
descriptor_volume_score_overlap_compute_method               analytical
descriptor_weight_fps_score                                  1
descriptor_weight_pharmacophore_score                        1
descriptor_weight_fingerprint_tanimoto                       -1
descriptor_weight_hms_score                                  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.defn 
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.out
write_footprints                                             yes
write_hbonds                                                 yes
write_orientations                                           no
num_scored_conformers                                        1
rank_ligands                                                 no

Now, copy the following into the file: 

#! /bin/bash 
#PBS -l walltime=48:00:00 
#PBS -l nodes=1:ppn=24
#PBS -q long
#PBS -N 1c87_rescored
#PBS -V
cd $PBS_O_WORKDIR
dock6 -i descriptor.in -o rescored.out


We are going to submit the job again. 

qsub submit.sh
.
Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2018_DOCK_tutorial_2_with_PDBID_1C87&oldid=11120"