Difference between revisions of "2022 DOCK tutorial 1 with PDBID 6ME2"

From Rizzo_Lab
Jump to: navigation, search
(→‎Basics of the Command Line/Terminal)
(→‎Generating DMS of receptor)
 
(90 intermediate revisions by the same user not shown)
Line 1: Line 1:
−
Introduction =  
+
===Required Software===
 +
'''UCSF Chimera, UCSF DOCK 6.9'''
  
 
===Chimera===
 
===Chimera===
−
UCSF Chimera is a software which aids the interactive 3D visualization, editing, and analysis of biological macromolecules.
+
'''UCSF Chimera''' is a software which aids the interactive 3D visualization, editing, and analysis of biological macromolecules.
  
 
[insert screenshot of chimera]
 
[insert screenshot of chimera]
Line 9: Line 10:
  
 
===DOCK===
 
===DOCK===
−
(some info about DOCK)
+
'''UCSF DOCK''' is a molecular docking program, used to find the optimal orientation of a ligand binding to a receptor. DOCK can be used for a number of different molecular modeling/design tasks, such as testing binding affinities, virtual screening, and de novo design.
  
 
===6ME2===
 
===6ME2===
−
(some info about the structure)
+
6ME2 is the modeled structure of receptor MT1 of melatonin complexed with ramelteon (the ligand). Further details on the structure from the original researchers can be found at https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6696938/
−
 
 
−
===Required Software===
 
−
(chimera, dock,...)
 
  
 
= Seawulf: Getting Started =  
 
= Seawulf: Getting Started =  
Line 94: Line 92:
 
'''Important note:''' the action of deleting files with rm '''cannot be undone'''. Use extreme caution when removing files/directories as there is no Ctrl + Z in the terminal.
 
'''Important note:''' the action of deleting files with rm '''cannot be undone'''. Use extreme caution when removing files/directories as there is no Ctrl + Z in the terminal.
  
−
===Logging in and navigating Seawulf===
+
===Logging in and navigating SeaWulf===
 +
 
 +
'''SeaWulf''' is the High Performance Computing (HPC) cluster which researchers and students at Stony Brook University can gain access to use for computationally intensive tasks. Note Stony Brook students often do not automatically have access to SeaWulf, but students can request access from IT or be given access by their labs/instructors. For help getting access to SeaWulf, or questions/problems encountered when using the cluster at any point, '''submit a ticket''' to IT at https://it.stonybrook.edu/help/kb/understanding-seawulf.
 +
 
 +
If one has access to SeaWulf, they can '''access the cluster''' by typing the following into their terminal
 +
 
 +
  ssh -X NETID@login.seawulf.stonybrook.edu
 +
 
 +
Note that the flag -X enables visualizations when running programs in SeaWulf, and NETID is the user's assigned NETID username, used to sign into Blackboard and all associated Stony Brook University domains. The terminal will then prompt the user for the password associated with the entered NETID, and then will give a list of options to use DUO Authentication to grant access. After cooperating with DUO, the user will be in their SeaWulf home directory.
 +
 
 +
To navigate to the '''AMS 536 space''' within SeaWulf, type
 +
 
 +
  cd /gpfs/projects/AMS536
 +
 
 +
To access a '''user's personal scratch space''', which is where users can keep files and test programs (but note all files in this directory will be '''automatically deleted if not modified in over 30 days'''), type
 +
 
 +
  cd /gpfs/scratch/NETID
 +
 
 +
===File transfers to/from SeaWulf===
 +
