Difference between revisions of "2019 DOCK tutorial 3 with PDBID 3JQZ"
From Rizzo_Lab
Stonybrook (talk | contribs) (→VI. Virtual Screen) |
(→Molecular Footprint) |
||
| (71 intermediate revisions by 3 users not shown) | |||
| Line 1: | Line 1: | ||
This tutorial contains, a step by step by approach to dock a known ligand to a known receptor. | This tutorial contains, a step by step by approach to dock a known ligand to a known receptor. | ||
| + | |||
| + | NOTE: JDB - Change all saved conformers to 50 for rigid/fad/flex and include a note on rank ordering the scores. | ||
= I. Introduction = | = I. Introduction = | ||
| Line 37: | Line 39: | ||
= II. Preparation of the ligand and receptor = | = II. Preparation of the ligand and receptor = | ||
| − | Download the | + | Download the PDB file 3JQZ from the PDB database save it in the '''000_files folder.''' |
===Checking the chosen structure=== | ===Checking the chosen structure=== | ||
| − | It can be helpful to read through the journal article associated with the PDB file to understand protonation states, charges, environmental conditions and other important information regarding the receptor and ligand under the | + | It can be helpful to read through the journal article associated with the PDB file to understand protonation states, charges, environmental conditions and other important information regarding the receptor and ligand under the crystallization conditions used to obtain the structure you are working with. |
To do this; | To do this; | ||
| Line 48: | Line 50: | ||
You can now explore the crystal structure. | You can now explore the crystal structure. | ||
| − | Hint: From the | + | Hint: From the toolbar at the top of the Chimera window, click Actions -> Focus, to put your structure back in the middle of the screen. |
Identify the main components of the model (receptor, ligand, solvent, surfactants, metal ions) | Identify the main components of the model (receptor, ligand, solvent, surfactants, metal ions) | ||
| Line 56: | Line 58: | ||
===Preparation of receptor=== | ===Preparation of receptor=== | ||
- Open the PDB file (3JQZ.pdb) using Chimera | - Open the PDB file (3JQZ.pdb) using Chimera | ||
| − | - Isolate the receptor using select tool and delete tool in Chimera. | + | - Isolate the receptor using the select tool and delete tool in Chimera. |
| − | - Save the isolated receptor as a mol2 file. ( | + | - Save the isolated receptor as a mol2 file. (3jqz_rec_noH.mol2) |
| − | - Open | + | - Open 3jqz_rec_noH.mol2 file again using Chimera and use the following instructions to prepare the receptor file to be used in DOCK. |
Tools -> Structure Editing -> Add H (To add Hydrogen atoms) | Tools -> Structure Editing -> Add H (To add Hydrogen atoms) | ||
Tools -> Structure Editing -> Add Charge (To add the charge use the latest AMBER force filed available for standard residues. Here we used AMBER ff14SB) | Tools -> Structure Editing -> Add Charge (To add the charge use the latest AMBER force filed available for standard residues. Here we used AMBER ff14SB) | ||
| − | Save as a mol2 file. ( | + | Save as a mol2 file. (3jqz_rec_withH.mol2) |
| − | - If you follow the step below all the above stated steps will automatically appear one after the other to prepare the receptor. | + | - If you follow the step below all the above-stated steps will automatically appear one after the other to prepare the receptor. |
Tools -> Structure/Binding Analysis -> DockPrep | Tools -> Structure/Binding Analysis -> DockPrep | ||
| Line 70: | Line 72: | ||
- Open the PDB file via Chimera. | - Open the PDB file via Chimera. | ||
- Using Chimera, isolate the ligand, add H atoms, add charge and save it as a mol2 file by following the same steps followed for the receptor. | - Using Chimera, isolate the ligand, add H atoms, add charge and save it as a mol2 file by following the same steps followed for the receptor. | ||
| + | - It is ''incredibly important'' both here and in all modeling experiments that you examine the ligand's protonation state and ring saturation. For this particular experiment, if you just take Chimera's protonation output you'll find it to be incorrect. Note that 6-membered ring structure will probably be protonated incorrectly - checking [https://www.sciencedirect.com/science/article/pii/S1047847710000912?via%3Dihub#aep-section-id28 the structure in the reference] shows that it should be aromatic, not saturated. You can fix this in Chimera using Tools > Structure Editing > Build Structure. | ||
| + | |||
| + | [[File:3jqz_good_bad_protonation.png|thumb|none|1000px|Left: The incorrect protonation of the ligand in 3JQZ. Right: The correct protonation.]] | ||
| + | |||
| − | Once all the files are prepared make sure to save the files in | + | Once all the files are prepared make sure to save the files in 001.dockprep folder. |
=III. Generating receptor surface and spheres= | =III. Generating receptor surface and spheres= | ||
| Line 85: | Line 91: | ||
===Generating spheres=== | ===Generating spheres=== | ||
| − | - Go to | + | - Go to 003_surface_spheres folder |
- Create a new input file to create spheres by typing vim INSPH and type the following lines inside the file. | - Create a new input file to create spheres by typing vim INSPH and type the following lines inside the file. | ||
| Line 117: | Line 123: | ||
This command will select all of the spheres within 10.0 angstroms of the ligand and output them to selected_spheres.sph. Visualize the selected spheres using Chimera to make sure the correct spheres are selected. Notice that, spheres around the ligand binding site are kept and all the other spheres are deleted in the image below. | This command will select all of the spheres within 10.0 angstroms of the ligand and output them to selected_spheres.sph. Visualize the selected spheres using Chimera to make sure the correct spheres are selected. Notice that, spheres around the ligand binding site are kept and all the other spheres are deleted in the image below. | ||
| − | [[File:Sphgen 3jqz bj3.png|thumb|center|1000px| | + | [[File:Sphgen 3jqz bj3.png|thumb|center|1000px|3jqz receptor and selected spheres]] |
=IV. Generating box and grid= | =IV. Generating box and grid= | ||
===Generating box=== | ===Generating box=== | ||
| − | Move to | + | '''Move to 003_gridbox directory''' |
| + | |||
Create a new file showbox.in and write the following lines in the file. | Create a new file showbox.in and write the following lines in the file. | ||
| Line 141: | Line 148: | ||
showbox < showbox.in | showbox < showbox.in | ||
| − | If this step is successful, you should see a new file (3jqz.box.pdb) in | + | If this step is successful, you should see a new file (3jqz.box.pdb) in 003_gridbox folder. |
[[File:Gridbox 3jqz bj1.png|thumb|center|1000px|Box grid (blue) generated for 3JQZ receptor (gray surface) at the ligand binding site (green surface). Selected spheres are shown in pink. ]] | [[File:Gridbox 3jqz bj1.png|thumb|center|1000px|Box grid (blue) generated for 3JQZ receptor (gray surface) at the ligand binding site (green surface). Selected spheres are shown in pink. ]] | ||
| Line 174: | Line 181: | ||
=V. Docking a single molecule for pose reproduction = | =V. Docking a single molecule for pose reproduction = | ||
| − | Under this section, the ligand for | + | Under this section, the ligand for 3jqz.pdb will be re-docked into the receptor. 3 Methods will be used to achieve this. |
1. rigid docking | 1. rigid docking | ||
| Line 199: | Line 206: | ||
read_mol_solvation no | read_mol_solvation no | ||
calculate_rmsd yes | calculate_rmsd yes | ||
| − | use_rmsd_reference_mol | + | use_rmsd_reference_mol yes |
| + | rmsd_reference_filename ../001.dockprep/3jqz_lig_withH.mol2 | ||
use_database_filter no | use_database_filter no | ||
orient_ligand no | orient_ligand no | ||
| Line 211: | Line 219: | ||
grid_score_vdw_scale 1 | grid_score_vdw_scale 1 | ||
grid_score_es_scale 1 | grid_score_es_scale 1 | ||
| − | grid_score_grid_prefix ../003. | + | grid_score_grid_prefix ../003.gridbox/grid |
multigrid_score_secondary no | multigrid_score_secondary no | ||
dock3.5_score_secondary no | dock3.5_score_secondary no | ||
| Line 245: | Line 253: | ||
If the process is successful a new file (3jqz.lig.min_scored.mol2) will be generated. You can compare how is it changed from the initial structure by analyzing the RMSD value generated in the file. Visualize the new mol2 file along with receptor and the initial ligand mol2 files using Chimera to see the differences. | If the process is successful a new file (3jqz.lig.min_scored.mol2) will be generated. You can compare how is it changed from the initial structure by analyzing the RMSD value generated in the file. Visualize the new mol2 file along with receptor and the initial ligand mol2 files using Chimera to see the differences. | ||
| − | [[File: | + | [[File:20190304 crystal and energy minimised ligand figure for wiki.jpg|thumb|center|1000px|3jqz with crystal pose ligand (blue) and the energy minimised ligand (magenta)]] |
==Rigid Docking== | ==Rigid Docking== | ||
| Line 263: | Line 271: | ||
read_mol_solvation no | read_mol_solvation no | ||
calculate_rmsd yes | calculate_rmsd yes | ||
| − | use_rmsd_reference_mol | + | use_rmsd_reference_mol yes |
| + | rmsd_reference_filename 3jqz.lig.min_scored.mol2 | ||
use_database_filter no | use_database_filter no | ||
orient_ligand yes | orient_ligand yes | ||
| Line 281: | Line 290: | ||
grid_score_vdw_scale 1 | grid_score_vdw_scale 1 | ||
grid_score_es_scale 1 | grid_score_es_scale 1 | ||
| − | grid_score_grid_prefix ../003. | + | grid_score_grid_prefix ../003.gridbox/grid |
multigrid_score_secondary no | multigrid_score_secondary no | ||
dock3.5_score_secondary no | dock3.5_score_secondary no | ||
| Line 353: | Line 362: | ||
internal_energy_rep_exp 12 | internal_energy_rep_exp 12 | ||
internal_energy_cutoff 100.0 | internal_energy_cutoff 100.0 | ||
| − | ligand_atom_file ../ | + | ligand_atom_file ../001.dockprep/3jqz_lig_withH.mol2 |
limit_max_ligands no | limit_max_ligands no | ||
skip_molecule no | skip_molecule no | ||
| Line 359: | Line 368: | ||
calculate_rmsd yes | calculate_rmsd yes | ||
use_rmsd_reference_mol yes | use_rmsd_reference_mol yes | ||
| − | rmsd_reference_filename ../ | + | rmsd_reference_filename ../001.dockprep/3jqz_lig_withH.mol2 |
use_database_filter no | use_database_filter no | ||
orient_ligand no | orient_ligand no | ||
| Line 441: | Line 450: | ||
read_mol_solvation no | read_mol_solvation no | ||
calculate_rmsd yes | calculate_rmsd yes | ||
| − | use_rmsd_reference_mol 3jqz.lig.min_scored.mol2 | + | use_rmsd_reference_mol yes |
| + | rmsd_reference_filename 3jqz.lig.min_scored.mol2 | ||
use_database_filter no | use_database_filter no | ||
orient_ligand yes | orient_ligand yes | ||
| Line 486: | Line 496: | ||
simplex_restraint_min no | simplex_restraint_min no | ||
atom_model all | atom_model all | ||
| − | vdw_defn_file | + | vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn |
| − | flex_defn_file | + | flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn |
| − | flex_drive_file | + | flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl |
ligand_outfile_prefix flex.out | ligand_outfile_prefix flex.out | ||
write_orientations no | write_orientations no | ||
| Line 512: | Line 522: | ||
conformer_search_type rigid | conformer_search_type rigid | ||
use_internal_energy no | use_internal_energy no | ||
| − | ligand_atom_file ../004.dock/ | + | ligand_atom_file ../004.dock/3jqz_lig_min_scored.mol2 |
limit_max_ligands no | limit_max_ligands no | ||
skip_molecule no | skip_molecule no | ||
| Line 534: | Line 544: | ||
footprint_similarity_score_secondary no | footprint_similarity_score_secondary no | ||
fps_score_use_footprint_reference_mol2 yes | fps_score_use_footprint_reference_mol2 yes | ||
| − | fps_score_footprint_reference_mol2_filename | + | fps_score_footprint_reference_mol2_filename 3jqz_lig_with.mol2 |
fps_score_foot_compare_type Euclidean | fps_score_foot_compare_type Euclidean | ||
fps_score_normalize_foot no | fps_score_normalize_foot no | ||
fps_score_foot_comp_all_residue yes | fps_score_foot_comp_all_residue yes | ||
| − | fps_score_receptor_filename ../1.dockprep/ | + | fps_score_receptor_filename ../1.dockprep/3jqz_rec_withH.mol2 |
fps_score_vdw_att_exp 6 | fps_score_vdw_att_exp 6 | ||
fps_score_vdw_rep_exp 12 | fps_score_vdw_rep_exp 12 | ||
| Line 558: | Line 568: | ||
flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn | flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn | ||
flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl | flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl | ||
| − | ligand_outfile_prefix | + | ligand_outfile_prefix 3jqz_footprint_min_cryst |
write_footprints yes | write_footprints yes | ||
write_hbonds yes | write_hbonds yes | ||
| Line 565: | Line 575: | ||
rank_ligands no | rank_ligands no | ||
| − | + | If successful, the following output files should be generated. (3jqz_footprint_min_cryst_footprint_scored.txt, 3jqz_footprint_min_cryst_hbond_scored.txt, 3jqz_footprint_min_cryst_scored.mol2) | |
| − | Use a python script to visualize the molecular footprint | + | Use [http://ringo.ams.sunysb.edu/~rizzo/StonyBrook/downloads/plot_footprint_single_magnitude.py a python script] to visualize the molecular footprint (read the txt files to pdf). Once the script has run, a pdf will be generated and the molecular footprint should be similar to the image below. Notice the deviations in energy at charged amino acid residues. Those residues contribute more electrostatically towards the interaction between the ligand and the receptor and therefore can be identified as important towards the binding of new ligands. In this case, the crystallographic ligand lidocaine has a +1 net charge. |
| − | [[File: | + | [[File:Correctly named footprint for wiki figure.jpg|thumb|center|1000px| Molecular Footprint of crystal pose ligand (blue) and the energy minimised ligand (red)]] |
= VI. Virtual Screen = | = VI. Virtual Screen = | ||
Virtual screening is the method of screening a ligand library (drug-like molecules) to filter the best ligands which can bind to the binding site of a specific receptor. Here we will be using a ligand library which contains 25000 molecules to select the best ligands which can be used to replace the original ligand that the PDB file contained. | Virtual screening is the method of screening a ligand library (drug-like molecules) to filter the best ligands which can bind to the binding site of a specific receptor. Here we will be using a ligand library which contains 25000 molecules to select the best ligands which can be used to replace the original ligand that the PDB file contained. | ||
| − | Move to the directory | + | '''Move to the directory 007_virtual_screen.''' |
| + | |||
| + | Copy the ligand library to the same directory. | ||
| + | Copy '''VS_library_25k.mol2''' from the zzz.programs directory in project space. The VS_library_25k.mol2 file is an in house library of small molecules curated from ZINC. | ||
Create a new input file for virtual screen. | Create a new input file for virtual screen. | ||
touch virtual.in | touch virtual.in | ||
Use the input file to perform virtual screen using DOCK6. | Use the input file to perform virtual screen using DOCK6. | ||
| − | dock6 -i virtual.in | + | dock6 -i virtual.in |
Use the following lines to answer the prompted questions. | Use the following lines to answer the prompted questions. | ||
| Line 597: | Line 610: | ||
internal_energy_rep_exp 12 | internal_energy_rep_exp 12 | ||
internal_energy_cutoff 100.0 | internal_energy_cutoff 100.0 | ||
| − | ligand_atom_file | + | ligand_atom_file VS_library_25K.mol2 |
limit_max_ligands no | limit_max_ligands no | ||
skip_molecule no | skip_molecule no | ||
read_mol_solvation no | read_mol_solvation no | ||
| − | calculate_rmsd | + | calculate_rmsd no |
| − | |||
| − | |||
use_database_filter no | use_database_filter no | ||
orient_ligand yes | orient_ligand yes | ||
automated_matching yes | automated_matching yes | ||
| − | receptor_site_file ../ | + | receptor_site_file ../002_surface_spheres/selected_spheres.sph |
max_orientations 1000 | max_orientations 1000 | ||
critical_points no | critical_points no | ||
| Line 621: | Line 632: | ||
grid_score_vdw_scale 1 | grid_score_vdw_scale 1 | ||
grid_score_es_scale 1 | grid_score_es_scale 1 | ||
| − | grid_score_grid_prefix ../ | + | grid_score_grid_prefix ../003_gridbox/grid |
multigrid_score_secondary no | multigrid_score_secondary no | ||
dock3.5_score_secondary no | dock3.5_score_secondary no | ||
| Line 656: | Line 667: | ||
rank_ligands no | rank_ligands no | ||
| − | Since the ligand library contains 25000 molecules, it is going to take a long time to complete the virtual screen. Therefore we are going to use the mpi version of DOCK to complete the virtual screen. | + | Since the ligand library contains 25000 molecules, it is going to take a long time to complete the virtual screen. Therefore we are going to use the mpi version of DOCK to complete the virtual screen. |
| + | |||
| + | You can terminate the virtual screen which is already running by pressing ctrl+c and follow the next steps to submit your job to the queue. | ||
=VII.Virtual Screen (MPI)= | =VII.Virtual Screen (MPI)= | ||
Until now we used the head node of the seawulf cluster. By using the mpi version of DOCK we will be using 4 processors that contain 28 nodes in each and it will complete the virtual screen quicker. | Until now we used the head node of the seawulf cluster. By using the mpi version of DOCK we will be using 4 processors that contain 28 nodes in each and it will complete the virtual screen quicker. | ||
| − | Move to a new directory | + | |
| − | Copy the input file and the ligand database file from | + | '''Move to a new directory (008.virtual_screen_mpi).''' |
| + | |||
| + | Copy the input file and the ligand database file from 007.virtual_screen directory. | ||
To submit the job to the seawulf cluster we are using a new file. (virtual.sh) | To submit the job to the seawulf cluster we are using a new file. (virtual.sh) | ||
| Line 670: | Line 685: | ||
#PBS -l nodes=4:ppn=28 | #PBS -l nodes=4:ppn=28 | ||
#PBS -q long | #PBS -q long | ||
| − | #PBS -N | + | #PBS -N 3jqz.virtual |
#PBS -V | #PBS -V | ||
cd $PBS_O_WORKDIR | cd $PBS_O_WORKDIR | ||
| − | mpirun -np 112 dock6.mpi -i virtual.in -o | + | mpirun -np 112 dock6.mpi -i virtual.in -o 3jqz.virtual.mpi.out |
The virtual screen job can be submitted to the cluster using the following command. | The virtual screen job can be submitted to the cluster using the following command. | ||
| Line 685: | Line 700: | ||
Here we will use the docked molecules and perform a cartesian minimization of them. | Here we will use the docked molecules and perform a cartesian minimization of them. | ||
| − | Move to the directory | + | '''Move to the directory 009.cartesianmin.''' |
Create a new input file for the minimization. | Create a new input file for the minimization. | ||
| Line 698: | Line 713: | ||
internal_energy_rep_exp 12 | internal_energy_rep_exp 12 | ||
internal_energy_cutoff 100.0 | internal_energy_cutoff 100.0 | ||
| − | ligand_atom_file | + | ligand_atom_file ../008.virtual_screen_mpi/3jqz.virtualscreen_scored.mol2 |
limit_max_ligands no | limit_max_ligands no | ||
skip_molecule no | skip_molecule no | ||
| Line 717: | Line 732: | ||
continuous_score_primary yes | continuous_score_primary yes | ||
continuous_score_secondary no | continuous_score_secondary no | ||
| − | cont_score_rec_filename ../ | + | cont_score_rec_filename ../001.dockprep/3jqz.rec.charged.mol2 SHOULD WE GENERATE A CHARGED STRUCTURE?? |
cont_score_att_exp 6 | cont_score_att_exp 6 | ||
cont_score_rep_exp 12 | cont_score_rep_exp 12 | ||
| Line 744: | Line 759: | ||
simplex_restraint_min no | simplex_restraint_min no | ||
atom_model all | atom_model all | ||
| − | vdw_defn_file | + | vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn |
| − | flex_defn_file | + | flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn |
| − | flex_drive_file | + | flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl |
| − | ligand_outfile_prefix | + | ligand_outfile_prefix 3jqz.virtualscreen.min |
write_orientations no | write_orientations no | ||
num_scored_conformers 1 | num_scored_conformers 1 | ||
rank_ligands no | rank_ligands no | ||
| + | |||
| + | You can terminate the cartesian minimization which is already running by pressing ctrl+c and follow the next steps to submit your job to the queue. | ||
Create a new submission script for the minimization. | Create a new submission script for the minimization. | ||
| − | + | vi min.sh | |
#!/bin/bash | #!/bin/bash | ||
| Line 763: | Line 780: | ||
cd $PBS_O_WORKDIR | cd $PBS_O_WORKDIR | ||
dock6 -i min.in -o min.out | dock6 -i min.in -o min.out | ||
| + | |||
| + | |||
| + | The minimization can be submitted to the cluster using the following command. | ||
| + | qsub min.sh | ||
| + | |||
| + | Type the following command to check the status of the job and other jobs that you submitted via your user login. | ||
| + | qstat -u username | ||
=IX.Rescoring Docked Molecules= | =IX.Rescoring Docked Molecules= | ||
We would also like to rank our docked ligands and extract the 100 best ligands (which have the most negative, lowest scores). Footprint similarity, pharmacophore score, tanimoto score, the hungarian and the volume overlap score will all be used by Dock to rescore the virtual screen. | We would also like to rank our docked ligands and extract the 100 best ligands (which have the most negative, lowest scores). Footprint similarity, pharmacophore score, tanimoto score, the hungarian and the volume overlap score will all be used by Dock to rescore the virtual screen. | ||
| − | Move to the directory | + | '''Move to the directory 010.rescore.''' |
Create a new input file for the rescoring. | Create a new input file for the rescoring. | ||
| Line 778: | Line 802: | ||
internal_energy_rep_exp 12 | internal_energy_rep_exp 12 | ||
internal_energy_cutoff 100.0 | internal_energy_cutoff 100.0 | ||
| − | ligand_atom_file ../ | + | ligand_atom_file ../009.cartesianmin/3jqz.virtualscreen.min_scored.mol2 |
limit_max_ligands no | limit_max_ligands no | ||
skip_molecule no | skip_molecule no | ||
| Line 784: | Line 808: | ||
calculate_rmsd no | calculate_rmsd no | ||
use_database_filter no | use_database_filter no | ||
| − | orient_ligand no | + | '''orient_ligand no''' |
bump_filter no | bump_filter no | ||
score_molecules yes | score_molecules yes | ||
contact_score_primary no | contact_score_primary no | ||
contact_score_secondary no | contact_score_secondary no | ||
| − | grid_score_primary no | + | '''grid_score_primary no''' |
grid_score_secondary no | grid_score_secondary no | ||
multigrid_score_primary no | multigrid_score_primary no | ||
| Line 801: | Line 825: | ||
pharmacophore_score_primary no | pharmacophore_score_primary no | ||
pharmacophore_score_secondary no | pharmacophore_score_secondary no | ||
| − | descriptor_score_primary yes | + | '''descriptor_score_primary yes''' |
descriptor_score_secondary no | descriptor_score_secondary no | ||
| − | descriptor_use_grid_score no | + | '''descriptor_use_grid_score no''' |
| − | descriptor_use_multigrid_score no | + | '''descriptor_use_multigrid_score no''' |
| − | descriptor_use_continuous_energy no | + | '''descriptor_use_continuous_energy no''' |
descriptor_use_footprint_similarity yes | descriptor_use_footprint_similarity yes | ||
| − | descriptor_use_pharmacophore_score yes | + | '''descriptor_use_pharmacophore_score yes''' |
| − | descriptor_use_tanimoto yes | + | '''descriptor_use_tanimoto yes''' |
| − | descriptor_use_hungarian yes | + | '''descriptor_use_hungarian yes''' |
| − | descriptor_use_volume_overlap yes | + | '''descriptor_use_volume_overlap yes''' |
| − | descriptor_fps_use_footprint_reference_mol2 yes | + | '''descriptor_fps_use_footprint_reference_mol2 yes''' |
| − | descriptor_fps_footprint_reference_mol2_filename ../ | + | descriptor_fps_footprint_reference_mol2_filename ../004.dock/3jqz.lig.min_scored.mol2 |
descriptor_fps_foot_compare_type Euclidean | descriptor_fps_foot_compare_type Euclidean | ||
descriptor_fps_normalize_foot no | descriptor_fps_normalize_foot no | ||
descriptor_fps_foot_comp_all_residue yes | descriptor_fps_foot_comp_all_residue yes | ||
| − | descriptor_fps_receptor_filename ../ | + | '''descriptor_fps_receptor_filename ../001.dockprep/3jqz.rec.withH.charged.mol2''' |
descriptor_fps_vdw_att_exp 6 | descriptor_fps_vdw_att_exp 6 | ||
descriptor_fps_vdw_rep_exp 12 | descriptor_fps_vdw_rep_exp 12 | ||
| Line 825: | Line 849: | ||
descriptor_fps_es_fp_scale 1 | descriptor_fps_es_fp_scale 1 | ||
descriptor_fps_hb_fp_scale 0 | descriptor_fps_hb_fp_scale 0 | ||
| − | descriptor_fms_score_use_ref_mol2 yes | + | '''descriptor_fms_score_use_ref_mol2 yes''' |
| − | descriptor_fms_score_ref_mol2_filename ../ | + | descriptor_fms_score_ref_mol2_filename ../004.dock/3jqz.lig.min_scored.mol2 |
descriptor_fms_score_write_reference_pharmacophore_mol2 no | descriptor_fms_score_write_reference_pharmacophore_mol2 no | ||
descriptor_fms_score_write_reference_pharmacophore_txt no | descriptor_fms_score_write_reference_pharmacophore_txt no | ||
| Line 837: | Line 861: | ||
descriptor_fms_score_match_proj_cutoff 0.7071 | descriptor_fms_score_match_proj_cutoff 0.7071 | ||
descriptor_fms_score_max_score 20 | descriptor_fms_score_max_score 20 | ||
| − | descriptor_fingerprint_ref_filename ../ | + | descriptor_fingerprint_ref_filename ../004.dock/3jqz.lig.min_scored.mol2 |
| − | descriptor_hungarian_ref_filename ../ | + | descriptor_hungarian_ref_filename ../004.dock/3jqz.lig.min_scored.mol2 |
descriptor_hungarian_matching_coeff -5 | descriptor_hungarian_matching_coeff -5 | ||
descriptor_hungarian_rmsd_coeff 1 | descriptor_hungarian_rmsd_coeff 1 | ||
| − | descriptor_volume_reference_mol2_filename ../ | + | descriptor_volume_reference_mol2_filename ../004.dock/3jqz.lig.min_scored.mol2 |
descriptor_volume_overlap_compute_method analytical | descriptor_volume_overlap_compute_method analytical | ||
descriptor_weight_fps_score 1 | descriptor_weight_fps_score 1 | ||
| Line 852: | Line 876: | ||
SASA_score_secondary no | SASA_score_secondary no | ||
amber_score_secondary no | amber_score_secondary no | ||
| − | minimize_ligand no | + | '''minimize_ligand no''' |
atom_model all | atom_model all | ||
vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn | vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn | ||
| Line 859: | Line 883: | ||
chem_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/chem.defn | chem_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/chem.defn | ||
pharmacophore_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/ph4.defn | pharmacophore_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/ph4.defn | ||
| − | ligand_outfile_prefix descriptor. | + | ligand_outfile_prefix descriptor.out |
| − | write_footprints yes | + | '''write_footprints yes''' |
| − | write_hbonds yes | + | '''write_hbonds yes''' |
write_orientations no | write_orientations no | ||
num_scored_conformers 1 | num_scored_conformers 1 | ||
rank_ligands no | rank_ligands no | ||
Latest revision as of 11:39, 7 March 2022
This tutorial contains, a step by step by approach to dock a known ligand to a known receptor.
NOTE: JDB - Change all saved conformers to 50 for rigid/fad/flex and include a note on rank ordering the scores.
Contents
- 1 I. Introduction
- 2 II. Preparation of the ligand and receptor
- 3 III. Generating receptor surface and spheres
- 4 IV. Generating box and grid
- 5 V. Docking a single molecule for pose reproduction
- 6 VI. Virtual Screen
- 7 VII.Virtual Screen (MPI)
- 8 VIII.Cartesian Minimization
- 9 IX.Rescoring Docked Molecules
I. Introduction
DOCK
DOCK is a molecular docking program used in drug discovery. It was developed by Irwin D. Kuntz, Jr. and colleagues at UCSF (see UCSF DOCK). This program, given a protein binding site and a small molecule, tries to predict the correct binding mode of the small molecule in the binding site, and the associated binding energy. Small molecules with highly favorable binding energies could be new drug leads. This makes DOCK a valuable drug discovery tool. DOCK is typically used to screen massive libraries of millions of compounds against a protein to isolate potential drug leads. These leads are then further studied, and could eventually result in a new, marketable drug. DOCK works well as a screening procedure for generating leads, but is not currently as useful for optimization of those leads.
DOCK 6 uses an incremental construction algorithm called anchor and grow. It is described by a three-step process:
- Rigid portion of ligand (anchor) is docked by geometric methods.
- Non-rigid segments added in layers; energy minimized.
- The resulting configurations are 'pruned' and energy re-minimized, yielding the docked configurations.
If you are interested in reading more about DOCK 6 and the background to these processes, interested readers may reference the DOCK User's Manual.
3JQZ
This tutorial is based on the PDB file 3JQZ downloaded from the PDB Database. 3JQZ is the 3.3 Angstrom resolution crystal structure of Human serum albumin complexed with Lidocaine.
Hint: for docking practice you should identify a high resolution crystal structure with a well defined ligand binding site, not situated at the interface of two chains and with minimal active site waters.
Organization of Directories
Maintaining a clearly organised set of folders will be helpful in finding specific files, calling files and most importantly keeping track of everything you do. We recommend maintaining the following directory structure throughout the tutorial.
000_files
001_dockprep
002_surface_spheres
003_gridbox
004_dock
005
006_footprint
007_virtual_screen
008_virtual_screen_mpi
009_cartesianmin
0010_rescore
II. Preparation of the ligand and receptor
Download the PDB file 3JQZ from the PDB database save it in the 000_files folder.
Checking the chosen structure
It can be helpful to read through the journal article associated with the PDB file to understand protonation states, charges, environmental conditions and other important information regarding the receptor and ligand under the crystallization conditions used to obtain the structure you are working with.
To do this;
- Open Chimera - File -> Fetch by ID -> Type 3JQZ into the 'PDB' box -> Fetch
You can now explore the crystal structure.
Hint: From the toolbar at the top of the Chimera window, click Actions -> Focus, to put your structure back in the middle of the screen.
Identify the main components of the model (receptor, ligand, solvent, surfactants, metal ions)
- Carefully look to identify if there are any missing residues or missing loops, this can be the case where regions of the protein were flexible. For the purposes of learning the software, we want to steer clear of any systems with missing loops. The PDB 3JQZ, chosen for this tutorial, does not contain any missing loops or residues.
Preparation of receptor
- Open the PDB file (3JQZ.pdb) using Chimera - Isolate the receptor using the select tool and delete tool in Chimera. - Save the isolated receptor as a mol2 file. (3jqz_rec_noH.mol2)
- Open 3jqz_rec_noH.mol2 file again using Chimera and use the following instructions to prepare the receptor file to be used in DOCK.
Tools -> Structure Editing -> Add H (To add Hydrogen atoms)
Tools -> Structure Editing -> Add Charge (To add the charge use the latest AMBER force filed available for standard residues. Here we used AMBER ff14SB)
Save as a mol2 file. (3jqz_rec_withH.mol2)
- If you follow the step below all the above-stated steps will automatically appear one after the other to prepare the receptor.
Tools -> Structure/Binding Analysis -> DockPrep
Preparation of ligand
- Open the PDB file via Chimera. - Using Chimera, isolate the ligand, add H atoms, add charge and save it as a mol2 file by following the same steps followed for the receptor. - It is incredibly important both here and in all modeling experiments that you examine the ligand's protonation state and ring saturation. For this particular experiment, if you just take Chimera's protonation output you'll find it to be incorrect. Note that 6-membered ring structure will probably be protonated incorrectly - checking the structure in the reference shows that it should be aromatic, not saturated. You can fix this in Chimera using Tools > Structure Editing > Build Structure.
Once all the files are prepared make sure to save the files in 001.dockprep folder.
III. Generating receptor surface and spheres
Preparation of DMS file
- Open 3jqz_rec_noH.mol2 using chimera. - Action -> Surface -> Show - Tools -> Structure Editing -> Write DMS - Save the 3jqz_rec_noH.dms into 3.surface_spheres folder
Reopen the file and make sure the surface was generated.
Transfer all the folders created so far to Seawulf Cluster to be used in DOCK.
Generating spheres
- Go to 003_surface_spheres folder - Create a new input file to create spheres by typing vim INSPH and type the following lines inside the file.
3jqz_rec_noH.dms R X 0.0 4.0 1.4 3jqz_rec.sph
The first line 3jqz_rec_noH.dms specifies the input file. R indicates that spheres generated will be outside of the receptor surface. X specifies all the points will be used. 0.0 is the distance in angstroms and it will avoid steric clashes. 4.0 is the maximum surface radius of the spheres and 1.4 is the minimum radius in angstroms. The last line 3jqz_spheres.sph creates the sph file that contains clustered spheres.
Once the INSPH file is ready, type the following command to generate the spheres.
sphgen -i INSPH -o OUTSPH
Once sphgen command is successful, 3jqz_spheres.sph file will be created. Open it up using Chimera along with 3jqz_rec_noH.mol2 file. You should get a similar output like the image below.
Some potential problems:
- If you are unable to run the sphgen command because the command isn't found, ensure that DOCK is in your PATH.
- If you run sphgen multiple times, be sure to delete any files generated by sphgen before running it again (ex. OUTSPH), as it isn't done automatically. You should just have '3jqz_rec_noH.dms' and 'INSPH' before running the command.
Selecting Spheres
Here we will be selecting the spheres which define the binding pocket of the ligand because we are trying to direct the ligand towards that binding site rather than all over the receptor. To select the spheres execute the following command.
sphere_selector 3jqz_rec.sph ../001_dockprep/3jqz_lig_withH.mol2 10.0
This command will select all of the spheres within 10.0 angstroms of the ligand and output them to selected_spheres.sph. Visualize the selected spheres using Chimera to make sure the correct spheres are selected. Notice that, spheres around the ligand binding site are kept and all the other spheres are deleted in the image below.
IV. Generating box and grid
Generating box
Move to 003_gridbox directory
Create a new file showbox.in and write the following lines in the file.
Y 8.0 ../002_surface_spheres/selected_spheres.sph 1 3jqz.box.pdb
Each of the above lines indicates that;
We intend to generate a box The box length should be 8 Angstroms Use the selected_spheres file in the designated location The name of the file that contains the generated box.
Use the following command to generate the box.
showbox < showbox.in
If this step is successful, you should see a new file (3jqz.box.pdb) in 003_gridbox folder.
Generating grid
Create a new file (grid.in)
Include the following lines in the grid.in file:
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_dockprep/3jqz_rec_withH.mol2 box_file 3jqz.box.pdb vdw_definition_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn score_grid_prefix grid
Use the following command to generate the grid.
grid -i grid.in -o gridinfo.out
If the command is successful, three new files will be generated. (gridinfo.out, grid.nrg, grid.bmp). Go through gridinfo.out file to make sure all the information about the receptor in the file matches with the original information of the receptor. (Eg:- Total charge, residues and their charges) If the information doesn't match, that means you have made an error in one of the steps that you followed so far.
V. Docking a single molecule for pose reproduction
Under this section, the ligand for 3jqz.pdb will be re-docked into the receptor. 3 Methods will be used to achieve this.
1. rigid docking
2. fixed anchor docking
3. flexible docking
Energy minimization
Before performing docking, here the ligand will be subjected to energy minimization in order to remove unfavorable clashes. These clashes will affect rigid docking because in rigid docking the ligand will be docked as the complete ligand, whereas in other docking methods the ligand will be broken into fragments and the ligand will be built step by step considering favorable orientations and torsion angles after each fragment addition.
Go to the directory 4.dock and a create a new file (touch min.in) and enter the command below.
dock6 -i min.in
Answer the prompted questions using the answers given below or include the following lines in the min.in file at before entering the above command to avoid answering the questions manually.
conformer_search_type rigid use_internal_energy yes internal_energy_rep_exp 12 internal_energy_cutoff 100.0 ligand_atom_file ../001.dockprep/3jqz_lig_withH.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd yes use_rmsd_reference_mol yes rmsd_reference_filename ../001.dockprep/3jqz_lig_withH.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_restraint_min yes simplex_coefficient_restraint 10.0 atom_model all vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix 3jqz.lig.min write_orientations no num_scored_conformers 1 rank_ligands no
If the process is successful a new file (3jqz.lig.min_scored.mol2) will be generated. You can compare how is it changed from the initial structure by analyzing the RMSD value generated in the file. Visualize the new mol2 file along with receptor and the initial ligand mol2 files using Chimera to see the differences.
Rigid Docking
Create an input file for rigid docking
touch rigid.in
Run dock using the created input file.
dock6 -i rigid.in -o rigid.out
Follow a similar approach as we did for minimization to answer the prompted questions by either answering them manually using the answers in the lines below or by including the following lines in the input file before running dock.
conformer_search_type rigid use_internal_energy yes ligand_atom_file 3jqz.lig.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 3jqz.lig.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 simplex_max_iterations 1000 simplex_tors_premin_iterations 0 simplex_max_cycles 1 simplex_score_converge 0.1 simplex_cycle_converge 1.0 simplex_trans_step 1.0 simplex_rot_step 0.1 simplex_tors_step 10.0 simplex_random_seed 0 simplex_restraint_min no atom_model all vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix rigid.out write_orientations no num_scored_conformers 1 rank_ligands no
Once rigid docking is successful, you will get an output file. (rigid.out_scored.mol2) Visualize the output file using Chimera by following steps to check the rigid docking success.
Open Chimera File -> Open -> 3jqz_rec_withH.mol2 File -> Open -> 3jqz_lig_withH.mol2 Tools -> Surface/binding Analysis -> ViewDock -> Select the Rigid Dock output file. (004.dock/rigid.out_scored.mol2)
Under File type, select Dock4,5 or 6
Once everything is loaded go to the ViewDock window and use it's menu to view all the calculated properties regarding the rigid docked ligand by following the steps below.
Column -> Show -> gridscore Column -> Show -> HA_RMSDs Follow the same steps to get all the properties
Your visualized structure should be similar to the image below.
Fixed Anchor Docking
Create an input file for fixed anchor docking.
touch fixed.in
Use the input file to perform fixed anchor docking
dock6 -i fixed.in -o fixed.out
Use the following lines to answer the prompted questions as we did in rigid docking.
conformer_search_type flex user_specified_anchor no limit_max_anchors no min_anchor_size 5 pruning_use_clustering yes pruning_max_orients 1000 pruning_clustering_cutoff 100 pruning_conformer_score_cutoff 100.0 pruning_conformer_score_scaling_factor 1.0 use_clash_overlap no write_growth_tree no write_fragment_libraries no use_internal_energy yes internal_energy_rep_exp 12 internal_energy_cutoff 100.0 ligand_atom_file ../001.dockprep/3jqz_lig_withH.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd yes use_rmsd_reference_mol yes rmsd_reference_filename ../001.dockprep/3jqz_lig_withH.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 ../002.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/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix fixed.out write_orientations no num_scored_conformers 100 write_conformations no cluster_conformations yes cluster_rmsd_threshold 2.0 rank_ligands no
Once docking is completed an output file will be generated. (fixed.out_scored.mol2) Follow the same method used in rigid docking to visualize the docked poses using Chimera. Once it is visualized, it should like the image below.
Flexible Docking
Create a new input file for flexible docking. (flex.in)
touch flex.in
Use the created input file to perform flexible docking using DOCK6.
dock6 -i flex.in -o flex.out
Answer the prompted questions using the following lines as we did in rigid and fixed anchor docking.
conformer_search_type flex user_specified_anchor no limit_max_anchors no min_anchor_size 5 pruning_use_clustering yes pruning_max_orients 1000 pruning_clustering_cutoff 100 pruning_conformer_score_cutoff 100.0 pruning_conformer_score_scaling_factor 1.0 use_clash_overlap no write_growth_tree no write_fragment_libraries no use_internal_energy yes internal_energy_rep_exp 12 internal_energy_cutoff 100.0 ligand_atom_file 3jqz.lig.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 3jqz.lig.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/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix flex.out write_orientations no num_scored_conformers 1 rank_ligands no
Once flexible docking is completed an output mol2 file will be generated. (flex.out_scored.mol2). Use the visualization steps used in rigid and fixed anchor docking and study the properties of the docking results.
Molecular Footprint
Molecular footprints can be used to determine how a ligand interacts with the receptor. Usually, the molecular footprint shows electrostatic interactions and Van der Waals interactions. Here, the molecular footprint will be used to determine how the ligand interacts with the receptor before and after minimization. To generate molecular footprints use following steps.
Go to directory 006.footprint
Generate an input file by typing;
touch footprint.in
Use DOCK6 to generate footprints
dock6 -i footprint.in
Use the following lines to answer the prompted questions.
conformer_search_type rigid use_internal_energy no ligand_atom_file ../004.dock/3jqz_lig_min_scored.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd no use_database_filter no orient_ligand no bump_filter no score_molecules yes contact_score_primary no contact_score_secondary no grid_score_primary no grid_score_secondary no multigrid_score_primary no multigrid_score_secondary no dock3.5_score_primary no dock3.5_score_secondary no continuous_score_primary no continuous_score_secondary no footprint_similarity_score_primary yes footprint_similarity_score_secondary no fps_score_use_footprint_reference_mol2 yes fps_score_footprint_reference_mol2_filename 3jqz_lig_with.mol2 fps_score_foot_compare_type Euclidean fps_score_normalize_foot no fps_score_foot_comp_all_residue yes fps_score_receptor_filename ../1.dockprep/3jqz_rec_withH.mol2 fps_score_vdw_att_exp 6 fps_score_vdw_rep_exp 12 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/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix 3jqz_footprint_min_cryst write_footprints yes write_hbonds yes write_orientations no num_scored_conformers 1 rank_ligands no
If successful, the following output files should be generated. (3jqz_footprint_min_cryst_footprint_scored.txt, 3jqz_footprint_min_cryst_hbond_scored.txt, 3jqz_footprint_min_cryst_scored.mol2)
Use a python script to visualize the molecular footprint (read the txt files to pdf). Once the script has run, a pdf will be generated and the molecular footprint should be similar to the image below. Notice the deviations in energy at charged amino acid residues. Those residues contribute more electrostatically towards the interaction between the ligand and the receptor and therefore can be identified as important towards the binding of new ligands. In this case, the crystallographic ligand lidocaine has a +1 net charge.
VI. Virtual Screen
Virtual screening is the method of screening a ligand library (drug-like molecules) to filter the best ligands which can bind to the binding site of a specific receptor. Here we will be using a ligand library which contains 25000 molecules to select the best ligands which can be used to replace the original ligand that the PDB file contained.
Move to the directory 007_virtual_screen.
Copy the ligand library to the same directory. Copy VS_library_25k.mol2 from the zzz.programs directory in project space. The VS_library_25k.mol2 file is an in house library of small molecules curated from ZINC.
Create a new input file for virtual screen.
touch virtual.in
Use the input file to perform virtual screen using DOCK6.
dock6 -i virtual.in
Use the following lines to answer the prompted questions.
conformer_search_type flex user_specified_anchor no limit_max_anchors no min_anchor_size 5 pruning_use_clustering yes pruning_max_orients 1000 pruning_clustering_cutoff 100 pruning_conformer_score_cutoff 100.0 pruning_conformer_score_scaling_factor 1.0 use_clash_overlap no write_growth_tree no write_fragment_libraries no use_internal_energy yes internal_energy_rep_exp 12 internal_energy_cutoff 100.0 ligand_atom_file VS_library_25K.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/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix virtual.out write_orientations no num_scored_conformers 1 rank_ligands no
Since the ligand library contains 25000 molecules, it is going to take a long time to complete the virtual screen. Therefore we are going to use the mpi version of DOCK to complete the virtual screen.
You can terminate the virtual screen which is already running by pressing ctrl+c and follow the next steps to submit your job to the queue.
VII.Virtual Screen (MPI)
Until now we used the head node of the seawulf cluster. By using the mpi version of DOCK we will be using 4 processors that contain 28 nodes in each and it will complete the virtual screen quicker.
Move to a new directory (008.virtual_screen_mpi).
Copy the input file and the ligand database file from 007.virtual_screen directory.
To submit the job to the seawulf cluster we are using a new file. (virtual.sh)
vim virtual.sh
Type the following lines in the new file.
#!/bin/bash #PBS -l walltime=48:00:00 #PBS -l nodes=4:ppn=28 #PBS -q long #PBS -N 3jqz.virtual #PBS -V cd $PBS_O_WORKDIR mpirun -np 112 dock6.mpi -i virtual.in -o 3jqz.virtual.mpi.out
The virtual screen job can be submitted to the cluster using the following command.
qsub virtual.sh
Type the following command to check the status of the job and other jobs that you submitted via your user login.
qstat -u username
VIII.Cartesian Minimization
Here we will use the docked molecules and perform a cartesian minimization of them.
Move to the directory 009.cartesianmin.
Create a new input file for the minimization.
touch min.in
Use the input file to perform the cartesian minimization using DOCK6.
dock6 -i min.in
Use the following lines to answer the prompted questions.
conformer_search_type rigid use_internal_energy yes internal_energy_rep_exp 12 internal_energy_cutoff 100.0 ligand_atom_file ../008.virtual_screen_mpi/3jqz.virtualscreen_scored.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd no use_database_filter no orient_ligand no bump_filter no score_molecules yes contact_score_primary no contact_score_secondary no grid_score_primary no grid_score_secondary no multigrid_score_primary no multigrid_score_secondary no dock3.5_score_primary no dock3.5_score_secondary no continuous_score_primary yes continuous_score_secondary no cont_score_rec_filename ../001.dockprep/3jqz.rec.charged.mol2 SHOULD WE GENERATE A CHARGED STRUCTURE?? cont_score_att_exp 6 cont_score_rep_exp 12 cont_score_rep_rad_scale 1 cont_score_use_dist_dep_dielectric yes cont_score_dielectric 4.0 cont_score_vdw_scale 1 cont_score_es_scale 1 footprint_similarity_score_secondary no pharmacophore_score_secondary no descriptor_score_secondary no gbsa_zou_score_secondary no gbsa_hawkins_score_secondary no SASA_score_secondary no amber_score_secondary no minimize_ligand yes simplex_max_iterations 1000 simplex_tors_premin_iterations 0 simplex_max_cycles 1 simplex_score_converge 0.1 simplex_cycle_converge 1.0 simplex_trans_step 1.0 simplex_rot_step 0.1 simplex_tors_step 10.0 simplex_random_seed 0 simplex_restraint_min no atom_model all vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix 3jqz.virtualscreen.min write_orientations no num_scored_conformers 1 rank_ligands no
You can terminate the cartesian minimization which is already running by pressing ctrl+c and follow the next steps to submit your job to the queue.
Create a new submission script for the minimization.
vi min.sh
#!/bin/bash #PBS -l walltime=48:00:00 #PBS -l nodes=1:ppn=28 #PBS -q long #PBS -N cartesian min #PBS -V cd $PBS_O_WORKDIR dock6 -i min.in -o min.out
The minimization can be submitted to the cluster using the following command.
qsub min.sh
Type the following command to check the status of the job and other jobs that you submitted via your user login.
qstat -u username
IX.Rescoring Docked Molecules
We would also like to rank our docked ligands and extract the 100 best ligands (which have the most negative, lowest scores). Footprint similarity, pharmacophore score, tanimoto score, the hungarian and the volume overlap score will all be used by Dock to rescore the virtual screen.
Move to the directory 010.rescore.
Create a new input file for the rescoring.
touch rescore.in
Use the input file to perform the cartesian minimization using DOCK6.
dock6 -i rescore.in
conformer_search_type rigid use_internal_energy yes internal_energy_rep_exp 12 internal_energy_cutoff 100.0 ligand_atom_file ../009.cartesianmin/3jqz.virtualscreen.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_energy no descriptor_use_footprint_similarity yes descriptor_use_pharmacophore_score yes descriptor_use_tanimoto yes descriptor_use_hungarian yes descriptor_use_volume_overlap yes descriptor_fps_use_footprint_reference_mol2 yes descriptor_fps_footprint_reference_mol2_filename ../004.dock/3jqz.lig.min_scored.mol2 descriptor_fps_foot_compare_type Euclidean descriptor_fps_normalize_foot no descriptor_fps_foot_comp_all_residue yes descriptor_fps_receptor_filename ../001.dockprep/3jqz.rec.withH.charged.mol2 descriptor_fps_vdw_att_exp 6 descriptor_fps_vdw_rep_exp 12 descriptor_fps_vdw_rep_rad_scale 1 descriptor_fps_use_distance_dependent_dielectric yes descriptor_fps_dielectric 4.0 descriptor_fps_vdw_fp_scale 1 descriptor_fps_es_fp_scale 1 descriptor_fps_hb_fp_scale 0 descriptor_fms_score_use_ref_mol2 yes descriptor_fms_score_ref_mol2_filename ../004.dock/3jqz.lig.min_scored.mol2 descriptor_fms_score_write_reference_pharmacophore_mol2 no descriptor_fms_score_write_reference_pharmacophore_txt no descriptor_fms_score_write_candidate_pharmacophore no descriptor_fms_score_write_matched_pharmacophore no descriptor_fms_score_compare_type overlap descriptor_fms_score_full_match yes descriptor_fms_score_match_rate_weight 5.0 descriptor_fms_score_match_dist_cutoff 1.0 descriptor_fms_score_match_proj_cutoff 0.7071 descriptor_fms_score_max_score 20 descriptor_fingerprint_ref_filename ../004.dock/3jqz.lig.min_scored.mol2 descriptor_hungarian_ref_filename ../004.dock/3jqz.lig.min_scored.mol2 descriptor_hungarian_matching_coeff -5 descriptor_hungarian_rmsd_coeff 1 descriptor_volume_reference_mol2_filename ../004.dock/3jqz.lig.min_scored.mol2 descriptor_volume_overlap_compute_method analytical descriptor_weight_fps_score 1 descriptor_weight_pharmacophore_score 1 descriptor_weight_fingerprint_tanimoto -1 descriptor_weight_hungarian 1 descriptor_weight_volume_overlap_score -1 gbsa_zou_score_secondary no gbsa_hawkins_score_secondary no SASA_score_secondary no amber_score_secondary no minimize_ligand no atom_model all vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl chem_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/chem.defn pharmacophore_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/ph4.defn ligand_outfile_prefix descriptor.out write_footprints yes write_hbonds yes write_orientations no num_scored_conformers 1 rank_ligands no
