Difference between revisions of "2022 DOCK tutorial 3 with PDBID 1X70"

From Rizzo_Lab
Jump to: navigation, search
(Docking & Virtual Screening)
("Rescore")
 
(26 intermediate revisions by 2 users not shown)
Line 3: Line 3:
  
 
===DOCK===
 
===DOCK===
 +
DOCK is a suite of software intended for molecular docking experiments. It was originally developed by Robert Sheridan, Renee DesJarlais, and Irwin Kuntz. It has been further developed by many additional authors and has many different versions and features. More information about DOCK or the features of DOCK6.9, the version used in this tutorial are available in the DOCK6 manual[https://dock.compbio.ucsf.edu/DOCK_6/dock6_manual.htm#History]
  
 +
===1X70===
 +
1X70 is a protease responsible for degrading a hormone peptide that has potential in mitigation the impact of type 2 diabetes[https://pubs.acs.org/doi/10.1021/jm0493156]. It is a 174.94 kDa homodimer with 728 residues per monomer with ligand bound in both active sites[https://www.rcsb.org/structure/1x70]. It has a resolution of 2.1 Å with an R value of 0.193 and an R free value of 0.228[https://www.rcsb.org/structure/1x70].
  
 
==System Preparation==
 
==System Preparation==
 +
'''Directory Setup'''
 +
Start by running the following command on two computers:
 +
  a computer which you can access the molecular visualization software Chimera
 +
  a computer capable of running DOCK calculations:
 +
 +
  mkdir 001.structure 002.spheres 003.gridbox 004.energy_min 005.footprint 006.vs 007.cart_min 008.rescore
 +
 +
  IMPORTANT:
 +
  Throughout this tutorial the paths of files and directories may vary.
 +
  It is important to change them to match what you have on your local machine when preparing input files or your calculations will fail!
 +
  Also consider path changes when shuttling files between computers for visualization using Chimera!
 +
 
'''Fetching 1X70'''
 
'''Fetching 1X70'''
 +
 +
Go into the first directory
 +
 +
  cd 001.structure
  
 
Open Chimera and do the following to grab the protein:
 
Open Chimera and do the following to grab the protein:
Line 40: Line 59:
 
Now that we the receptor by itself, we have to clean up the rest of the receptor by adding any missing side chains, dealing with multiple occupancies and mutated residues, and protonating and calculating partial charges.
 
Now that we the receptor by itself, we have to clean up the rest of the receptor by adding any missing side chains, dealing with multiple occupancies and mutated residues, and protonating and calculating partial charges.
  
For 1X70, there are several residues with multiple occupancies.
+
For 1X70 one of the problems you will run into is that there are several residues with multiple occupancies.
  
 
One way to check for these residues is to grep the number of alpha carbons from the pdb in the command line.  
 
One way to check for these residues is to grep the number of alpha carbons from the pdb in the command line.  
Line 47: Line 66:
 
[[File:1X70_multiple_occupancy.png|thumb|center|800px|Two potential occupancies shown for CYS649]]
 
[[File:1X70_multiple_occupancy.png|thumb|center|800px|Two potential occupancies shown for CYS649]]
  
To clean this up and protonate/charge the receptor:
+
To clean this up, avoid other potential issues including mutated residues and missing atoms/chains, and protonate/charge the receptor:
 
   Tools > Structure Editing > Dock Prep
 
   Tools > Structure Editing > Dock Prep
  
Line 80: Line 99:
  
 
'''Surface Generation'''
 
'''Surface Generation'''
 +
 +
Go into the second Directory
 +
 +
  cd 002.spheres
  
 
   Load 1X70 w/o Hydrogens in Chimera
 
   Load 1X70 w/o Hydrogens in Chimera
Line 127: Line 150:
 
===Making the box===
 
===Making the box===
  
In order to make the grid we first have to determine how big our grid will be. To do this:
+
In order to make the grid we first have to determine what residues to include when making the grid. To do this:
  
 
   cd 03.grid
 
   cd 03.grid
Line 161: Line 184:
 
   bump_filter                              yes
 
   bump_filter                              yes
 
   bump_overlap                              0.75
 
   bump_overlap                              0.75
   receptor_file                            ../01.structures/1X70_rec.mol2
+
   receptor_file                            ../01.structures/1X70_rec.mol2 (this is the protonated and charged ligand)
 
   box_file                                  ./1X70.box.pdb
 
   box_file                                  ./1X70.box.pdb
 
   vdw_definition_file                      /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/vdw_AMBER_parm99.defn
 
   vdw_definition_file                      /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/vdw_AMBER_parm99.defn
Line 171: Line 194:
 
   grid -i grid.in -o gridinfo.out
 
   grid -i grid.in -o gridinfo.out
  
Once it is done, vi into gridinfo.out and make sure the charges are all integer and that the calculation finished. If they are not, that likely means there is something wrong with your 1X70_H.mol2, and you may want to go back to the dockprep step and remake this file.
+
Once it is done, vi into gridinfo.out and make sure the charges are all integer and that the calculation finished. If they are not, that likely means there is something wrong with your 1X70_H.mol2, and you may want to go back and double check your receptor prep step and remake this file.
  
 
You should get two files:
 
You should get two files:
 
   grid.nrg (energy grid)
 
   grid.nrg (energy grid)
 
   grid.bmp (bump grid)
 
   grid.bmp (bump grid)
 +
 +
The energy grid is used in the dock calculations to give you energies of interaction between the ligand and protein.
 +
 +
The bump grid is used to make sure the ligand doesn't overlap with any receptor atoms and also contains information of the other grids.
  
 
Try loading them into chimera over your protein.
 
Try loading them into chimera over your protein.
Line 191: Line 218:
 
   internal_energy_rep_exp                                      12
 
   internal_energy_rep_exp                                      12
 
   internal_energy_cutoff                                      100.0
 
   internal_energy_cutoff                                      100.0
   ligand_atom_file                                            ./../01.structures/715.mol2
+
   ligand_atom_file                                            ./../01.structures/715.mol2 (this should be the charged, protonated ligand)
 
   limit_max_ligands                                            no
 
   limit_max_ligands                                            no
 
   skip_molecule                                                no
 
   skip_molecule                                                no
Line 247: Line 274:
  
 
[[File:1X70_lig_minimized.png|thumb|center|800px|The minimized ligand (tan) overlayed on the unminimized ligand (blue)]]
 
[[File:1X70_lig_minimized.png|thumb|center|800px|The minimized ligand (tan) overlayed on the unminimized ligand (blue)]]
 +
 +
==Footprint Generation==
 +
 +
Now if we wanted to see any changes between the interactions of the unminimized ligand and the minized pose, we can create the following file "footprint.in":
 +
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          no
 +
  ligand_atom_file                                            ./../04.energy_min/1X70.lig.min_scored.mol2
 +
  limit_max_ligands                                            no
 +
  skip_molecule                                                no
 +
  read_mol_solvation                                          no
 +
  calculate_rmsd                                              no
 +
  use_database_filter                                          no
 +
  orient_ligand                                                no
 +
  bump_filter                                                  no
 +
  score_molecules                                              yes
 +
  contact_score_primary                                        no
 +
  contact_score_secondary                                      no
 +
  grid_score_primary                                          no
 +
  grid_score_secondary                                        no
 +
  multigrid_score_primary                                      no
 +
  multigrid_score_secondary                                    no
 +
  dock3.5_score_primary                                        no
 +
  dock3.5_score_secondary                                      no
 +
  continuous_score_primary                                    no
 +
  continuous_score_secondary                                  no
 +
  footprint_similarity_score_primary                          yes
 +
  footprint_similarity_score_secondary                        no
 +
  fps_score_use_footprint_reference_mol2                      yes
 +
  fps_score_footprint_reference_mol2_filename                  ./../01.structures/715.mol2
 +
  fps_score_foot_compare_type                                  Euclidean
 +
  fps_score_normalize_foot                                    no
 +
  fps_score_foot_comp_all_residue                              yes
 +
  fps_score_receptor_filename                                  ./../01.structures/1X70_rec.mol2
 +
  fps_score_vdw_att_exp                                        6
 +
  fps_score_vdw_rep_exp                                        9
 +
  fps_score_vdw_rep_rad_scale                                  1
 +
  fps_score_use_distance_dependent_dielectric                  yes
 +
  fps_score_dielectric                                        4.0
 +
  fps_score_vdw_fp_scale                                      1
 +
  fps_score_es_fp_scale                                        1
 +
  fps_score_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.9_release/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 +
  flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        footprint.out
 +
  write_footprints                                            yes
 +
  write_hbonds                                                yes
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 +
 +
Then we run it using dock:
 +
  dock6 -i footprint.in -o footprint.out
 +
 +
We can then plot the 50 most significant residues using the following python script:
 +
  python plot_footprint_single_magnitude.py footprint.out_footprint_scored.txt 50
 +
 +
[[File:1X70_min_fp.png|thumb|center|800px|The energies of the minimized and unminimized ligand]]
  
 
=Docking & Virtual Screening=
 
=Docking & Virtual Screening=
Line 501: Line 596:
 
   #SBATCH --output=1x70_flex.txt
 
   #SBATCH --output=1x70_flex.txt
 
   #SBATCH -p short-24core
 
   #SBATCH -p short-24core
 
+
 
 
   cd $SLURM_SUBMIT_DIR
 
   cd $SLURM_SUBMIT_DIR
 
+
 
 
   dock6 -i flex.in -o 1x70_flex.out
 
   dock6 -i flex.in -o 1x70_flex.out
  
Line 516: Line 611:
 
If the entry under the st (status) says pd, then that means your code is currently waiting to start.  If it says r, then it is currently being run and it will show the elapsed time as well.  Once the job has finished, all of the output files should be left in your current directory.  You can load the .mol2 file in Chimera to see if it matches your crystal structure.
 
If the entry under the st (status) says pd, then that means your code is currently waiting to start.  If it says r, then it is currently being run and it will show the elapsed time as well.  Once the job has finished, all of the output files should be left in your current directory.  You can load the .mol2 file in Chimera to see if it matches your crystal structure.
  
[[File:.png|thumb|center|800px| placeholder text]]
+
= Virtual Screen =
 +
Opposed to looking at a single molecule, it may be desired to look at an entire library of molecules through a virtual screen. Fortunately, DOCK will read multiple separate mol2s from the same input file.
 +
 
 +
  cd 07.vs
 +
  vi vs.in
 +
 
 +
    search_type                                                  flex
 +
    write_fragment_libraries                                    no
 +
    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                                      9
 +
    internal_energy_cutoff                                      100.0
 +
    ligand_atom_file                                            /gpfs/projects/AMS536/zzz.programs/VS_library_25K.mol2 #this is where the name of your library goes
 +
    limit_max_ligands                                            no
 +
    skip_molecule                                                no
 +
    read_mol_solvation                                          no
 +
    calculate_rmsd                                              no
 +
    use_database_filter                                          no
 +
    orient_ligand                                                yes
 +
    automated_matching                                          yes
 +
    receptor_site_file                                          ../02.spheres/selected_spheres.sph
 +
    max_orientations                                            1000
 +
    critical_points                                              no
 +
    chemical_matching                                            no
 +
    use_ligand_spheres                                          no
 +
    bump_filter                                                  no
 +
    score_molecules                                              yes
 +
    contact_score_primary                                        no
 +
    contact_score_secondary                                      no
 +
    grid_score_primary                                          yes
 +
    grid_score_secondary                                        no
 +
    grid_score_rep_rad_scale                                    1
 +
    grid_score_vdw_scale                                        1
 +
    grid_score_es_scale                                          1
 +
    grid_score_grid_prefix                                      ../03.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.9_release/parameters/vdw_AMBER_parm99.defn
 +
    flex_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 +
    flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 +
    ligand_outfile_prefix                                        virtual_25k.out
 +
    write_orientations                                          no
 +
    num_scored_conformers                                        1
 +
    rank_ligands                                                no
 +
 
 +
Running this using a single node will take a long time, to run it on multiple submit using the following script.
 +
 
 +
  #!/bin/sh
 +
  #SBATCH --partition=long-24core
 +
  #SBATCH --time=7-00:00:00
 +
  #SBATCH --nodes=4
 +
  #SBATCH --ntasks=96
 +
  #SBATCH --job-name=1X70_25k_vs
 +
  #SBATCH --output=%x-%j.o
 +
 
 +
  echo "starting DOCK6.9 Simulation"
 +
  module load intel/mpi/64/2018/18.0.3
 +
 
 +
  mpirun -np 32 /gpfs/projects/AMS536/zzz.programs/dock6.9_release/bin/dock6.mpi -i virtual.in -o virtual.out
 +
  echo "done"
 +
 
 +
You can then visualize the results using chimera and the dock viewing tool.
 +
 
 +
== Cartesian Minimization ==
 +
 
 +
We can't look at the results for the virtual screen quite yet.  First, we have to minimize the energy for all of the screened ligands.  To start, make the following directory.
  
= Placeholder =
+
  mkdir 008.cartesianmin
 +
  cd ./008.cartesianmin
  
[[File:|thumb|center|800px|image placeholder]]
+
Then make the input file
  
=Placeholder=
+
  vim min_scored.in
  
=Placeholder=
+
And use the following parameters
 +
 
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100
 +
  ligand_atom_file                                            ../006.virtual_screen_mpi/virtual.out_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
 +
  vcontact_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.structure/1x70_receptor_wH.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.0
 +
  cont_score_es_scale                                          1.0
 +
  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.0
 +
  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.9_release/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 +
  flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 +
  ligand_outfile_prefix                                        1x70.virtual_screen.min
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 
 +
Since this will take a while to run, we need to submit a slurm script to run this job.  We can make the following file.
 +
 
 +
  vim min.sh
 +
 
 +
And add the following options.
 +
 
 +
  #!/bin/bash
 +
  #SBATCH --time=48:00:00
 +
  #SBATCH --nodes=2
 +
  #SBATCH --ntasks=28
 +
  #SBATCH --job-name=mpi_1HW9_VS_AR
 +
  #SBATCH --output=1HW9_VS_mpi.out
 +
  #SBATCH -p long-28core
 +
 
 +
  cd $SLURM_SUBMIT_DIR
 +
  mpirun -np 56 /gpfs/projects/AMS536/zzz.programs/dock6.9_release/bin/dock6.mpi -i min_scored.in -o min_scored.out
 +
 
 +
== Rescore ==
 +
 
 +
This is the final step of the virtual screen.  After this we'll be able to look at various parameters and rank order the molecules in the virtual screen based on each one.  To start make the a new directory and cd into it.
 +
 
 +
  mkdir 009.rescore
 +
  cd ./009.rescore
 +
 
 +
And we can make the following input file.
 +
 
 +
  vim rescore.in
 +
 
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100
 +
  ligand_atom_file                                            ../006.virtual_screen_mpi/virtual.out_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_hungarian                                    yes
 +
  descriptor_use_volume_overlap                                yes
 +
  descriptor_fps_score_use_footprint_reference_mol2            yes
 +
  descriptor_fps_score_footprint_reference_mol2_filename      ../004.dock/1x70.ligand.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                      ../001.structure/1x70_receptor_wH.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                      ../004.dock/1HW9.lig.min_scored.mol2
 +
  descriptor_fms_score_write_reference_pharmacophore_mol2      no
 +
  descriptor_fms_score_write_reference_pharmacophore_txt      no
 +
  descriptor_fms_score_write_candidate_pharmacophore          no
 +
  descriptor_fms_score_write_matched_pharmacophore            no
 +
  descriptor_fms_score_compare_type                            overlap
 +
  descriptor_fms_score_full_match                              yes
 +
  descriptor_fms_score_match_rate_weight                      5.0
 +
  descriptor_fms_score_match_dist_cutoff                      1.0
 +
  descriptor_fms_score_match_proj_cutoff                      .7071
 +
  descriptor_fms_score_max_score                              20
 +
  descriptor_fingerprint_ref_filename                          ../004.dock/1x70.ligand.min_scored.mol2
 +
  descriptor_hms_score_ref_filename                            ../004.dock/1x70.ligand.min_scored.mol2
 +
  descriptor_hms_score_matching_coeff                          -5
 +
  descriptor_hms_score_rmsd_coeff                              1
 +
  descriptor_volume_score_reference_mol2_filename              ../004.dock/1HW9.lig.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.9_release/parameters/vdw_AMBER_parm99.defn
 +
  flex_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 +
  flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 +
  chem_defn_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/chem.defn
 +
  pharmacophore_defn_file                                      /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/ph4.defn
 +
  ligand_outfile_prefix                                        descriptor.output
 +
  write_footprints                                            yes
 +
  write_hbonds                                                yes
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 
 +
Like before, we should submit this to the cluster using slurm.  This can be done with the following .sh file.
 +
 
 +
  vim rescore.sh
 +
 
 +
  #!/bin/sh
 +
  #SBATCH --partition=long-40core
 +
  #SBATCH --time=48:00:00
 +
  #SBATCH --nodes=4
 +
  #SBATCH --ntasks=32
 +
  #SBATCH --job-name=tutorial_run
 +
  #SBATCH --output=%x-%j.o
 +
 
 +
  module load intel/mpi/64/2018/18.0.3
 +
  mpirun -np 32 dock6.mpi -i rescore.in -o rescore.out
  
 +
After this is finished running.  The finals results can be loaded into Chimera.  It's best to visualize it using ViewDock.  This can be found by
  
[[File:_mpi_minimized.png|thumb|center|800px|the image shows ]]
+
  Tools -> Structure/Binding Analysis -> ViewDock
  
=placeholder=
+
This allows you to sort the ligands by any of the listed parameters and cycle through them to see how they fit in the receptor.

Latest revision as of 14:45, 9 March 2022

Introduction

DOCK

DOCK is a suite of software intended for molecular docking experiments. It was originally developed by Robert Sheridan, Renee DesJarlais, and Irwin Kuntz. It has been further developed by many additional authors and has many different versions and features. More information about DOCK or the features of DOCK6.9, the version used in this tutorial are available in the DOCK6 manual[1]

1X70

1X70 is a protease responsible for degrading a hormone peptide that has potential in mitigation the impact of type 2 diabetes[2]. It is a 174.94 kDa homodimer with 728 residues per monomer with ligand bound in both active sites[3]. It has a resolution of 2.1 Å with an R value of 0.193 and an R free value of 0.228[4].

System Preparation

Directory Setup Start by running the following command on two computers:

 a computer which you can access the molecular visualization software Chimera
 a computer capable of running DOCK calculations:
 mkdir 001.structure 002.spheres 003.gridbox 004.energy_min 005.footprint 006.vs 007.cart_min 008.rescore
 IMPORTANT: 
 Throughout this tutorial the paths of files and directories may vary. 
 It is important to change them to match what you have on your local machine when preparing input files or your calculations will fail! 
 Also consider path changes when shuttling files between computers for visualization using Chimera!

Fetching 1X70

Go into the first directory

 cd 001.structure

Open Chimera and do the following to grab the protein:

 File > Fetch By ID > 1X70
Biological assembly 1 of 1X70.

The first thing to notice is that this is a dimer and the ligand, 715, is not bound at the dimer interface. Thus, one of the monomers is entirely redundant and should be deleted.

 Select > Chain > A
 Actions > Atoms > Delete

Next it is important to remove cofactors, ions, and water molecules not involved in the binding interactions. This can be checked by reading through the paper associated with the PDB.

 Select > Residue > NAG > Actions > Atoms > Delete
 Select > Residue > NDG > Actions > Atoms > Delete
 Select > Residue > HOH > Actions > Atoms > Delete

This will leave us with just the ligand and receptor.

1 monomer of 1X70 with all waters and non standard residues deleted except for the ligand 715, which is colored light blue for clarity.

Receptor Prep

Receptor 1X70 without Hydrogens

Now that we the receptor by itself, we have to clean up the rest of the receptor by adding any missing side chains, dealing with multiple occupancies and mutated residues, and protonating and calculating partial charges.

For 1X70 one of the problems you will run into is that there are several residues with multiple occupancies.

One way to check for these residues is to grep the number of alpha carbons from the pdb in the command line.

 In the terminal, type: grep -e CA 1X70_noH.pdb
Two potential occupancies shown for CYS649

To clean this up, avoid other potential issues including mutated residues and missing atoms/chains, and protonate/charge the receptor:

 Tools > Structure Editing > Dock Prep
Using dock prep to fill in missing atoms/sidechains, add Hydrogens and charges.

Make sure that the protonation makes sense for residues in the active site or coordinated with metals (none here), especially histidines, by checking the paper and chimera for any nitrogens coordinating with metals are not protonated or residues in the active site with differing protonation states. Also make sure that all residues have integral charge (i.e. the charge is an integer). Chimera should give a warning after the dockprep if any residue charges are non-integral.

Ligand Prep

First you have to save the ligand as a separate file. You can do this in Chimera by deleting all of the protein and saving that as a separate file: "715_noH.pdb".

 Select > Residues > 715 > Invert Selection
 Actions > Atoms > delete
 File > save PDB 

Now you should just have the ligand.

715 without Hydrogens added

Next Hydrogens and partial atomic charges need to be added and saved as "715_H.mol2". It will ask for the ligands overall charge, which you should verify using chemical knowledge.

 tools > structure editing > addH 
 tools > structure editing > addCharge > Select the correct charge for your ligand, Use AM1-BCC. 
 File > save Mol2
715 with Hydrogens

Surface Generation & Spheres

This section details the generation of sphere files which will be used to describe where you are trying to DOCK to on your protein.

Surface Generation

Go into the second Directory

 cd 002.spheres
 Load 1X70 w/o Hydrogens in Chimera
 actions > show > surface


The VDW surface for 1X70 that will be used to generate the sphere files.

Then you will write a DMS (Molecular Surface) File With the surface generated in Chimera:

 Tools > Structure Editing > Write DMS

Save this as 1X70_dms.dms, and now you should have a DMS file for the next step.

Sphere Generation To generate spheres make the following input file: "INSPH" -

 1X70_dms.dms      #Molecular Surface File 
 R                 #Whether to generate spheres outside of surface (R) or inside (L) 
 X                 #Surface points from the DMS file to use in sphere generation 
 0                 #Minimum radius between spheres
 4.0               #Maximum radius of sphere
 1.4               #Minimum radius of sphere
 1X70_wo_H.sph     #Output sphere file

For more information on sphere generation see: https://dock.compbio.ucsf.edu/DOCK_6/tutorials/sphere_generation/generating_spheres.htm

The .sph file should give you something similar to the following image if you load it up over your protein in Chimera:

The spheres representing the empty space within the protein.


Sphere Selection

Now we need to select the subset of all those spheres that will be used in the docking experiment by excluding those too far away from the ligand. We can do this through the following command:

 #sphere_selector <receptor> <ligand> <radius in angstroms>
 sphere_selector 1X70_wo_H.sph 715_H.mol2 10.0
 
The spheres we just selected to be using in docking experiments.

Making The Infamous Grid

Making the box

In order to make the grid we first have to determine what residues to include when making the grid. To do this:

 cd 03.grid
 showbox 
   #automatically construct box to enclose spheres [Y/N]
   Y
   #extra margin to also be enclosed (angstroms)?
    #(this will be added in all 6 directions)
   8.0
   #sphere file-
   ./../02.spheres/selected_spheres.sph
   #cluster number-
   1
   #output filename?
   1X70.box.pdb

You can also copy these inputs into "showbox.in" (none of the commented lines) and then type "showbox < showbox.in"

Making the grid

Making the following input file "grid.in":

 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                             ../01.structures/1X70_rec.mol2 (this is the protonated and charged ligand)
 box_file                                  ./1X70.box.pdb
 vdw_definition_file                       /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/vdw_AMBER_parm99.defn
 score_grid_prefix                         grid


Now generate the grid. It should take several minutes.

 grid -i grid.in -o gridinfo.out

Once it is done, vi into gridinfo.out and make sure the charges are all integer and that the calculation finished. If they are not, that likely means there is something wrong with your 1X70_H.mol2, and you may want to go back and double check your receptor prep step and remake this file.

You should get two files:

 grid.nrg (energy grid)
 grid.bmp (bump grid)

The energy grid is used in the dock calculations to give you energies of interaction between the ligand and protein.

The bump grid is used to make sure the ligand doesn't overlap with any receptor atoms and also contains information of the other grids.

Try loading them into chimera over your protein.

The energy grid for 1X70.
The bump grid for 1X70.

Energy Minimization for the ligand

Now cd into directory 04.energy_min and make the following input file "min.in":

 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100.0
 ligand_atom_file                                             ./../01.structures/715.mol2 (this should be the charged, protonated ligand)
 limit_max_ligands                                            no
 skip_molecule                                                no
 read_mol_solvation                                           no
 calculate_rmsd                                               yes
 use_rmsd_reference_mol                                       yes
 rmsd_reference_filename                                      ./../01.structures/715.mol2
 use_database_filter                                          no
 orient_ligand                                                no
 bump_filter                                                  no
 score_molecules                                              yes
 contact_score_primary                                        no
 contact_score_secondary                                      no
 grid_score_primary                                           yes
 grid_score_secondary                                         no
 grid_score_rep_rad_scale                                     1
 grid_score_vdw_scale                                         1
 grid_score_es_scale                                          1
 grid_score_grid_prefix                                       ./../03.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.9_release/parameters/vdw_AMBER_parm99.defn
 flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 ligand_outfile_prefix                                        1X70.lig.min
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no

Then run this through dock

 dock6 -i min.in -o min.out

You should get an energy minimized version of the ligand. Try pulling it up over the original ligand in chimera.

The minimized ligand (tan) overlayed on the unminimized ligand (blue)

Footprint Generation

Now if we wanted to see any changes between the interactions of the unminimized ligand and the minized pose, we can create the following file "footprint.in":

  conformer_search_type                                        rigid
 use_internal_energy                                          no
 ligand_atom_file                                             ./../04.energy_min/1X70.lig.min_scored.mol2
 limit_max_ligands                                            no
 skip_molecule                                                no
 read_mol_solvation                                           no
 calculate_rmsd                                               no
 use_database_filter                                          no
 orient_ligand                                                no
 bump_filter                                                  no
 score_molecules                                              yes
 contact_score_primary                                        no
 contact_score_secondary                                      no
 grid_score_primary                                           no
 grid_score_secondary                                         no
 multigrid_score_primary                                      no
 multigrid_score_secondary                                    no
 dock3.5_score_primary                                        no
 dock3.5_score_secondary                                      no
 continuous_score_primary                                     no
 continuous_score_secondary                                   no
 footprint_similarity_score_primary                           yes
 footprint_similarity_score_secondary                         no
 fps_score_use_footprint_reference_mol2                       yes
 fps_score_footprint_reference_mol2_filename                  ./../01.structures/715.mol2
 fps_score_foot_compare_type                                  Euclidean
 fps_score_normalize_foot                                     no
 fps_score_foot_comp_all_residue                              yes
 fps_score_receptor_filename                                  ./../01.structures/1X70_rec.mol2
 fps_score_vdw_att_exp                                        6
 fps_score_vdw_rep_exp                                        9
 fps_score_vdw_rep_rad_scale                                  1
 fps_score_use_distance_dependent_dielectric                  yes
 fps_score_dielectric                                         4.0
 fps_score_vdw_fp_scale                                       1
 fps_score_es_fp_scale                                        1
 fps_score_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.9_release/parameters/vdw_AMBER_parm99.defn
 flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 ligand_outfile_prefix                                        footprint.out
 write_footprints                                             yes
 write_hbonds                                                 yes
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no


Then we run it using dock:

 dock6 -i footprint.in -o footprint.out

We can then plot the 50 most significant residues using the following python script:

 python plot_footprint_single_magnitude.py footprint.out_footprint_scored.txt 50
The energies of the minimized and unminimized ligand

Docking & Virtual Screening

Now that the ligand and receptor files are set up, we can begin docking the ligand back into the receptor. This makes sure that files are set up correctly, since we should expect the docking results to closely match the crystal structure from the PDB.

Rigid Docking

The first docking method we will be using is rigid docking. With this method, the ligand is kept in a fixed position while the algorithm tries to move/rotate the entire molecule to reduce energy. This is the least computationally intensive method for docking. We can start by creating the input file with the command

 vim rigid.in

And using the following parameters.

 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100
 ligand_atom_file                                             1x70.ligand.min_scored.mol2
 limit_max_ligands                                            no
 skip_molecule                                                no
 read_mol_solvation                                           no
 calculate_rmsd                                               yes
 use_rmsd_reference_mol                                       yes
 rmsd_reference_filename                                      1x70.ligand.min_scored.mol2
 use_database_filter                                          no
 orient_ligand                                                yes
 automated_matching                                           yes
 receptor_site_file                                           ../002.surface_spheres/selected_spheres.sph
 max_orientations                                             2000
 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                                       ../003.gridbox/grid
 multigrid_score_secondary                                    no
 dock3.5_score_secondary                                      no
 continuous_score_secondary                                   no
 footprint_similarity_score_secondary                         no
 pharmacophore_score_secondary                                no
 descriptor_score_secondary                                   no
 gbsa_zou_score_secondary                                     no
 gbsa_hawkins_score_secondary                                 no
 SASA_score_secondary                                         no
 amber_score_secondary                                        no
 minimize_ligand                                              yes
 simplex_max_iterations                                       1000
 simplex_tors_premin_iterations                               0
 simplex_max_cycles                                           1
 simplex_score_converge                                       0.1
 simplex_cycle_converge                                       1.0
 simplex_trans_step                                           1.0
 simplex_rot_step                                             0.1
 simplex_tors_step                                            10.0
 simplex_random_seed                                          0
 simplex_restraint_min                                        no
 atom_model                                                   all
 vdw_defn_file                                                /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/vdw_AMBER_parm99.defn
 flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 ligand_outfile_prefix                                        rigid.out
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no
File:Placeholder.png
placeholder text

Fixed Anchor Docking

In fixed anchor docking, part of the ligand remains fixed, and the rest of it is given some rotational degrees of freedom. This makes it slightly more computationally intensive than rigid docking, but still less intensive than flex docking. So it serves as a middle ground between the two. We can start by creating the input file.

 vim fixed.in

Then adding the following parameters.

 conformer_search_type                                        flex
 write_fragment_libraries                                     no
 user_specified_anchor                                        no
 limit_max_anchors                                            no
 min_anchor_size                                              5
 pruning_use_clustering                                       yes
 pruning_max_orients                                          10000
 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.structure/1x70_ligand_wH.mol2
 limit_max_ligands                                            no
 skip_molecule                                                no
 read_mol_solvation                                           no
 calculate_rmsd                                               yes
 use_rmsd_reference_mol                                       yes
 rmsd_reference_filename                                      ../001.structure/1x70_ligand_wH.mol2
 use_database_filter                                          no
 orient_ligand                                                no
 bump_filter                                                  no
 score_molecules                                              yes
 contact_score_primary                                        no
 contact_score_secondary                                      no
 grid_score_primary                                           yes
 grid_score_secondary                                         no
 grid_score_rep_rad_scale                                     1
 grid_score_vdw_scale                                         1
 grid_score_es_scale                                          1
 grid_score_grid_prefix                                       ../003.gridbox/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
 simplex_trans_step                                           1
 simplex_rot_step                                             0.1
 simplex_tors_step                                            10
 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.9_release/parameters/vdw_AMBER_parm99.defn
 flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 ligand_outfile_prefix                                        1x70_fad
 write_orientations                                           no
 num_scored_conformers                                        100
 write_conformations                                          no
 cluster_conformations                                        yes
 cluster_rmsd_threshold                                       2.0
 rank_ligands                                                 no


File:Place holder.png
there is no figure here

Flexible Docking

In flex docking, the ligand is given full rotational freedom, so the software can put it in whatever conformation it wants. Because of this, flex docking is the most computationally intensive docking method. We can start by creating the input file.

 vim flex.in

Then adding the following parameters

 conformer_search_type                                        flex
 write_fragment_libraries                                     no
 user_specified_anchor                                        no
 limit_max_anchors                                            no
 min_anchor_size                                              5
 pruning_use_clustering                                       yes
 pruning_max_orients                                          5000
 pruning_clustering_cutoff                                    2500
 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                                             1x70.ligand.min_scored.mol2
 limit_max_ligands                                            no
 skip_molecule                                                no
 read_mol_solvation                                           no
 calculate_rmsd                                               yes
 use_rmsd_reference_mol                                       yes
 rmsd_reference_filename                                      1x70.ligand.min_scored.mol2
 use_database_filter                                          no
 orient_ligand                                                yes
 automated_matching                                           yes
 receptor_site_file                                           ../002.surface_spheres/selected_spheres.sph
 max_orientations                                             1000
 critical_points                                              no
 chemical_matching                                            no
 use_ligand_spheres                                           no
 bump_filter                                                  no
 score_molecules                                              yes
 contact_score_primary                                        no
 contact_score_secondary                                      no
 grid_score_primary                                           yes
 grid_score_secondary                                         no
 grid_score_rep_rad_scale                                     1
 grid_score_vdw_scale                                         1
 grid_score_es_scale                                          1
 grid_score_grid_prefix                                       ../003.gridbox/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.9_release/parameters/vdw_AMBER_parm99.defn
 flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 ligand_outfile_prefix                                        flex.out
 write_orientations                                           yes
 num_scored_conformers                                        25
 write_conformations                                          yes
 cluster_conformations                                        yes
 cluster_rmsd_threshold                                       2
 rank_ligands                                                 no

Because this method is so intensive, we do not want to run the code on the login node of seawulf. Instead, we should submit a job to run it on one of the compute nodes. To do this, we can create the following slurm script.

 vim slurm.sh

Then fill it in with the following commands.

 #!/bin/bash
 #SBATCH --time=3:00:00
 #SBATCH --nodes=1
 #SBATCH --ntasks=24
 #SBATCH --job-name=1x70_flex
 #SBATCH --output=1x70_flex.txt
 #SBATCH -p short-24core
 
 cd $SLURM_SUBMIT_DIR
 
 dock6 -i flex.in -o 1x70_flex.out

It is worth it to note that short-24core processor may not be the best one to use, as other people may be using all of the nodes. You can use the command

 sinfo

to see which nodes are currently in use and pick one of the short processors that has some nodes idle. After submitting the slurm job, you can check the status of it using the command

 squeue -u netid

If the entry under the st (status) says pd, then that means your code is currently waiting to start. If it says r, then it is currently being run and it will show the elapsed time as well. Once the job has finished, all of the output files should be left in your current directory. You can load the .mol2 file in Chimera to see if it matches your crystal structure.

Virtual Screen

Opposed to looking at a single molecule, it may be desired to look at an entire library of molecules through a virtual screen. Fortunately, DOCK will read multiple separate mol2s from the same input file.

 cd 07.vs
 vi vs.in
    search_type                                                  flex
    write_fragment_libraries                                     no
    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                                      9
    internal_energy_cutoff                                       100.0
    ligand_atom_file                                             /gpfs/projects/AMS536/zzz.programs/VS_library_25K.mol2 #this is where the name of your library goes
    limit_max_ligands                                            no
    skip_molecule                                                no
    read_mol_solvation                                           no
    calculate_rmsd                                               no
    use_database_filter                                          no
    orient_ligand                                                yes
    automated_matching                                           yes
    receptor_site_file                                           ../02.spheres/selected_spheres.sph
    max_orientations                                             1000
    critical_points                                              no
    chemical_matching                                            no
    use_ligand_spheres                                           no
    bump_filter                                                  no
    score_molecules                                              yes
    contact_score_primary                                        no
    contact_score_secondary                                      no
    grid_score_primary                                           yes
    grid_score_secondary                                         no
    grid_score_rep_rad_scale                                     1
    grid_score_vdw_scale                                         1
    grid_score_es_scale                                          1
    grid_score_grid_prefix                                       ../03.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.9_release/parameters/vdw_AMBER_parm99.defn
    flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
    flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
    ligand_outfile_prefix                                        virtual_25k.out 
    write_orientations                                           no
    num_scored_conformers                                        1
    rank_ligands                                                 no
 

Running this using a single node will take a long time, to run it on multiple submit using the following script.

 #!/bin/sh 
 #SBATCH --partition=long-24core
 #SBATCH --time=7-00:00:00
 #SBATCH --nodes=4
 #SBATCH --ntasks=96
 #SBATCH --job-name=1X70_25k_vs
 #SBATCH --output=%x-%j.o
 
 echo "starting DOCK6.9 Simulation"
 module load intel/mpi/64/2018/18.0.3
 
 mpirun -np 32 /gpfs/projects/AMS536/zzz.programs/dock6.9_release/bin/dock6.mpi -i virtual.in -o virtual.out
 echo "done"

You can then visualize the results using chimera and the dock viewing tool.

Cartesian Minimization

We can't look at the results for the virtual screen quite yet. First, we have to minimize the energy for all of the screened ligands. To start, make the following directory.

 mkdir 008.cartesianmin
 cd ./008.cartesianmin

Then make the input file

 vim min_scored.in

And use the following parameters

 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100
 ligand_atom_file                                             ../006.virtual_screen_mpi/virtual.out_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
 vcontact_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.structure/1x70_receptor_wH.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.0
 cont_score_es_scale                                          1.0
 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.0
 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.9_release/parameters/vdw_AMBER_parm99.defn
 flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 ligand_outfile_prefix                                        1x70.virtual_screen.min
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no

Since this will take a while to run, we need to submit a slurm script to run this job. We can make the following file.

 vim min.sh

And add the following options.

 #!/bin/bash
 #SBATCH --time=48:00:00
 #SBATCH --nodes=2
 #SBATCH --ntasks=28
 #SBATCH --job-name=mpi_1HW9_VS_AR
 #SBATCH --output=1HW9_VS_mpi.out
 #SBATCH -p long-28core
 
 cd $SLURM_SUBMIT_DIR
 mpirun -np 56 /gpfs/projects/AMS536/zzz.programs/dock6.9_release/bin/dock6.mpi -i min_scored.in -o min_scored.out

Rescore

This is the final step of the virtual screen. After this we'll be able to look at various parameters and rank order the molecules in the virtual screen based on each one. To start make the a new directory and cd into it.

 mkdir 009.rescore
 cd ./009.rescore

And we can make the following input file.

 vim rescore.in
 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100
 ligand_atom_file                                             ../006.virtual_screen_mpi/virtual.out_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_hungarian                                     yes
 descriptor_use_volume_overlap                                yes
 descriptor_fps_score_use_footprint_reference_mol2            yes
 descriptor_fps_score_footprint_reference_mol2_filename       ../004.dock/1x70.ligand.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                       ../001.structure/1x70_receptor_wH.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                       ../004.dock/1HW9.lig.min_scored.mol2
 descriptor_fms_score_write_reference_pharmacophore_mol2      no
 descriptor_fms_score_write_reference_pharmacophore_txt       no
 descriptor_fms_score_write_candidate_pharmacophore           no
 descriptor_fms_score_write_matched_pharmacophore             no
 descriptor_fms_score_compare_type                            overlap
 descriptor_fms_score_full_match                              yes
 descriptor_fms_score_match_rate_weight                       5.0
 descriptor_fms_score_match_dist_cutoff                       1.0
 descriptor_fms_score_match_proj_cutoff                       .7071
 descriptor_fms_score_max_score                               20
 descriptor_fingerprint_ref_filename                          ../004.dock/1x70.ligand.min_scored.mol2
 descriptor_hms_score_ref_filename                            ../004.dock/1x70.ligand.min_scored.mol2
 descriptor_hms_score_matching_coeff                          -5
 descriptor_hms_score_rmsd_coeff                              1
 descriptor_volume_score_reference_mol2_filename              ../004.dock/1HW9.lig.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.9_release/parameters/vdw_AMBER_parm99.defn
 flex_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex.defn
 flex_drive_file                                              /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/flex_drive.tbl
 chem_defn_file                                               /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/chem.defn
 pharmacophore_defn_file                                      /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/ph4.defn
 ligand_outfile_prefix                                        descriptor.output
 write_footprints                                             yes
 write_hbonds                                                 yes
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no

Like before, we should submit this to the cluster using slurm. This can be done with the following .sh file.

 vim rescore.sh
 #!/bin/sh
 #SBATCH --partition=long-40core
 #SBATCH --time=48:00:00
 #SBATCH --nodes=4
 #SBATCH --ntasks=32
 #SBATCH --job-name=tutorial_run
 #SBATCH --output=%x-%j.o
 
 module load intel/mpi/64/2018/18.0.3
 mpirun -np 32 dock6.mpi -i rescore.in -o rescore.out

After this is finished running. The finals results can be loaded into Chimera. It's best to visualize it using ViewDock. This can be found by

 Tools -> Structure/Binding Analysis -> ViewDock

This allows you to sort the ligands by any of the listed parameters and cycle through them to see how they fit in the receptor.