To '''transfer files''' between the user's local machine and the SeaWulf cluster, this can be done using '''scp'''. For example, to transfer the file "file_1" from the directory Directory_1 on the user's local directory to the user's scratch space on SeaWulf, type the following into the user's local terminal window (not connected to SeaWulf):
 +
 
 +
  scp Directory_1/file_1 NETID@login.seawulf.stonybrook.edu:/gpfs/scratch/NETID
 +
 
 +
Likewise, to transfer file_1 from the scratch space in SeaWulf back to a local directory, type the following into the local terminal window:
  
−
===File transfers to/from Seawulf===
+
  scp NETID@login.seawulf.stonybrook.edu:/gpfs/scratch/NETID/file_1 /Directory_1
  
 
= Chimera: Getting Started =
 
= Chimera: Getting Started =
Line 105: Line 126:
 
(insert picture of PDB w/ search bar, results below after searching)
 
(insert picture of PDB w/ search bar, results below after searching)
  
−
After finding the page for the structure of interest, one can download it in pdb or mol2 (or other) format by selecting the dark blue "Download Files" drop-down menu to the right of the page. Select the desired format, and the download will begin.
+
After finding the page for the structure of interest, one can download it in PDB or Mol2 (or other) format by selecting the dark blue '''Download Files''' drop-down menu to the right of the page. Select the desired format, and the download will begin.
−
 
 
−
===Opening and viewing PDB files===
 
−
(show some screenshots and explain file types)
 
  
 
===First-steps for visualizations in Chimera===
 
===First-steps for visualizations in Chimera===
  
−
====Opening files in Chimera====
+
===Opening files in Chimera===
−
Once one opens Chimera, they are met with a blue screen, and at the bottom right corner are two options called "Browse..." and "Fetch...". If the structure of interest is downloaded already to a local computer, choose "Browse..."; it is possible to download the structure from the PDB directly in Chimera using the "Fetch..." option, but note this requires knowledge of the 4-digit code for the structure listed in the PDB.
+
Once one opens Chimera, they are met with a blue screen, and at the bottom right corner are two options called '''Browse...''' and '''Fetch...'''. If the structure of interest is downloaded already to a local computer, choose '''Browse...'''; it is possible to download the structure from the PDB directly in Chimera using the '''Fetch...''' option, but note this requires knowledge of the 4-digit code for the structure listed in the PDB.
  
 
====Visual configurations====
 
====Visual configurations====
−
Chimera can be a helpful tool for producing publication-ready visualizations of molecules. When Chimera is first opened, the structure will be drawn against a black screen; if a different appearance is desired, click the "Presets" tab at the top, where there are a number of options. A recommended preset is "Publication 1 (silhouette, rounded ribbons)".  
+
Chimera can be a helpful tool for producing publication-ready visualizations of molecules. When Chimera is first opened, the structure will be drawn against a black screen; if a different appearance is desired, click the '''Presets''' tab at the top, where there are a number of options. A recommended preset is '''Publication 1 (silhouette, rounded ribbons)'''.
  
 
====Selecting substructures====
 
====Selecting substructures====
−
Note that it is possible to select portions of the current structure. To do so, click the "Select" tab at the top. Options will become visible to select by chain, residue, structure, and so on. One can also hit "Control" + "Click" to interactively select a small portion of the structure. If one presses the "up" arrow on their keyboard, this will select a larger substructure containing the initial selection. Hitting "up" again will select more, until the entire structure is selected. Alternatively, hitting the "down" arrow will go backwards in selection scope, until the initial small selection is restored.  
+
Note that it is possible to select portions of the current structure. To do so, click the '''Select''' tab at the top. Options will become visible to select by chain, residue, structure, and so on. One can also hit '''Control''' + '''Click''' to interactively select a small portion of the structure. If one presses the "up" arrow on their keyboard, this will select a larger substructure containing the initial selection. Hitting "up" again will select more, until the entire structure is selected. Alternatively, hitting the "down" arrow will go backwards in selection scope, until the initial small selection is restored.
  
 
====Modifying selected substructures====
 
====Modifying selected substructures====
−
One of the easiest ways to make a visualization in Chimera more clear is through color editing. If one is dealing with a protein-ligand complex, for example, it may be useful to color the two structures differently to contrast the two. To do this, select one of the structures (this can be done by interactively clicking and hitting the "up" arrow, or more easily by selecting by structure), then, while the substructure is selected, click the "Actions" tab, click "Color", and make a color choice.  
+
One of the easiest ways to make a visualization in Chimera more clear is through '''color editing'''. If one is dealing with a protein-ligand complex, for example, it may be useful to color the two structures differently to contrast the two. To do this, select one of the structures (this can be done by interactively clicking and hitting the "up" arrow, or more easily by selecting by structure), then, while the substructure is selected, click the '''Actions''' tab, click '''Color''', and make a color choice.  
  
−
If one wishes to isolate specific parts of a structure from the rest, it is necessary to understand how to delete substructures. This can be done by selecting a substructure of interest, clicking "Actions", hovering over "Atoms/Bonds", then under the shown options clicking "delete" at the bottom. This will remove the selected substructure, and cannot be immediately undone, so make sure that parts of the substructure have been colored to make their separation as clear as possible to avoid deleting any parts by accident.
+
If one wishes to isolate specific parts of a structure from the rest, it is necessary to understand how to delete substructures. This can be done by selecting a substructure of interest, clicking "Actions", hovering over '''Atoms/Bonds''', then under the shown options clicking '''delete''' at the bottom. This will remove the selected substructure, and cannot be immediately undone, so make sure that parts of the substructure have been colored to make their separation as clear as possible to avoid deleting any parts by accident.
  
 
===Saving sessions===
 
===Saving sessions===
  
−
If one wishes to save an entire session, with all configurations (such as visual orientation, focus, coloring, etc. which cannot be captures by just saving as a Mol2 or PDB file), one can go to "File" -> "Save Session As..." and this will save the current configurations.  
+
If one wishes to save an entire session, with all configurations (such as visual orientation, focus, coloring, etc. which cannot be captures by just saving as a Mol2 or PDB file), one can go to '''File''' -> '''Save Session As...''' and this will save the current configurations.  
  
 
(put in a picture?)
 
(put in a picture?)
  
 
= Organizing File Directories =  
 
= Organizing File Directories =  
−
It is imperative to one's ability to keep and access records, for any purpose, by having a detailed and organized directory list wherein any necessary files may be kept. While following this docking tutorial, it is recommended to have the following directories:
+
It is imperative to one's ability to keep and access records, for any purpose, by having a '''detailed and organized directory list''' wherein any necessary files may be kept. While following this docking tutorial, it is recommended to have the following directories on both the local machine as well as in the user's designated space on SeaWulf:
  
−
(insert directory list)
+
  000_files
 +
  001_structure
 +
  002_surface_spheres
 +
  003_gridbox
 +
  004_dock
 +
  005_virtual_screen
 +
  006_virtual_screen_mpi
 +
  007_cartesian_min
 +
  008_rescore
  
 
= Dock Prep =
 
= Dock Prep =
  
 
===Initial investigation/preparation===
 
===Initial investigation/preparation===
−
(put in here anything you notice when first opening up a sometimes messy structure--in 6ME2's case we had to remove detergents and water that were far from the of-interest binding site, handle non-standard amino acids, etc. Compare to the structure described in the literature to see if everything is accurate)
+
First, download the 6ME2 file from the PDB, using either Mol2 or PDB format. Save this in '''000_files'''. Then, open this structure using Chimera. It should look as follows:
 +
 
 +
[insert image of PL complex]
 +
 
 +
Inspect the structure of the complex. Note that other non-standard residues besides ramelteon (code "JEV") are present, and will be deleted. Note that to begin the docking preparation process, all solvent must be removed from the file. There is water present in this file, so go to '''Select''' -> '''Residue''' -> '''HOH''', and click '''Actions''' -> '''Atoms/Bonds''' -> '''sphere''' then '''Actions''' -> '''Atoms/Bonds''' -> '''show''', which will show the selected as red spheres, shown below
 +
 
 +
[insert image of HOH selected]
 +
 
 +
Next, hit '''Select''' -> '''Residue''' and select all non-standard residues outside of JEV. While these are selected, go to '''Actions''' -> '''Atoms/Bonds''' -> '''delete''', leaving only the JEV ligand and receptor. Save this as a Mol2 file called '''6me2_noH_noDetergents.mol2''' in the '''000_files''' directory.
 +
 
 +
[insert photo of ligand-protein complex without extra residues]
  
 
==Isolating the ligand and receptor==
 
==Isolating the ligand and receptor==
−
(show how to create two separate files, one with just the ligand and one with just the receptor)
+
To prepare the ligand and receptor for docking, it is necessary to have two separate files--one for the prepared ligand, and one for the prepared receptor, which will be used separately at points.
 +
 
 +
Isolating the receptor:
 +
Load the 6ME2 pdb into Chimera. Select and delete the OLA and PEG residues. ''There are non-standard residues (YCM) which need to be mutated to the closest structural analog, cysteine.'' Using '''Tools > Structure Editing > Rotamers''', change YCM to the most probable rotamers of CYS. Use '''Tools > Surface/Binding analysis > DockPrep''' to add hydrogens to the receptor. Save this structure as a mol2 with '''File > Save mol2'''
 +
 
 +
Isolating ligand:
 +
Select the ligand (JEV) in Chimera. Initially, the structure will present incorrectly with an sp2 carbon after adding hydrogens with the '''Tools > Structure Editing > AddH.
 +
'''
 +
[[File:badjev.png]]
 +
Caption: Incorrect structure of JEV ligand.
  
−
===Prepping the ligand===
+
Select the incorrect atom. Go to '''Tools > Structure Editing > Build Structure'''. In the pop-up window, select “Modify structure”. Under the “Geometry” tab, select tetrahedral, and choose to leave the residue name unchanged, then choose “Apply”. Save this structure as jev_wH.mol2 with '''File > Save mol2'''.
−
(show how to add hydrogens/charges in the isolated ligand)
 
  
−
===Prepping the receptor===
+
[[File:goodjev.png]]
−
(show how to add hydrogens/charges to the isolated receptor)
+
Caption: Correct structure of JEV ligand
  
 
==Surface Spheres==
 
==Surface Spheres==
  
 
===Generating DMS of receptor===
 
===Generating DMS of receptor===
−
(show how to get the surface/save as a DMS file)
+
Open the receptor_woH.mol2 in Chimera. In Chimera, go to '''Actions > Surface > Show''' to visualize the surface of the receptor.
 +
[[File:surfacesssssss.png]]
 +
Save this surface (without the ligand present) using '''Tools > Structure Editing > Write DMS > receptor_woH.dms'''. To visualize the binding pocket, load in the jev_woH.mol2 file, and show the ligand as a surface. For 6ME2, the binding pocket is in the interior of the protein.
  
 
===Creating Surface Spheres===
 
===Creating Surface Spheres===
−
(show how to create the INSPH file in the correct directory, how to run sphgen using the INSPH file, show output).
+
Next, the available space near the active site needs to be identified. These spaces represent possible binding sites for the ligand. DOCK6.9 provides many options here; the first line is the input .dms, the second line, '''R''', tells DOCK to only look for space on the exterior surface of the protein. '''X''' indicates to DOCK to use all surface points. The 4th line, '''0.0''' is the extent of acceptable steric clash. Here, we don't accept any steric clashing, and set it to zero. The 5th and 6th lines, '''4.0''' and '''1.4''' are the maximum and minimum sphere radii, respectively. The final line in the input file is the output .sph, which contains all of the spheres identified by DOCK.
 +
 
 +
'''vi INSPH'''
 +
 
 +
    '''6ME2_rec_surface.dms'''
 +
    '''R'''
 +
    '''X'''
 +
    '''0.0'''
 +
    '''4.0'''
 +
    '''1.4'''
 +
    '''6ME2_receptor_noH.sph'''
 +
 
 +
The following command will call spheregen and generate the spheres
 +
 
 +
'''sphgen -i INSPH -o OUTSPH'''
 +
 
 +
 
 +
If we look at the output in Chimera, we see that there are many spheres generated.
 +
 
 +
[[File:allthebubbles.png]]
  
 
===Selecting Relevant Spheres===
 
===Selecting Relevant Spheres===
−
(show how to run sphere_selector and create necessary input file)
+
Most of these spheres aren't relevant for the ligand in question. Only the spheres that overlap well (within 10 Å) with the docked structure of the ligand are important in this context. To isolate these specific spheres,
 +
 
 +
'''sphere_selector 6ME2_receptor_noH.sph jev_wH.mol2 10.0'''
 +
 
 +
When opening selected_spheres.sph in Chimera, it can be see that most of the spheres have been erased, leaving only the spheres that are within 10 Å of the ligand.
 +
[[File:afewbubbles.png]]
  
 
==Generating Box/Grid==
 
==Generating Box/Grid==
 +
In order to score the ligands upon docking, DOCK uses a grid - with its size defined by a box - to generate a series of points. Energetic contributions will be calculated at each of these points.
 +
===Box Generation===
 +
To create the box, run
 +
'''vi showbox.in '''
 +
The first line, '''Y''' tell the showbox "Yes, create a box". The second line is the distance (Å) from the spheres that the box should be created. The third line is the location of the selected_spheres.sph file, and the last line is the output, which contains the box.
 +
    '''Y'''
 +
    '''8.0'''
 +
    '''../002_surface_spheres/selected_spheres.sph'''
 +
    '''1'''
 +
    '''6me2.box.pdb'''
  
−
===Box Generation===
+
To create the box, run
−
(show how to create showbox.in, how to run showbox)
+
'''showbox < showbox.in'''
  
 
===Grid Generation===
 
===Grid Generation===
−
(show how to create the file grid.in, how to run grid)
+
To make the grid,
 +
'''vi 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                            ../001_structure/jev_wH.mol2'''
 +
  '''box_file                                  6me2.box.pdb'''
 +
  '''vdw_definition_file                      /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/vdw_AMBER_parm99.defn'''
 +
  '''score_grid_prefix                        grid'''
 +
 
 +
Then run,
 +
'''grid -i grid.in -o gridinfo.out'''
 +
If the grid generation was successful, then there will now be gridinfo.out, grid.bmp, and grid.nrg in the directory.
  
 
= Docking =
 
= Docking =
  
 
===Energy Minimization===
 
===Energy Minimization===
 +
 +
The purpose of using DOCK for this tutorial is to find the most energetically favorable position of the ligand in the target site pocket. To get as good of a DOCK score as possible, it is pertinent to minimize the potential energy of the ligand itself.
 +
 +
Switch to the '''004_dock''' directory by typing
 +
 +
  cd ../004_dock
 +
 +
Then create an input file called '''min.in''' by typing
 +
 +
  vi min.in
 +
 +
Then copy and paste the following into this file:
 +
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100.0
 +
  ligand_atom_file                                            ../001_structure/6me2_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/6me2_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
 +
  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_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                                        6me2.ligand.min
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        1
 +
  rank_ligands                                                no
 +
 +
Next, run the energy minimization by typing
 +
 +
  dock6 -i min.in -o min.out
 +
 +
When the produced '''6me2.ligand.min_scored.mol2''' is visualized in Chimera, it looks as follows:
 +
 +
[[File:Lig min score.png||1000px]]
 +
 +
To view the RMSD of the minimized ligand conformation, vi into this file:
 +
 +
  vi 6me2.ligand.min_scored.mol2
  
 
==Rigid Docking==
 
==Rigid Docking==
 +
 +
Note that there are multiple different docking protocols that this tutorial will utilize--'''rigid, fixed anchor,''' and '''flexible'''. The first this tutorial will utilize is '''rigid''' docking.
 +
 +
Make sure to still be in the '''004_dock''' directory, and create a rigid docking input file as follows:
 +
 +
  vi rigid.in
 +
 +
Copy and paste the following input parameters into the file:
 +
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          yes
 +
  internal_energy_rep_exp                                      12
 +
  internal_energy_cutoff                                      100
 +
  ligand_atom_file                                            6me2.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                                      6me2.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
 +
 +
Now run the rigid docking protocol using this input file by typing:
 +
 +
  dock6 -i rigid.in
 +
 +
Note this will produce the output file '''rigid.out_scored.mol2'''.
 +
 +
Be sure to open the file and inspect the RMSD by typing:
 +
 +
  vi rigid.out_scored.mol2
  
 
==Fixed Anchor Docking==
 
==Fixed Anchor Docking==
 +
 +
The next docking protocol used is the '''fixed anchor''' protocol, wherein the location of the ligand in the binding pocket is fixed but the rotatable bonds are moved during energy minimization. Make sure to still be in the '''004_dock''' directory, and create the following input file
 +
 +
  vi fixed.in
 +
 +
And copy and paste the following parameters into this file:
 +
 +
  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/6me2_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/6me2_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                                        6me2_fad
 +
  write_orientations                                          no
 +
  num_scored_conformers                                        100
 +
  write_conformations                                          no
 +
  cluster_conformations                                        yes
 +
  cluster_rmsd_threshold                                      2.0
 +
  rank_ligands                                                no
 +
 +
And now run the docking protocol by typing:
 +
 +
  dock6 -i fixed.in -o fixed.out
 +
 +
This produces the output file '''6me2_fad_scored.mol2'''.
 +
 +
And be sure to check the RMSD value of the docked ligand by typing:
 +
 +
  vi 6me2_fad_scored.mol2
  
 
==Flexible Docking==
 
==Flexible Docking==
 +
 +
The final docking protocol that will be performed in this tutorial is '''flexible''' docking. In this method, the ligand is completely free to move with all degrees of freedom being taken into account; because of the added number of possible conformations in this protocol, it is the most computationally intensive out of the three tested herein.
 +
 +
Be sure to remain in the '''004_dock''' folder, and create the following input file:
 +
 +
  vi flex.in
 +
 +
And paste the following parameters into this file:
 +
 +
  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                                            6me2.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                                      6me2.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
 +
 +
Now run the flexible docking protocol with the input file by typing the following:
 +
 +
  dock6 -i flex.in
 +
 +
When the output file '''flex.out_scored.mol2''' is opened in Chimera it appears as follows:
 +
 +
[insert image]
 +
 +
Be sure to open up the file and check the RMSD by typing:
 +
 +
  vi flex.out_scored.mol2
  
 
==Calculating Molecular Footprint==
 
==Calculating Molecular Footprint==
  
−
= Virtual Screening =  
+
The '''molecular footprint''' is a helpful tool to visually inspect and analyze the van der Waals and electrostatic interactions between the ligand and local amino acid residues, which can imply affinity of the ligand to the receptor on a per-residue basis. The footprint can also give clues to structural characteristics, such as the appearance of hydrogen bonds at peaks in the footprint graph, and thus can be used to find errors in the ligand conformation during energy minimization.
 +
 
 +
To calculate the molecular footprint, create an input file in the '''004_dock''' directory as follows:
 +
 
 +
  vi footprint.in
 +
 
 +
And paste the following into this file:
 +
 
 +
  conformer_search_type                                        rigid
 +
  use_internal_energy                                          no
 +
  ligand_atom_file                                            6me2.ligand.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                  ../001.structure/6me2_ligand_wH.mol2
 +
  fps_score_foot_compare_type                                  Euclidean
 +
  fps_score_normalize_foot                                    no
 +
  fps_score_foot_comp_all_residue                              yes
 +
  fps_score_receptor_filename                                  ../001.structure/6me2_receptor_dockprep.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
 +
 
 +
Now run the footprint calculation by typing the following:
 +
 
 +
  dock6 -i footprint.in -o footprint.out
 +
 
 +
There will now be a python script created in the '''004_dock''' folder, called '''plot_footprint_single_magnitude.py''', which must be run to produce a visualization of the footprint calculation:
 +
 
 +
  python plot_footprint_single_magnitude.py footprint.out_footprint_scored.txt  50
 +
 
 +
And this will produce a PDF containing a visualization of the footprint, called '''footprint.out.pdf'''. Open this in a PDF viewer, and it appears as follows:
 +
 
 +
[[File:Mol foot 6me2.png||1000px]]
 +
 
 +
= Virtual Screening =
 +
 
 +
Now that we have finished preparing both the ligand and receptor and verified that the ligand can be reasonably docked back into the receptor, we can now begin with the virtual screen. The idea behind virtual screening is that you can take a library of known molecules, dock them to a receptor, and then score how well those molecules are binding to the receptor. In this tutorial, our library is a file named "VS_library_5K.mol2", which contains 5,000 small molecules for us to dock to 6ME2. Let's begin by generating that library file of 5,000 molecules. Cd into the directory "005_virtual_screen" and copy the "VS_library_25K.mol2" file from "/gpfs/projects/AMS536/zzz.programs" into that directory with the following command:
 +
 
 +
cp /gpfs/projects/AMS536/zzz.programs/VS_library_25K.mol2 ./
 +
 
 +
Now we need to truncate this library for the top 5,000 molecules. We can do this with a very simple python script. Vi into "truncate_library.py" and write the following to that file (you can delete the comments, but they will give context to what is happening line-by-line if you don't know):
 +
 
 +
#Keep track of how many molecules have been copied over.
 +
moleculenumber=0                                                       
 +
#Create a list to hold each line being copied.
 +
libraryfile = []                                                 
 +
#Open the 25K library file, and number each line so it can be read line-by-line.               
 +
read_file = open("VS_library_25K.mol2")
 +
lines = read_file.readlines()
 +
#Add each individual line to the list.
 +
for i, line in enumerate(lines):
 +
        libraryfile.append(line)
 +
#Create the new library file so that molecules can be written to it.
 +
with open("VS_library_5K.mol2", "w") as new_f:
 +
#Copy over each line in the list for the first 5,000 molecules
 +
        for d in libraryfile:
 +
                new_f.write(d)
 +
                if "MOLECULE" in d:
 +
                        moleculenumber += 1
 +
                if moleculenumber == 5001:
 +
                        break
 +
 
 +
Once you have finished doing that, save and exit the file. To run this, type the following command into the terminal:
 +
 
 +
python truncate_library.py
 +
 
 +
Once you have your library of molecules, the next thing we want to do is write the input file. Vi into virtual.in and write the following to that file:
 +
 
 +
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                                          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                                            VS_library_5K.mol2
 +
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                                          ../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                                        virtual.out
 +
write_orientations                                          no
 +
num_scored_conformers                                        1
 +
rank_ligands                                                no
 +
 
 +
Once you have finished doing that, save and exit the file. Now, we could try running the virtual screen similarly to previous steps, but it would be really slow. The best way to do this is to parallelize the task as many times as possible, and to do this we are going to use MPI.
 +
 
 +
= Virtual Screening Using MPI =
 +
 
 +
Running the virtual screen using MPI does not require us to change anything about in input file, but the submit file will look a little different than the previous one you made for the flexible docking step. To start, cd into the "006_virtual_screen_mpi" directory and copy over the 5,000 molecule library we generated. Vi into "virtualscreen.sh" and write the following to that file:
 +
 
 +
#!/bin/bash
 +
#SBATCH --time=48:00:00
 +
#SBATCH --nodes=1
 +
#SBATCH --ntasks=28
 +
#SBATCH --job-name=6ME2_VS_min
 +
#SBATCH --output=6ME2_VS_min
 +
#SBATCH -p long-28core
 +
 
 +
cd $SLURM_SUBMIT_DIR
 +
 
 +
mpirun -np 28 dock6.mpi -i ../005_virtual_screen/virtual.in -o 6ME2.virtual.mpi.out
 +
 
 +
This will run the virtual screen across 28 cores which will significantly speed up the virtual screen process. To run this even faster, you can increase the number of nodes from one to up to 4, which would increase the total number of cores to either 56, 84, or 112. If you elect to try with more nodes, make sure to adjust both the number of nodes (line 3) and the total number of cores (line 11) in the submit file. To submit this simulation to the queue, type the following command into the terminal:
 +
 
 +
sbatch virtualscreen.sh
 +
 
 +
This will take awhile to finish, even with the MPI helping to speed things up. Once it's done, make sure to check that everything ran to completion per your output files. You should then, using Chimera, load the receptor and ligand files and see how the molecules from the virtual screen align with the ligand. The easiest way to do that is not to load in every molecule at once but to use a program in Chimera called ViewDock. ViewDock will let you visualize each molecule one at a time and order them based on parameters scuh as molecular weight or docking score. You can find ViewDock in (Tools -> Surface/Binding Analysis -> ViewDock). Select the output file from the virtual screen and from the next menu select "DOCK 4, 5, or 6" as we are using DOCK6. You can then scroll through the different molecules and see how well they align.
 +
 
 +
= Cartesian Minimization of Docked Molecules =
 +
 
 +
The structures you just looked at were the raw output from the virtual screen. It would be better, before making a conclusion on which of the 5,000 molecules bind the best, to run an energy minimization on each of those molecules and then rescore them based on any changes in their interactions with the receptor. Starting with the energy minimization, change directories to 007_cartesian_min. We will need to write an input and submit file for this simulation, so let's start with the input file and vi into cartmin.in. Write the following to that file:
  
−
= Virtual Screening Using MPI =
+
conformer_search_type                                        rigid
 +
use_internal_energy                                          yes
 +
internal_energy_rep_exp                                      12
 +
internal_energy_cutoff                                      100
 +
ligand_atom_file                                            ../006_virtualscreenmpi/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                                    yes
 +
continuous_score_secondary                                  no
 +
cont_score_rec_filename                                      ../001_structure/6ME2.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                                        6ME2.virtual_screen.min
 +
write_orientations                                          no
 +
num_scored_conformers                                        1
 +
rank_ligands                                                no
  
−
= Cartesian Minimization of Docked Molecules =
+
Once you finish save and exit that file and then we can start to write the submit file. Vi into cartmin.sh and write the following to that file:
 +
 
 +
#!/bin/bash
 +
#SBATCH --time=48:00:00
 +
#SBATCH --nodes=1
 +
#SBATCH --ntasks=28
 +
#SBATCH --job-name=mpi_6ME2_VS_AR
 +
#SBATCH --output=6ME2_VS_mpi.out
 +
#SBATCH -p long-28core
 +
 
 +
cd $SLURM_SUBMIT_DIR
 +
 
 +
mpirun -np 28 /gpfs/projects/AMS536/zzz.programs/dock6.9_release/bin/dock6.mpi -i cartmin.in -o cartmin.out
 +
 
 +
We are using mpi again to speed up the energy minimization, so if you choose to adjust the number of nodes being used to speed it up even further, adjust the file accordingly as mentioned in the virtual screen using MPI section. Once you are happy with the submit file, you can start it by typing the following command into the terminal:
 +
 
 +
sbatch cartmin.sh
 +
 
 +
Once this finishes running, check the output files to make sure everything ran smoothly. Once that is done, we can now move to rescore these molecules and then visualize them.
  
 
= Rescoring Docked Molecules =
 
= Rescoring Docked Molecules =
 +
 +
Now that we have minimized all the molecules docked in the virtual screen, we can rescore them based on their new conformations. In order to do this, change directories to 008_rescore. We will need to start by creating the input and submit files, so vi into rescore.in and write the following to that file:
 +
 +
conformer_search_type                                        rigid
 +
use_internal_energy                                          yes
 +
internal_energy_rep_exp                                      12
 +
internal_energy_cutoff                                      100
 +
ligand_atom_file                                            ../007_cartesian_min/6ME2.virtual_screen.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                          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                                      no
 +
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/6ME2.lig.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/6ME2.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/6ME2.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_hms_score_ref_filename                            ../004_dock/6ME2.lig.min_scored.mol2
 +
descriptor_hms_score_matching_coeff                          -5
 +
descriptor_hms_score_rmsd_coeff                              1
 +
descriptor_volume_score_reference_mol2_filename              ../004_dock/6ME2.lig.min_scored.mol2
 +
descriptor_volume_score_overlap_compute_method              analytical
 +
descriptor_weight_fps_score                                  1
 +
descriptor_weight_pharmacophore_score                        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
 +
 +
Once you have finished, save and exit the file. Then vi rescore.sh and write the following to that file:
 +
 +
#!/bin/sh
 +
#SBATCH --partition=long-28core
 +
#SBATCH --time=48:00:00
 +
#SBATCH --nodes=1
 +
#SBATCH --ntasks=28
 +
#SBATCH --job-name=6ME2-rescore
 +
#SBATCH --output=%x-%j.o
 +
 +
echo "starting DOCK6.9 Simulation"
 +
module load intel/mpi/64/2018/18.0.3
 +
 +
mpirun -np 28 dock6.mpi -i rescore.in -o rescore.out
 +
echo "done"
 +
 +
Again, we are using mpi so if you want to speed things up use for than one node and then adjust the submit file accordingly. Once you have finished, save and exit the file. Submit the job to the queue as described above. Once this has finished, check your input files to make sure everything ran smoothly. You will have generated multiple files named rescore.out.x, but every molecule described in the output files can also be found in the descriptor.output_scored.mol2 file.
 +
 +
Now, we can properly visualize the minimized and rescored molecules. Transfer the descriptor.output_scored.mol2 file to your computer and open Chimera. Load into Chimera the ligand and receptor files (6ME2.receptor_wH.mol2 and 6ME2.ligand_wH.mol2). Additionally, load in the the gridbox generated in showbox (6ME2.gridbox.pdb). Your screen should look as follows:
 +
 +
[[File:Receptor+ligand+showbox.png]]
 +
 +
Now, using DockView as described above: load in the descriptor.output_scored.mol2 file and begin moving through different ligands one at a time. Which ligands have the highest molecular weight, or the best score from DOCK? Using the ViewDock panel, go to (Column -> Show -> Molecular Weight) and a new column will appear in the panel which you can use to sort the molecules on any number of properties listed in that submenu. Once you have finished analyzing your different ligands, you have finished this part of the DOCK tutorial for 6ME2. Congratulations!

Latest revision as of 13:30, 20 February 2023

Required Software

UCSF Chimera, UCSF DOCK 6.9

Chimera

UCSF Chimera is a software which aids the interactive 3D visualization, editing, and analysis of biological macromolecules.

[insert screenshot of chimera]

Chimera will be used in this tutorial to prepare molecules for docking, as well as for simply exploring the structures of molecules of interest.

DOCK

UCSF DOCK is a molecular docking program, used to find the optimal orientation of a ligand binding to a receptor. DOCK can be used for a number of different molecular modeling/design tasks, such as testing binding affinities, virtual screening, and de novo design.

6ME2

6ME2 is the modeled structure of receptor MT1 of melatonin complexed with ramelteon (the ligand). Further details on the structure from the original researchers can be found at https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6696938/

Seawulf: Getting Started

Basics of the Command Line/Terminal

There are some basic commands for the command line/terminal that are necessary to know in order to complete this tutorial effectively.

When one first opens a command line window, it will begin in their home directory. To change this directory to, say, Desktop, type: 

 cd ~/Desktop

To get to a folder, say folder_1, in Desktop, type 

 cd ~/Desktop/folder_1

To list out all of the directories/files contained within this current directory, type 

 ls

To change the directory to one folder in the hierarchy above the current, type 

 cd ../

to change directories a folder, say folder_2 inside the current directory, type 

 cd ./folder_2

To create a file in the current directory, called "file_example", using vim, type 

 vi filename_example

This will open up a prompt where one can type in, paste, or edit (if the file already exists in the current directory) the contents of the file. To edit the file once opened, type "i", which will change the environment to "INSERT" mode. "INSERT" is displayed at the bottom of the window when in this mode. To exit "INSERT" mode, press "Escape". If in "INSERT" mode, to save the file, press "Escape" and then type 

 :wq

which will save the file and open the terminal window back up to commands, in the current directory. If one types 

 ls

the newly created "file_example" should be visible in the outputted list of files.

To move files between directories, say "file_example" is the file of interest, from Directory_1 to Directory_2, type 

 mv /Directory_1/file_example /Directory_2/file_example

The name of file_example may also be changed (to "file_new_name") during this move: 

 mv /Directory_1/file_example /Directory_2/file_new_name

Note this should be used to rename directories--move the file to the same directory and change the name to rename the file.

To copy files between directories, type the following 

 cp /Directory_1/file_example /Directory_2/file_new_name

And one can copy directories as well, by typing 

 cp -R /Directory_1/sub_directory /Directory_2/sub_directory_new_name

Note that the "-R" flag specifies to copy over the sub_directory folder and all of its contents from Directory_1 to Directory_2.

To create a new directory, called "new_directory" for example, type 

 mkdir new_directory

To delete a file in the current directory, "file_example", type 

 rm file_example

To delete a directory called "directory_example", including all of its subdirectories and associated files, type 

 rm -r

Note that if one wishes to perform the same action on multiple files at once, this can be done on the same line by separating the filenames with a space. For example to delete the files filename_1 in Directory_1 and filename_2 in Directory_2, type 

 rm Directory_1/filename_1 Directory_2/filename_2

Important note: the action of deleting files with rm cannot be undone. Use extreme caution when removing files/directories as there is no Ctrl + Z in the terminal.

Logging in and navigating SeaWulf

SeaWulf is the High Performance Computing (HPC) cluster which researchers and students at Stony Brook University can gain access to use for computationally intensive tasks. Note Stony Brook students often do not automatically have access to SeaWulf, but students can request access from IT or be given access by their labs/instructors. For help getting access to SeaWulf, or questions/problems encountered when using the cluster at any point, submit a ticket to IT at https://it.stonybrook.edu/help/kb/understanding-seawulf.

If one has access to SeaWulf, they can access the cluster by typing the following into their terminal 

 ssh -X NETID@login.seawulf.stonybrook.edu

Note that the flag -X enables visualizations when running programs in SeaWulf, and NETID is the user's assigned NETID username, used to sign into Blackboard and all associated Stony Brook University domains. The terminal will then prompt the user for the password associated with the entered NETID, and then will give a list of options to use DUO Authentication to grant access. After cooperating with DUO, the user will be in their SeaWulf home directory.

To navigate to the AMS 536 space within SeaWulf, type 

 cd /gpfs/projects/AMS536

To access a user's personal scratch space, which is where users can keep files and test programs (but note all files in this directory will be automatically deleted if not modified in over 30 days), type 

 cd /gpfs/scratch/NETID

File transfers to/from SeaWulf

To transfer files between the user's local machine and the SeaWulf cluster, this can be done using scp. For example, to transfer the file "file_1" from the directory Directory_1 on the user's local directory to the user's scratch space on SeaWulf, type the following into the user's local terminal window (not connected to SeaWulf): 

 scp Directory_1/file_1 NETID@login.seawulf.stonybrook.edu:/gpfs/scratch/NETID

Likewise, to transfer file_1 from the scratch space in SeaWulf back to a local directory, type the following into the local terminal window: 

 scp NETID@login.seawulf.stonybrook.edu:/gpfs/scratch/NETID/file_1 /Directory_1

Chimera: Getting Started

Searching for/Downloading files from the PDB

The RCSB Protein Data Bank is a useful tool for searching for and downloading files containing the 3D structures of biological macromolecules. One can search for a molecule by using an exact 4 letter code (such as "6ME2") or by searching for a short description (such as "melatonin with ramelteon complex"). Note that there may be many molecules that fit a given description, so one may need to scroll through results shown after the search to find the one of interest.

(insert picture of PDB w/ search bar, results below after searching)

After finding the page for the structure of interest, one can download it in PDB or Mol2 (or other) format by selecting the dark blue Download Files drop-down menu to the right of the page. Select the desired format, and the download will begin.

First-steps for visualizations in Chimera

Opening files in Chimera

Once one opens Chimera, they are met with a blue screen, and at the bottom right corner are two options called Browse... and Fetch.... If the structure of interest is downloaded already to a local computer, choose Browse...; it is possible to download the structure from the PDB directly in Chimera using the Fetch... option, but note this requires knowledge of the 4-digit code for the structure listed in the PDB.

Visual configurations

Chimera can be a helpful tool for producing publication-ready visualizations of molecules. When Chimera is first opened, the structure will be drawn against a black screen; if a different appearance is desired, click the Presets tab at the top, where there are a number of options. A recommended preset is Publication 1 (silhouette, rounded ribbons).

Selecting substructures

Note that it is possible to select portions of the current structure. To do so, click the Select tab at the top. Options will become visible to select by chain, residue, structure, and so on. One can also hit Control + Click to interactively select a small portion of the structure. If one presses the "up" arrow on their keyboard, this will select a larger substructure containing the initial selection. Hitting "up" again will select more, until the entire structure is selected. Alternatively, hitting the "down" arrow will go backwards in selection scope, until the initial small selection is restored.

Modifying selected substructures

One of the easiest ways to make a visualization in Chimera more clear is through color editing. If one is dealing with a protein-ligand complex, for example, it may be useful to color the two structures differently to contrast the two. To do this, select one of the structures (this can be done by interactively clicking and hitting the "up" arrow, or more easily by selecting by structure), then, while the substructure is selected, click the Actions tab, click Color, and make a color choice.

If one wishes to isolate specific parts of a structure from the rest, it is necessary to understand how to delete substructures. This can be done by selecting a substructure of interest, clicking "Actions", hovering over Atoms/Bonds, then under the shown options clicking delete at the bottom. This will remove the selected substructure, and cannot be immediately undone, so make sure that parts of the substructure have been colored to make their separation as clear as possible to avoid deleting any parts by accident.

Saving sessions

If one wishes to save an entire session, with all configurations (such as visual orientation, focus, coloring, etc. which cannot be captures by just saving as a Mol2 or PDB file), one can go to File -> Save Session As... and this will save the current configurations.

(put in a picture?)

Organizing File Directories

It is imperative to one's ability to keep and access records, for any purpose, by having a detailed and organized directory list wherein any necessary files may be kept. While following this docking tutorial, it is recommended to have the following directories on both the local machine as well as in the user's designated space on SeaWulf: 

 000_files
 001_structure
 002_surface_spheres
 003_gridbox
 004_dock
 005_virtual_screen
 006_virtual_screen_mpi
 007_cartesian_min
 008_rescore

Dock Prep

Initial investigation/preparation

First, download the 6ME2 file from the PDB, using either Mol2 or PDB format. Save this in 000_files. Then, open this structure using Chimera. It should look as follows:

[insert image of PL complex]

Inspect the structure of the complex. Note that other non-standard residues besides ramelteon (code "JEV") are present, and will be deleted. Note that to begin the docking preparation process, all solvent must be removed from the file. There is water present in this file, so go to Select -> Residue -> HOH, and click Actions -> Atoms/Bonds -> sphere then Actions -> Atoms/Bonds -> show, which will show the selected as red spheres, shown below

[insert image of HOH selected]

Next, hit Select -> Residue and select all non-standard residues outside of JEV. While these are selected, go to Actions -> Atoms/Bonds -> delete, leaving only the JEV ligand and receptor. Save this as a Mol2 file called 6me2_noH_noDetergents.mol2 in the 000_files directory.

[insert photo of ligand-protein complex without extra residues]

Isolating the ligand and receptor

To prepare the ligand and receptor for docking, it is necessary to have two separate files--one for the prepared ligand, and one for the prepared receptor, which will be used separately at points.

Isolating the receptor: Load the 6ME2 pdb into Chimera. Select and delete the OLA and PEG residues. There are non-standard residues (YCM) which need to be mutated to the closest structural analog, cysteine. Using Tools > Structure Editing > Rotamers, change YCM to the most probable rotamers of CYS. Use Tools > Surface/Binding analysis > DockPrep to add hydrogens to the receptor. Save this structure as a mol2 with File > Save mol2

Isolating ligand: Select the ligand (JEV) in Chimera. Initially, the structure will present incorrectly with an sp2 carbon after adding hydrogens with the Tools > Structure Editing > AddH. Badjev.png Caption: Incorrect structure of JEV ligand.

Select the incorrect atom. Go to Tools > Structure Editing > Build Structure. In the pop-up window, select “Modify structure”. Under the “Geometry” tab, select tetrahedral, and choose to leave the residue name unchanged, then choose “Apply”. Save this structure as jev_wH.mol2 with File > Save mol2.

Goodjev.png Caption: Correct structure of JEV ligand

Surface Spheres

Generating DMS of receptor

Open the receptor_woH.mol2 in Chimera. In Chimera, go to Actions > Surface > Show to visualize the surface of the receptor. Surfacesssssss.png Save this surface (without the ligand present) using Tools > Structure Editing > Write DMS > receptor_woH.dms. To visualize the binding pocket, load in the jev_woH.mol2 file, and show the ligand as a surface. For 6ME2, the binding pocket is in the interior of the protein.

Creating Surface Spheres

Next, the available space near the active site needs to be identified. These spaces represent possible binding sites for the ligand. DOCK6.9 provides many options here; the first line is the input .dms, the second line, R, tells DOCK to only look for space on the exterior surface of the protein. X indicates to DOCK to use all surface points. The 4th line, 0.0 is the extent of acceptable steric clash. Here, we don't accept any steric clashing, and set it to zero. The 5th and 6th lines, 4.0 and 1.4 are the maximum and minimum sphere radii, respectively. The final line in the input file is the output .sph, which contains all of the spheres identified by DOCK.

vi INSPH 

   6ME2_rec_surface.dms
   R
   X
   0.0
   4.0
   1.4
   6ME2_receptor_noH.sph

The following command will call spheregen and generate the spheres

sphgen -i INSPH -o OUTSPH


If we look at the output in Chimera, we see that there are many spheres generated.

Allthebubbles.png

Selecting Relevant Spheres

Most of these spheres aren't relevant for the ligand in question. Only the spheres that overlap well (within 10 Å) with the docked structure of the ligand are important in this context. To isolate these specific spheres,

sphere_selector 6ME2_receptor_noH.sph jev_wH.mol2 10.0

When opening selected_spheres.sph in Chimera, it can be see that most of the spheres have been erased, leaving only the spheres that are within 10 Å of the ligand. Afewbubbles.png

Generating Box/Grid

In order to score the ligands upon docking, DOCK uses a grid - with its size defined by a box - to generate a series of points. Energetic contributions will be calculated at each of these points.

Box Generation

To create the box, run vi showbox.in The first line, Y tell the showbox "Yes, create a box". The second line is the distance (Å) from the spheres that the box should be created. The third line is the location of the selected_spheres.sph file, and the last line is the output, which contains the box. 

   Y
   8.0
   ../002_surface_spheres/selected_spheres.sph
   1
   6me2.box.pdb

To create the box, run showbox < showbox.in

Grid Generation

To make the grid, vi 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                             ../001_structure/jev_wH.mol2
 box_file                                  6me2.box.pdb
 vdw_definition_file                       /gpfs/projects/AMS536/zzz.programs/dock6.9_release/parameters/vdw_AMBER_parm99.defn
 score_grid_prefix                         grid

Then run, grid -i grid.in -o gridinfo.out If the grid generation was successful, then there will now be gridinfo.out, grid.bmp, and grid.nrg in the directory.

Docking

Energy Minimization

The purpose of using DOCK for this tutorial is to find the most energetically favorable position of the ligand in the target site pocket. To get as good of a DOCK score as possible, it is pertinent to minimize the potential energy of the ligand itself.

Switch to the 004_dock directory by typing 

 cd ../004_dock

Then create an input file called min.in by typing 

 vi min.in

Then copy and paste the following into this file: 

 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100.0
 ligand_atom_file                                             ../001_structure/6me2_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/6me2_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
 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_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                                        6me2.ligand.min
 write_orientations                                           no
 num_scored_conformers                                        1
 rank_ligands                                                 no

Next, run the energy minimization by typing 

 dock6 -i min.in -o min.out

When the produced 6me2.ligand.min_scored.mol2 is visualized in Chimera, it looks as follows:

Lig min score.png

To view the RMSD of the minimized ligand conformation, vi into this file: 

 vi 6me2.ligand.min_scored.mol2

Rigid Docking

Note that there are multiple different docking protocols that this tutorial will utilize--rigid, fixed anchor, and flexible. The first this tutorial will utilize is rigid docking.

Make sure to still be in the 004_dock directory, and create a rigid docking input file as follows: 

 vi rigid.in

Copy and paste the following input parameters into the file: 

 conformer_search_type                                        rigid
 use_internal_energy                                          yes
 internal_energy_rep_exp                                      12
 internal_energy_cutoff                                       100
 ligand_atom_file                                             6me2.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                                      6me2.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

Now run the rigid docking protocol using this input file by typing: 

 dock6 -i rigid.in

Note this will produce the output file rigid.out_scored.mol2.

Be sure to open the file and inspect the RMSD by typing: 

 vi rigid.out_scored.mol2

Fixed Anchor Docking

The next docking protocol used is the fixed anchor protocol, wherein the location of the ligand in the binding pocket is fixed but the rotatable bonds are moved during energy minimization. Make sure to still be in the 004_dock directory, and create the following input file 

 vi fixed.in

And copy and paste the following parameters into this file: 

 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/6me2_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/6me2_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                                        6me2_fad
 write_orientations                                           no
 num_scored_conformers                                        100
 write_conformations                                          no
 cluster_conformations                                        yes
 cluster_rmsd_threshold                                       2.0
 rank_ligands                                                 no

And now run the docking protocol by typing: 

 dock6 -i fixed.in -o fixed.out

This produces the output file 6me2_fad_scored.mol2.

And be sure to check the RMSD value of the docked ligand by typing: 

 vi 6me2_fad_scored.mol2

Flexible Docking

The final docking protocol that will be performed in this tutorial is flexible docking. In this method, the ligand is completely free to move with all degrees of freedom being taken into account; because of the added number of possible conformations in this protocol, it is the most computationally intensive out of the three tested herein.

Be sure to remain in the 004_dock folder, and create the following input file: 

 vi flex.in

And paste the following parameters into this file: 

 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                                             6me2.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                                      6me2.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

Now run the flexible docking protocol with the input file by typing the following: 

 dock6 -i flex.in

When the output file flex.out_scored.mol2 is opened in Chimera it appears as follows:

[insert image]

Be sure to open up the file and check the RMSD by typing: 

 vi flex.out_scored.mol2

Calculating Molecular Footprint

The molecular footprint is a helpful tool to visually inspect and analyze the van der Waals and electrostatic interactions between the ligand and local amino acid residues, which can imply affinity of the ligand to the receptor on a per-residue basis. The footprint can also give clues to structural characteristics, such as the appearance of hydrogen bonds at peaks in the footprint graph, and thus can be used to find errors in the ligand conformation during energy minimization.

To calculate the molecular footprint, create an input file in the 004_dock directory as follows: 

 vi footprint.in

And paste the following into this file: 

 conformer_search_type                                        rigid
 use_internal_energy                                          no
 ligand_atom_file                                             6me2.ligand.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                  ../001.structure/6me2_ligand_wH.mol2
 fps_score_foot_compare_type                                  Euclidean
 fps_score_normalize_foot                                     no
 fps_score_foot_comp_all_residue                              yes
 fps_score_receptor_filename                                  ../001.structure/6me2_receptor_dockprep.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

Now run the footprint calculation by typing the following: 

 dock6 -i footprint.in -o footprint.out

There will now be a python script created in the 004_dock folder, called plot_footprint_single_magnitude.py, which must be run to produce a visualization of the footprint calculation: 

 python plot_footprint_single_magnitude.py footprint.out_footprint_scored.txt  50

And this will produce a PDF containing a visualization of the footprint, called footprint.out.pdf. Open this in a PDF viewer, and it appears as follows:

Mol foot 6me2.png

Virtual Screening

Now that we have finished preparing both the ligand and receptor and verified that the ligand can be reasonably docked back into the receptor, we can now begin with the virtual screen. The idea behind virtual screening is that you can take a library of known molecules, dock them to a receptor, and then score how well those molecules are binding to the receptor. In this tutorial, our library is a file named "VS_library_5K.mol2", which contains 5,000 small molecules for us to dock to 6ME2. Let's begin by generating that library file of 5,000 molecules. Cd into the directory "005_virtual_screen" and copy the "VS_library_25K.mol2" file from "/gpfs/projects/AMS536/zzz.programs" into that directory with the following command: 

cp /gpfs/projects/AMS536/zzz.programs/VS_library_25K.mol2 ./

Now we need to truncate this library for the top 5,000 molecules. We can do this with a very simple python script. Vi into "truncate_library.py" and write the following to that file (you can delete the comments, but they will give context to what is happening line-by-line if you don't know): 

#Keep track of how many molecules have been copied over.
moleculenumber=0                                                         
#Create a list to hold each line being copied.
libraryfile = []                                                   
#Open the 25K library file, and number each line so it can be read line-by-line.                 
read_file = open("VS_library_25K.mol2")
lines = read_file.readlines()
#Add each individual line to the list.
for i, line in enumerate(lines):
        libraryfile.append(line)
#Create the new library file so that molecules can be written to it.
with open("VS_library_5K.mol2", "w") as new_f:
#Copy over each line in the list for the first 5,000 molecules
        for d in libraryfile:
                new_f.write(d)
                if "MOLECULE" in d:
                        moleculenumber += 1
                if moleculenumber == 5001:
                        break

Once you have finished doing that, save and exit the file. To run this, type the following command into the terminal: 

python truncate_library.py

Once you have your library of molecules, the next thing we want to do is write the input file. Vi into virtual.in and write the following to that file: 

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                                          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                                             VS_library_5K.mol2
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                                           ../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                                        virtual.out 
write_orientations                                           no
num_scored_conformers                                        1
rank_ligands                                                 no

Once you have finished doing that, save and exit the file. Now, we could try running the virtual screen similarly to previous steps, but it would be really slow. The best way to do this is to parallelize the task as many times as possible, and to do this we are going to use MPI.

Virtual Screening Using MPI

Running the virtual screen using MPI does not require us to change anything about in input file, but the submit file will look a little different than the previous one you made for the flexible docking step. To start, cd into the "006_virtual_screen_mpi" directory and copy over the 5,000 molecule library we generated. Vi into "virtualscreen.sh" and write the following to that file: 

#!/bin/bash
#SBATCH --time=48:00:00
#SBATCH --nodes=1
#SBATCH --ntasks=28
#SBATCH --job-name=6ME2_VS_min
#SBATCH --output=6ME2_VS_min
#SBATCH -p long-28core

cd $SLURM_SUBMIT_DIR

mpirun -np 28 dock6.mpi -i ../005_virtual_screen/virtual.in -o 6ME2.virtual.mpi.out

This will run the virtual screen across 28 cores which will significantly speed up the virtual screen process. To run this even faster, you can increase the number of nodes from one to up to 4, which would increase the total number of cores to either 56, 84, or 112. If you elect to try with more nodes, make sure to adjust both the number of nodes (line 3) and the total number of cores (line 11) in the submit file. To submit this simulation to the queue, type the following command into the terminal: 

sbatch virtualscreen.sh

This will take awhile to finish, even with the MPI helping to speed things up. Once it's done, make sure to check that everything ran to completion per your output files. You should then, using Chimera, load the receptor and ligand files and see how the molecules from the virtual screen align with the ligand. The easiest way to do that is not to load in every molecule at once but to use a program in Chimera called ViewDock. ViewDock will let you visualize each molecule one at a time and order them based on parameters scuh as molecular weight or docking score. You can find ViewDock in (Tools -> Surface/Binding Analysis -> ViewDock). Select the output file from the virtual screen and from the next menu select "DOCK 4, 5, or 6" as we are using DOCK6. You can then scroll through the different molecules and see how well they align.

Cartesian Minimization of Docked Molecules

The structures you just looked at were the raw output from the virtual screen. It would be better, before making a conclusion on which of the 5,000 molecules bind the best, to run an energy minimization on each of those molecules and then rescore them based on any changes in their interactions with the receptor. Starting with the energy minimization, change directories to 007_cartesian_min. We will need to write an input and submit file for this simulation, so let's start with the input file and vi into cartmin.in. Write the following to that file: 

conformer_search_type                                        rigid
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
internal_energy_cutoff                                       100
ligand_atom_file                                             ../006_virtualscreenmpi/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                                     yes
continuous_score_secondary                                   no
cont_score_rec_filename                                      ../001_structure/6ME2.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                                        6ME2.virtual_screen.min
write_orientations                                           no
num_scored_conformers                                        1
rank_ligands                                                 no

Once you finish save and exit that file and then we can start to write the submit file. Vi into cartmin.sh and write the following to that file: 

#!/bin/bash
#SBATCH --time=48:00:00
#SBATCH --nodes=1
#SBATCH --ntasks=28
#SBATCH --job-name=mpi_6ME2_VS_AR
#SBATCH --output=6ME2_VS_mpi.out
#SBATCH -p long-28core

cd $SLURM_SUBMIT_DIR

mpirun -np 28 /gpfs/projects/AMS536/zzz.programs/dock6.9_release/bin/dock6.mpi -i cartmin.in -o cartmin.out

We are using mpi again to speed up the energy minimization, so if you choose to adjust the number of nodes being used to speed it up even further, adjust the file accordingly as mentioned in the virtual screen using MPI section. Once you are happy with the submit file, you can start it by typing the following command into the terminal: 

sbatch cartmin.sh

Once this finishes running, check the output files to make sure everything ran smoothly. Once that is done, we can now move to rescore these molecules and then visualize them.

Rescoring Docked Molecules

Now that we have minimized all the molecules docked in the virtual screen, we can rescore them based on their new conformations. In order to do this, change directories to 008_rescore. We will need to start by creating the input and submit files, so vi into rescore.in and write the following to that file: 

conformer_search_type                                        rigid
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
internal_energy_cutoff                                       100
ligand_atom_file                                             ../007_cartesian_min/6ME2.virtual_screen.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                           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                                      no
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/6ME2.lig.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/6ME2.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/6ME2.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_hms_score_ref_filename                            ../004_dock/6ME2.lig.min_scored.mol2
descriptor_hms_score_matching_coeff                          -5
descriptor_hms_score_rmsd_coeff                              1
descriptor_volume_score_reference_mol2_filename              ../004_dock/6ME2.lig.min_scored.mol2
descriptor_volume_score_overlap_compute_method               analytical
descriptor_weight_fps_score                                  1
descriptor_weight_pharmacophore_score                        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

Once you have finished, save and exit the file. Then vi rescore.sh and write the following to that file: 

#!/bin/sh
#SBATCH --partition=long-28core
#SBATCH --time=48:00:00
#SBATCH --nodes=1
#SBATCH --ntasks=28
#SBATCH --job-name=6ME2-rescore
#SBATCH --output=%x-%j.o

echo "starting DOCK6.9 Simulation"
module load intel/mpi/64/2018/18.0.3

mpirun -np 28 dock6.mpi -i rescore.in -o rescore.out
echo "done"

Again, we are using mpi so if you want to speed things up use for than one node and then adjust the submit file accordingly. Once you have finished, save and exit the file. Submit the job to the queue as described above. Once this has finished, check your input files to make sure everything ran smoothly. You will have generated multiple files named rescore.out.x, but every molecule described in the output files can also be found in the descriptor.output_scored.mol2 file.

Now, we can properly visualize the minimized and rescored molecules. Transfer the descriptor.output_scored.mol2 file to your computer and open Chimera. Load into Chimera the ligand and receptor files (6ME2.receptor_wH.mol2 and 6ME2.ligand_wH.mol2). Additionally, load in the the gridbox generated in showbox (6ME2.gridbox.pdb). Your screen should look as follows:

Receptor+ligand+showbox.png

Now, using DockView as described above: load in the descriptor.output_scored.mol2 file and begin moving through different ligands one at a time. Which ligands have the highest molecular weight, or the best score from DOCK? Using the ViewDock panel, go to (Column -> Show -> Molecular Weight) and a new column will appear in the panel which you can use to sort the molecules on any number of properties listed in that submenu. Once you have finished analyzing your different ligands, you have finished this part of the DOCK tutorial for 6ME2. Congratulations!

Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2022_DOCK_tutorial_1_with_PDBID_6ME2&oldid=15334"