Difference between revisions of "2019 DOCK tutorial 2 with PDBID 2P16"
From Rizzo_Lab
Stonybrook (talk | contribs) (→Prepare the Receptor File) |
(→VI. Virtual Screen) |
||
| (88 intermediate revisions by 2 users not shown) | |||
| Line 14: | Line 14: | ||
# The resulting configurations are 'pruned' and energy re-minimized, yielding the docked configurations. | # The resulting configurations are 'pruned' and energy re-minimized, yielding the docked configurations. | ||
| + | Users of this tutorial may want to refer to a complete DOCK [http://dock.compbio.ucsf.edu/DOCK_6/dock6_manual.htm User Manual] for more detailed instructions. | ||
===2P16=== | ===2P16=== | ||
| Line 33: | Line 34: | ||
Second, use mkdir command to create the following folders: | Second, use mkdir command to create the following folders: | ||
| − | mkdir 001.files 002.surface_spheres 003.gridbox 004.dock 005. | + | mkdir 001.files 002.surface_spheres 003.gridbox 004.dock 005.virtual_screen 006.virtual_screen_mpi 007.cartesianmin 008.rescore |
Then, you can use the following command to check the folders you just created | Then, you can use the following command to check the folders you just created | ||
| Line 53: | Line 54: | ||
or simply using the graphical user interface of chimera. | or simply using the graphical user interface of chimera. | ||
| − | [[File:wiki_2p16_noh_com.png|thumb|center| | + | [[File:wiki_2p16_noh_com.png|thumb|center|500px|2p16.pdb]] |
Identify the main components of the model (receptor, ligand, solvent, surfactants, metal ions). In our case, the protein contains two separate chains. You can also find a Mg ion shown in green sphere, and some water molecules. Because there is no hydrogen atom in the pdb file, the water molecules are only shown in oxygen atoms. By default, only 3 waters are shown. You can use Actions -> Atoms/Bonds -> Show to view all the waters. | Identify the main components of the model (receptor, ligand, solvent, surfactants, metal ions). In our case, the protein contains two separate chains. You can also find a Mg ion shown in green sphere, and some water molecules. Because there is no hydrogen atom in the pdb file, the water molecules are only shown in oxygen atoms. By default, only 3 waters are shown. You can use Actions -> Atoms/Bonds -> Show to view all the waters. | ||
| Line 68: | Line 69: | ||
Then do Select -> Structure -> protein to selected the protein. Hold Ctrl and Shift button, using mouse to left-click the ligand. This would include a atom or a bond of the ligand in to the selection. Then press up arrow to select to whole ligand. With Ctrl and Shift holding, left-click the water (oxygen atom) which is in the binding pocket. Now the selected structure should include the whole protein, the ligand and one oxygen atom. Do Select -> Invert (all models) to invert the selection. Now all the waters except one will be selected. Delete those molecules. Do File -> Save MOL2 to save the structure as 2p16_noh_com.mol2. | Then do Select -> Structure -> protein to selected the protein. Hold Ctrl and Shift button, using mouse to left-click the ligand. This would include a atom or a bond of the ligand in to the selection. Then press up arrow to select to whole ligand. With Ctrl and Shift holding, left-click the water (oxygen atom) which is in the binding pocket. Now the selected structure should include the whole protein, the ligand and one oxygen atom. Do Select -> Invert (all models) to invert the selection. Now all the waters except one will be selected. Delete those molecules. Do File -> Save MOL2 to save the structure as 2p16_noh_com.mol2. | ||
| − | [[File:wiki_2p16_noh_complex.png|thumb|center| | + | [[File:wiki_2p16_noh_complex.png|thumb|center|500px|2p16_noh_com.mol2]] |
====Prepare the Ligand File==== | ====Prepare the Ligand File==== | ||
Do Select -> Structure -> ligand to select the ligand. Do Select -> Invert (all models), then delete the selected atoms. This will left the ligand molecule only. Save it as 2p16_noh_lig.mol2. | Do Select -> Structure -> ligand to select the ligand. Do Select -> Invert (all models), then delete the selected atoms. This will left the ligand molecule only. Save it as 2p16_noh_lig.mol2. | ||
| − | [[File:wiki_2p16_noh_lig.png|thumb|center| | + | [[File:wiki_2p16_noh_lig.png|thumb|center|500px|2p16_noh_lig.mol2]] |
====Prepare the Receptor File==== | ====Prepare the Receptor File==== | ||
Do Select -> Structure -> ligand to select the ligand. Then delete the selected atoms. This will left the protein and one water. Save it as 2p16_noh_rec.mol2. | Do Select -> Structure -> ligand to select the ligand. Then delete the selected atoms. This will left the protein and one water. Save it as 2p16_noh_rec.mol2. | ||
| − | [[File: | + | [[File:wiki_2p16_noh_rec 2.png|thumb|center|500px|2p16_noh_rec.mol2]] |
====Adding hydrogen and charge==== | ====Adding hydrogen and charge==== | ||
| Line 84: | Line 85: | ||
Tools -> surface/Binding Analysis -> Dock Prep | Tools -> surface/Binding Analysis -> Dock Prep | ||
De-check the first row "Delete Solvent" (it's already done by yourself manually). Then click 'OK'. Keep pressing 'OK' in the next window. In the following window, make sure we are using AMBER ff14SB and AM1-BCC, then click 'OK'. It will add hydrogen atoms and partial charges to the receptor. Save it as 2p16_wh_rec.mol2. | De-check the first row "Delete Solvent" (it's already done by yourself manually). Then click 'OK'. Keep pressing 'OK' in the next window. In the following window, make sure we are using AMBER ff14SB and AM1-BCC, then click 'OK'. It will add hydrogen atoms and partial charges to the receptor. Save it as 2p16_wh_rec.mol2. | ||
| + | |||
| + | [[File:wiki_2p16_wh_rec.png|thumb|center|500px|2p16_wh_rec.mol2]] | ||
| + | |||
Open 2p16_noh_lig.mol2, follow the same steps for rec file, add hydrogen and charges, and save the file as 2p16_wh_lig.mol2. | Open 2p16_noh_lig.mol2, follow the same steps for rec file, add hydrogen and charges, and save the file as 2p16_wh_lig.mol2. | ||
| + | |||
| + | [[File:wiki_2p16_wh_lig.png|thumb|center|500px|2p16_wh_lig.mol2]] | ||
=III. Generating receptor surface and spheres= | =III. Generating receptor surface and spheres= | ||
| + | |||
| + | ===Prepare the environment=== | ||
| + | |||
| + | vi ~/.bashrc | ||
| + | |||
| + | add the following lines to the bashrc file | ||
| + | |||
| + | export DOCKHOME="/gpfs/projects/AMS536/zzz.programs/dock6/" | ||
| + | export PATH=$DOCKHOME/bin\:$PATH | ||
| + | |||
| + | save the file and quit. Then do | ||
| + | |||
| + | source ~/.bashrc | ||
| + | |||
===Preparation of DMS file=== | ===Preparation of DMS file=== | ||
Open 2p16_noh_rec.mol2 using chimera. Do | Open 2p16_noh_rec.mol2 using chimera. Do | ||
| Line 98: | Line 118: | ||
Transfer all the files created so far to SeaWulf cluster to be used in DOCK. You can use the following command: | Transfer all the files created so far to SeaWulf cluster to be used in DOCK. You can use the following command: | ||
| − | scp 2p16* username@login.seawulf.stonybrook.edu:/ | + | scp 2p16* username@login.seawulf.stonybrook.edu:/gpfs/projects/AMS536/your_directory/ |
From now on, we will be working on the SeaWulf clusters instead of your own machine. | From now on, we will be working on the SeaWulf clusters instead of your own machine. | ||
===Generating spheres=== | ===Generating spheres=== | ||
| − | + | Go to 002.surface_spheres folder | |
| − | + | Create a new input file to create spheres by | |
| + | vi INSPH | ||
| + | |||
| + | then type the following lines inside the file. | ||
| − | + | ../001.files/2p16_noh_rec.dms | |
| − | + | R | |
| − | + | X | |
| − | + | 0.0 | |
| − | + | 4.0 | |
| − | + | 1.4 | |
| − | + | 2p16_rec.sph | |
The first line 2nnq_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 2nnq_spheres.sph creates the sph file that contains clustered spheres. | The first line 2nnq_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 2nnq_spheres.sph creates the sph file that contains clustered spheres. | ||
| Line 120: | Line 143: | ||
sphgen -i INSPH -o OUTSPH | sphgen -i INSPH -o OUTSPH | ||
| − | Once sphgen command is successful, | + | Once sphgen command is successful, 2p16_rec.sph file will be created. Open it up using Chimera along with 2p16_noh_rec.mol2 file. You should get a similar output like the image below. |
| − | [[File: | + | [[File:wiki_2p16_noh_rec_sphere.png|thumb|center|500px|2p16 sphere]] |
===Selecting Spheres=== | ===Selecting Spheres=== | ||
| Line 128: | Line 151: | ||
Here we will be selecting the spheres which defines 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 type the following command. | Here we will be selecting the spheres which defines 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 type the following command. | ||
| − | sphere_selector | + | sphere_selector 2p16_rec.sph ../001.files/2p16_noh_lig.mol2 10.0 |
This command will select all of the spheres within 10.0 angstroms of the ligand and output them to selected_spheres.sph. 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: | + | [[File:wiki_2p16_noh_rec_selectedsphere.png|thumb|center|500px|2p16 selected sphere]] |
=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. | ||
Y | Y | ||
8.0 | 8.0 | ||
| − | ../ | + | ../002.surface_spheres/selected_spheres.sph |
1 | 1 | ||
| − | + | 2p16.box.pdb | |
| + | |||
| + | Each of the above lines means: | ||
| + | |||
| + | Say 'Yes' to make a box | ||
| + | The boundary of the box is at least 8 Angstroms from any spheres | ||
| + | The location and file name of the sphere file | ||
| + | Use the first cluster of spheres (Only 1 cluster in this case) | ||
| + | The output filename | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
Use the following command to generate the box. | Use the following command to generate the box. | ||
| Line 156: | Line 182: | ||
showbox < showbox.in | showbox < showbox.in | ||
| − | If this step is successful, you should see a new file ( | + | If this step is successful, you should see a new file (2p16.box.pdb) in 003.gridbox folder. |
| + | |||
| + | [[File:wiki_2p16_noh_rec_selectedsphere_box.png|thumb|center|500px|2p16 selected sphere with box]] | ||
===Generating grid=== | ===Generating grid=== | ||
| − | + | In the same folder, create a new file | |
| − | + | vi grid.in | |
| − | |||
| − | + | Copy the following lines in to the file. | |
compute_grids yes | compute_grids yes | ||
| Line 174: | Line 201: | ||
atom_model a | atom_model a | ||
attractive_exponent 6 | attractive_exponent 6 | ||
| − | repulsive_exponent | + | repulsive_exponent 9 |
distance_dielectric yes | distance_dielectric yes | ||
dielectric_factor 4 | dielectric_factor 4 | ||
bump_filter yes | bump_filter yes | ||
bump_overlap 0.75 | bump_overlap 0.75 | ||
| − | receptor_file ../ | + | receptor_file ../001.files/2p16_wh_rec.mol2 |
| − | box_file | + | box_file 2p16.box.pdb |
vdw_definition_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn | vdw_definition_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn | ||
score_grid_prefix grid | score_grid_prefix grid | ||
| − | If the command is successful, three new files will be generated. ( | + | Use the following command to generate the grid. |
| + | |||
| + | grid -i grid.in -o grid.out | ||
| + | |||
| + | If the command is successful, three new files will be generated. (grid.out, grid.nrg, grid.bmp). Go through grid.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. | ||
| + | |||
| + | All the generated grid files are binary files. It's hard for us to visualize them. | ||
| + | |||
| + | Please note that we use 9 instead of 12 as the repulsive exponent in Lennard-Jones potential. This gives a 'softer' repulsion. | ||
=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 2p16.pdb will be re-docked into the receptor. 3 Methods will be used to achieve this. |
1. rigid docking | 1. rigid docking | ||
| Line 195: | Line 230: | ||
3. flexible 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. | 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 | + | Go to the directory 004.dock and a create a new file |
| + | |||
| + | touch min.in | ||
| + | |||
| + | and enter the command below. | ||
| + | |||
dock6 -i min.in | dock6 -i min.in | ||
| Line 207: | Line 247: | ||
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.files/2p16_wh_lig.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 yes | calculate_rmsd yes | ||
| − | use_rmsd_reference_mol ../ | + | use_rmsd_reference_mol ../001.files/2p16_wh_lig.mol2 |
use_database_filter no | use_database_filter no | ||
orient_ligand no | orient_ligand no | ||
| Line 224: | Line 264: | ||
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 251: | Line 291: | ||
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 2p16.lig.min |
write_orientations no | write_orientations no | ||
num_scored_conformers 1 | num_scored_conformers 1 | ||
rank_ligands no | rank_ligands no | ||
| − | If the | + | If you made any errors during the entering of the above lines, just keep doing. After you have finished every questions, you can use vi to correct your min.in file. Then you are ready for the minimization. Do the following: |
| − | [[File: | + | dock6 -i min.in -o min.out |
| + | |||
| + | If the process is successful a new file (2p16.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:wiki_2p16_wh_lig_min_mol2.png|thumb|center|500px|2p16_receptor with the crystalized ligand(green) and the minimized ligand(pink).]] | ||
| + | |||
| + | We can also use footprint to visualize the changes after minimization. Under the same directory, create a new file | ||
| + | |||
| + | touch footprint.in | ||
| + | |||
| + | then | ||
| + | |||
| + | dock6 -i footprint.in | ||
| + | |||
| + | Answer the following questions interactively | ||
| + | |||
| + | conformer_search_type rigid | ||
| + | use_internal_energy no | ||
| + | ligand_atom_file ./2p16.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 ../001.files/2p16_wh_lig.mol2 | ||
| + | fps_score_foot_compare_type Euclidean | ||
| + | fps_score_normalize_foot no | ||
| + | fps_score_foot_comp_all_residue yes | ||
| + | fps_score_receptor_filename ../001.files/2p16_wh_rec.mol2 | ||
| + | fps_score_vdw_att_exp 6 | ||
| + | fps_score_vdw_rep_exp 9 | ||
| + | fps_score_vdw_rep_rad_scale 1 | ||
| + | fps_score_use_distance_dependent_dielectric yes | ||
| + | fps_score_dielectric 4.0 | ||
| + | fps_score_vdw_fp_scale 1 | ||
| + | fps_score_es_fp_scale 1 | ||
| + | fps_score_hb_fp_scale 0 | ||
| + | pharmacophore_score_secondary no | ||
| + | descriptor_score_secondary no | ||
| + | gbsa_zou_score_secondary no | ||
| + | gbsa_hawkins_score_secondary no | ||
| + | SASA_score_secondary no | ||
| + | amber_score_secondary no | ||
| + | minimize_ligand no | ||
| + | atom_model all | ||
| + | vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn | ||
| + | flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn | ||
| + | flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl | ||
| + | ligand_outfile_prefix fps.min.output | ||
| + | write_footprints yes | ||
| + | write_hbonds yes | ||
| + | write_orientations no | ||
| + | num_scored_conformers 1 | ||
| + | rank_ligands no | ||
| + | |||
| + | If dock runs successfully, you can find the following files generated | ||
| + | |||
| + | fps.min.output_footprint_scored.txt | ||
| + | fps.min.output_hbond_scored.txt | ||
| + | |||
| + | In order to visualize the footprint, first we need to copy a python script into our directory, then run it: | ||
| + | |||
| + | cp /gpfs/projects/AMS536/zzz.programs/plot_footprint_single_magnitude.py ./ | ||
| + | python plot_footprint_single_magnitude.py fps.min.output_footprint_scored.txt 50 | ||
| + | |||
| + | A pdf document named fps.min.output_footprint_scored.txt.pdf will be generated. It looks like this: | ||
| + | |||
| + | [[File:Wiki 2p16 Footprint.png |thumb|center|1000px|footprint comparison with energy minimized (red) and the un-minimized ligand(blue).]] | ||
==Rigid Docking== | ==Rigid Docking== | ||
| Line 271: | Line 394: | ||
conformer_search_type rigid | conformer_search_type rigid | ||
use_internal_energy yes | use_internal_energy yes | ||
| − | ligand_atom_file | + | internal_energy_rep_exp 9 |
| + | internal_energy_cutoff 100.0 | ||
| + | ligand_atom_file ../001.files/2p16_wh_lig.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 yes | calculate_rmsd yes | ||
| − | use_rmsd_reference_mol | + | use_rmsd_reference_mol yes |
| + | rmsd_reference_filename ../001.files/2p16_wh_lig.mol2 | ||
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 294: | Line 420: | ||
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 307: | Line 433: | ||
minimize_ligand yes | minimize_ligand yes | ||
simplex_max_iterations 1000 | simplex_max_iterations 1000 | ||
| − | simplex_tors_premin_iterations | + | simplex_tors_premin_iterations 2 |
simplex_max_cycles 1 | simplex_max_cycles 1 | ||
simplex_score_converge 0.1 | simplex_score_converge 0.1 | ||
| Line 322: | Line 448: | ||
ligand_outfile_prefix rigid.out | ligand_outfile_prefix rigid.out | ||
write_orientations no | write_orientations no | ||
| − | num_scored_conformers | + | num_scored_conformers 40 |
| + | write_conformations yes | ||
| + | cluster_conformations yes | ||
| + | cluster_rmsd_threshold 2.0 | ||
rank_ligands no | rank_ligands no | ||
| Line 328: | Line 457: | ||
Open Chimera | Open Chimera | ||
| − | File -> Open -> | + | File -> Open -> 2p16_rec_withH.mol2 |
| − | File -> Open -> | + | File -> Open -> 2p16_lig_withH.mol2 |
Tools -> Surface/binding Analysis -> ViewDock -> Select the Rigid Dock output file. (rigid.out_scored.mol2) | Tools -> Surface/binding Analysis -> ViewDock -> Select the Rigid Dock output file. (rigid.out_scored.mol2) | ||
In the loaded dialog box select Dock4,5 or 6 | In the loaded dialog box select Dock4,5 or 6 | ||
| Line 339: | Line 468: | ||
Follow the same steps to get all the properties | Follow the same steps to get all the properties | ||
| − | Your visualized structure should be similar to the image below. | + | You can choose different conformations from the ViewDock window. The top one should be the conformation with lowest energy. Your visualized structure should be similar to the image below. |
| − | [[File: | + | [[File:2p16_rigid_3.png|thumb|center|500px|Rigid docking results for 2p16, Reference and docked ligand is colored with sky blue and pink, respectively]] |
| + | [[File:2p16_rigid_2_2.png|thumb|center|500px|Rigid docking results for 2p16]] | ||
==Fixed Anchor Docking== | ==Fixed Anchor Docking== | ||
Create an input file for fixed anchor docking. | Create an input file for fixed anchor docking. | ||
| − | touch | + | |
| + | touch fad.in | ||
| + | |||
Use the input file to perform fixed anchor docking | Use the input file to perform fixed anchor docking | ||
| − | dock6 -i | + | |
| + | dock6 -i fad.in | ||
Use the following lines to answer the prompted questions as we did in rigid docking. | Use the following lines to answer the prompted questions as we did in rigid docking. | ||
conformer_search_type flex | conformer_search_type flex | ||
| + | write_fragment_libraries no | ||
user_specified_anchor no | user_specified_anchor no | ||
limit_max_anchors no | limit_max_anchors no | ||
| Line 362: | Line 496: | ||
use_clash_overlap no | use_clash_overlap no | ||
write_growth_tree no | write_growth_tree no | ||
| − | |||
use_internal_energy yes | use_internal_energy yes | ||
| − | internal_energy_rep_exp | + | internal_energy_rep_exp 9 |
internal_energy_cutoff 100.0 | internal_energy_cutoff 100.0 | ||
| − | ligand_atom_file ../ | + | ligand_atom_file ../001.files/2p16_wh_lig.mol2 |
limit_max_ligands no | limit_max_ligands no | ||
skip_molecule no | skip_molecule no | ||
| Line 372: | Line 505: | ||
calculate_rmsd yes | calculate_rmsd yes | ||
use_rmsd_reference_mol yes | use_rmsd_reference_mol yes | ||
| − | rmsd_reference_filename ../ | + | rmsd_reference_filename ../001.files/2p16_wh_lig.mol2 |
use_database_filter no | use_database_filter no | ||
| − | orient_ligand 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 | bump_filter no | ||
score_molecules yes | score_molecules yes | ||
| Line 384: | Line 523: | ||
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 409: | Line 548: | ||
simplex_grow_tors_premin_iterations 0 | simplex_grow_tors_premin_iterations 0 | ||
simplex_random_seed 0 | simplex_random_seed 0 | ||
| − | simplex_restraint_min | + | simplex_restraint_min yes |
| + | simplex_coefficient_restraint 10.0 | ||
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 | ||
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 2p16_fad |
write_orientations no | write_orientations no | ||
| − | num_scored_conformers | + | num_scored_conformers 40 |
write_conformations no | write_conformations no | ||
cluster_conformations yes | cluster_conformations yes | ||
| Line 422: | Line 562: | ||
rank_ligands no | rank_ligands no | ||
| − | Once docking is completed an output file will be generated. ( | + | Once docking is completed an output file will be generated. (2p16_fad_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. |
| − | [[File: | + | [[File:2p16_fixed.png|thumb|center|500px| Poses generated for fixed anchor docking, Reference and docked ligand is colored with sky blue and pink, respectively]] |
| + | [[File:2p16_fixed_2.png|thumb|center|500px| Poses generated for fixed anchor docking]] | ||
==Flexible Docking== | ==Flexible Docking== | ||
| Line 449: | Line 590: | ||
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 2p16.lig.min_scored.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 yes | calculate_rmsd yes | ||
| − | use_rmsd_reference_mol | + | use_rmsd_reference_mol yes |
| + | rmsd_reference_filename 2p16.lig.min_scored.mol2 | ||
use_database_filter no | use_database_filter no | ||
orient_ligand yes | orient_ligand yes | ||
| Line 510: | Line 652: | ||
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. | 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. | ||
| − | [[File: | + | [[File:2p16_flex.png|thumb|center|500px| Flexible docking results for 2p16, Reference and docked ligand is colored with sky blue and pink, respectively]] |
| + | [[File:2p16_flex_2.png|thumb|center|500px| Flexible docking results for 2p16]] | ||
| − | = | + | = 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 005.virtual_screen. | |
| − | + | Copy the ligand library to the same directory. | |
| − | + | cp /gpfs/projects/AMS536/zzz.programs/VS_library_5K.mol2 ./ | |
| − | + | Create a new input file for virtual screen. | |
| − | |||
| − | + | 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. | ||
| + | |||
conformer_search_type flex | conformer_search_type flex | ||
| + | write_fragment_libraries no | ||
user_specified_anchor no | user_specified_anchor no | ||
limit_max_anchors no | limit_max_anchors no | ||
| Line 606: | Line 686: | ||
use_clash_overlap no | use_clash_overlap no | ||
write_growth_tree no | write_growth_tree no | ||
| − | |||
use_internal_energy yes | use_internal_energy yes | ||
| − | internal_energy_rep_exp | + | internal_energy_rep_exp 9 |
internal_energy_cutoff 100.0 | internal_energy_cutoff 100.0 | ||
| − | ligand_atom_file | + | ligand_atom_file VS_library_5K.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 | ||
chemical_matching no | chemical_matching no | ||
| − | use_ligand_spheres no | + | use_ligand_spheres no |
| − | bump_filter no | + | bump_filter no |
score_molecules yes | score_molecules yes | ||
contact_score_primary no | contact_score_primary no | ||
| Line 634: | Line 711: | ||
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 643: | Line 720: | ||
gbsa_zou_score_secondary no | gbsa_zou_score_secondary no | ||
gbsa_hawkins_score_secondary no | gbsa_hawkins_score_secondary no | ||
| − | SASA_score_secondary no | + | SASA_score_secondary no |
amber_score_secondary no | amber_score_secondary no | ||
minimize_ligand yes | minimize_ligand yes | ||
| Line 654: | Line 731: | ||
simplex_trans_step 1.0 | simplex_trans_step 1.0 | ||
simplex_rot_step 0.1 | simplex_rot_step 0.1 | ||
| − | simplex_tors_step 10.0 | + | simplex_tors_step 10.0 |
simplex_anchor_max_iterations 500 | simplex_anchor_max_iterations 500 | ||
simplex_grow_max_iterations 500 | simplex_grow_max_iterations 500 | ||
| Line 665: | Line 742: | ||
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 virtual.out | ligand_outfile_prefix virtual.out | ||
| − | write_orientations no | + | write_orientations no |
num_scored_conformers 1 | num_scored_conformers 1 | ||
rank_ligands no | rank_ligands no | ||
| − | + | ||
| − | + | You can run this job on head node for a short while, just to make sure everything is fine. Then you can kill the job by Ctrl-C. We will use the mpi version to finish this calculation. See next section. | |
=VII.Virtual Screen (MPI)= | =VII.Virtual Screen (MPI)= | ||
| Line 677: | Line 754: | ||
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) | ||
| − | + | ||
| + | vi virtual.sh | ||
| + | |||
Type the following lines in the new file. | Type the following lines in the new file. | ||
| + | |||
#!/bin/bash | #!/bin/bash | ||
#PBS -l walltime=48:00:00 | #PBS -l walltime=48:00:00 | ||
#PBS -l nodes=4:ppn=28 | #PBS -l nodes=4:ppn=28 | ||
#PBS -q long | #PBS -q long | ||
| − | #PBS -N | + | #PBS -N 2p16.virtual |
#PBS -V | #PBS -V | ||
cd $PBS_O_WORKDIR | cd $PBS_O_WORKDIR | ||
| − | mpirun -np 112 dock6.mpi -i virtual.in -o | + | mpirun -v -np 112 dock6.mpi -i virtual.in -o 2p16.virtual.mpi.out |
| + | |||
| + | Before submitting the job, please check the mpirun you are using | ||
| + | |||
| + | which mpirun | ||
| + | |||
| + | make sure you are using an intel mpi. | ||
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. | ||
Latest revision as of 11:35, 12 December 2019
This tutorial teaches you how to dock a drug molecule to a receptor.
/*This page is under construction*/
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.
Users of this tutorial may want to refer to a complete DOCK User Manual for more detailed instructions.
2P16
2P16 is the crystal structure of factor Xa in complex with the inhibitor APIXABAN (tradename Eliquis). You can get the pdb file from here 2P16. Click Download Files --> PDB format to download the PDB file to your directory. You may want to use a file editor like vi to check the content of a pdb file.
vi 2p16.pdb
There are usually many useful informations other than structure data in a pdb file.
To quit vi, use
:q
Organization of Directories
In this part, we are going to unify our nomenclature of directories. It's always helpful to maintain a good order of files.
First, cd into your working directory on SeaWulf.
Second, use mkdir command to create the following folders:
mkdir 001.files 002.surface_spheres 003.gridbox 004.dock 005.virtual_screen 006.virtual_screen_mpi 007.cartesianmin 008.rescore
Then, you can use the following command to check the folders you just created
ls -l
Files created in different steps will be saved to different folders accordingly in this tutorial.
II. Preparation of the ligand and receptor
Download the pdb file 2P16 from PDB database save it in 001.files folder.
Checking the structure
Read the article related to the PDB file [1] to understand protonation states, charges, environmental conditions and other important information regarding the receptor and the ligand.
Open the pdb file through chimera and look at the structure. You can do it by using the following command in shell
chimera 2p16.pdb
or simply using the graphical user interface of chimera.
Identify the main components of the model (receptor, ligand, solvent, surfactants, metal ions). In our case, the protein contains two separate chains. You can also find a Mg ion shown in green sphere, and some water molecules. Because there is no hydrogen atom in the pdb file, the water molecules are only shown in oxygen atoms. By default, only 3 waters are shown. You can use Actions -> Atoms/Bonds -> Show to view all the waters.
Carefully look to identify if there are any missing residues or missing loops. (This particular PDB file didn't contain any missing loops or missing residues)
Preparation of receptor
In this part we are going to generate the receptor file and the ligand file that will going to be used in docking by chimera.
Clean the structure
As you have seen, the pdb file is not very clean. We need to do some modifications to make it suitable for docking studies. With the pdb file loaded by chimera, do Select -> Chain -> L to select the smaller chain. Use Actions -> Atoms/Bonds -> Delete to delete it. Then do Select -> Structure -> Ions to select the Mg ion. Use the same method to delete it.
Then do Select -> Structure -> protein to selected the protein. Hold Ctrl and Shift button, using mouse to left-click the ligand. This would include a atom or a bond of the ligand in to the selection. Then press up arrow to select to whole ligand. With Ctrl and Shift holding, left-click the water (oxygen atom) which is in the binding pocket. Now the selected structure should include the whole protein, the ligand and one oxygen atom. Do Select -> Invert (all models) to invert the selection. Now all the waters except one will be selected. Delete those molecules. Do File -> Save MOL2 to save the structure as 2p16_noh_com.mol2.
Prepare the Ligand File
Do Select -> Structure -> ligand to select the ligand. Do Select -> Invert (all models), then delete the selected atoms. This will left the ligand molecule only. Save it as 2p16_noh_lig.mol2.
Prepare the Receptor File
Do Select -> Structure -> ligand to select the ligand. Then delete the selected atoms. This will left the protein and one water. Save it as 2p16_noh_rec.mol2.
Adding hydrogen and charge
Now we are going to add hydrogen atoms and charges to our receptor and ligand. Open 2p16_noh_rec.mol2 file again using Chimera and use the following instructions to prepare the receptor file to be used in DOCK.
Tools -> surface/Binding Analysis -> Dock Prep
De-check the first row "Delete Solvent" (it's already done by yourself manually). Then click 'OK'. Keep pressing 'OK' in the next window. In the following window, make sure we are using AMBER ff14SB and AM1-BCC, then click 'OK'. It will add hydrogen atoms and partial charges to the receptor. Save it as 2p16_wh_rec.mol2.
Open 2p16_noh_lig.mol2, follow the same steps for rec file, add hydrogen and charges, and save the file as 2p16_wh_lig.mol2.
III. Generating receptor surface and spheres
Prepare the environment
vi ~/.bashrc
add the following lines to the bashrc file
export DOCKHOME="/gpfs/projects/AMS536/zzz.programs/dock6/" export PATH=$DOCKHOME/bin\:$PATH
save the file and quit. Then do
source ~/.bashrc
Preparation of DMS file
Open 2p16_noh_rec.mol2 using chimera. Do
Action -> Surface -> Show Tools -> Structure Editing -> Write DMS
Save the structure as 2p16_noh_rec.dms.
Reopen the file using chimera and make sure the surface was generated.
Transfer all the files created so far to SeaWulf cluster to be used in DOCK. You can use the following command:
scp 2p16* username@login.seawulf.stonybrook.edu:/gpfs/projects/AMS536/your_directory/
From now on, we will be working on the SeaWulf clusters instead of your own machine.
Generating spheres
Go to 002.surface_spheres folder Create a new input file to create spheres by
vi INSPH
then type the following lines inside the file.
../001.files/2p16_noh_rec.dms R X 0.0 4.0 1.4 2p16_rec.sph
The first line 2nnq_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 2nnq_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, 2p16_rec.sph file will be created. Open it up using Chimera along with 2p16_noh_rec.mol2 file. You should get a similar output like the image below.
Selecting Spheres
Here we will be selecting the spheres which defines 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 type the following command.
sphere_selector 2p16_rec.sph ../001.files/2p16_noh_lig.mol2 10.0
This command will select all of the spheres within 10.0 angstroms of the ligand and output them to selected_spheres.sph. 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 2p16.box.pdb
Each of the above lines means:
Say 'Yes' to make a box The boundary of the box is at least 8 Angstroms from any spheres The location and file name of the sphere file Use the first cluster of spheres (Only 1 cluster in this case) The output filename
Use the following command to generate the box.
showbox < showbox.in
If this step is successful, you should see a new file (2p16.box.pdb) in 003.gridbox folder.
Generating grid
In the same folder, create a new file
vi grid.in
Copy the following lines in to the 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.files/2p16_wh_rec.mol2 box_file 2p16.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 grid.out
If the command is successful, three new files will be generated. (grid.out, grid.nrg, grid.bmp). Go through grid.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.
All the generated grid files are binary files. It's hard for us to visualize them.
Please note that we use 9 instead of 12 as the repulsive exponent in Lennard-Jones potential. This gives a 'softer' repulsion.
V. Docking a single molecule for pose reproduction
Under this section, the ligand for 2p16.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 004.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.files/2p16_wh_lig.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd yes use_rmsd_reference_mol ../001.files/2p16_wh_lig.mol2 use_database_filter no orient_ligand no bump_filter no score_molecules yes contact_score_primary no contact_score_secondary no grid_score_primary yes grid_score_secondary no grid_score_rep_rad_scale 1 grid_score_vdw_scale 1 grid_score_es_scale 1 grid_score_grid_prefix ../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 2p16.lig.min write_orientations no num_scored_conformers 1 rank_ligands no
If you made any errors during the entering of the above lines, just keep doing. After you have finished every questions, you can use vi to correct your min.in file. Then you are ready for the minimization. Do the following:
dock6 -i min.in -o min.out
If the process is successful a new file (2p16.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.
We can also use footprint to visualize the changes after minimization. Under the same directory, create a new file
touch footprint.in
then
dock6 -i footprint.in
Answer the following questions interactively
conformer_search_type rigid use_internal_energy no ligand_atom_file ./2p16.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 ../001.files/2p16_wh_lig.mol2 fps_score_foot_compare_type Euclidean fps_score_normalize_foot no fps_score_foot_comp_all_residue yes fps_score_receptor_filename ../001.files/2p16_wh_rec.mol2 fps_score_vdw_att_exp 6 fps_score_vdw_rep_exp 9 fps_score_vdw_rep_rad_scale 1 fps_score_use_distance_dependent_dielectric yes fps_score_dielectric 4.0 fps_score_vdw_fp_scale 1 fps_score_es_fp_scale 1 fps_score_hb_fp_scale 0 pharmacophore_score_secondary no descriptor_score_secondary no gbsa_zou_score_secondary no gbsa_hawkins_score_secondary no SASA_score_secondary no amber_score_secondary no minimize_ligand no atom_model all vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix fps.min.output write_footprints yes write_hbonds yes write_orientations no num_scored_conformers 1 rank_ligands no
If dock runs successfully, you can find the following files generated
fps.min.output_footprint_scored.txt fps.min.output_hbond_scored.txt
In order to visualize the footprint, first we need to copy a python script into our directory, then run it:
cp /gpfs/projects/AMS536/zzz.programs/plot_footprint_single_magnitude.py ./ python plot_footprint_single_magnitude.py fps.min.output_footprint_scored.txt 50
A pdf document named fps.min.output_footprint_scored.txt.pdf will be generated. It looks like this:
Rigid Docking
Create an input file for rigid docking
touch rigid.in
Run dock using the created input file.
dock6 -i rigid.in
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 internal_energy_rep_exp 9 internal_energy_cutoff 100.0 ligand_atom_file ../001.files/2p16_wh_lig.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd yes use_rmsd_reference_mol yes rmsd_reference_filename ../001.files/2p16_wh_lig.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 2 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 40 write_conformations yes cluster_conformations yes cluster_rmsd_threshold 2.0 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 -> 2p16_rec_withH.mol2 File -> Open -> 2p16_lig_withH.mol2 Tools -> Surface/binding Analysis -> ViewDock -> Select the Rigid Dock output file. (rigid.out_scored.mol2) In the loaded dialog box 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
You can choose different conformations from the ViewDock window. The top one should be the conformation with lowest energy. Your visualized structure should be similar to the image below.
Fixed Anchor Docking
Create an input file for fixed anchor docking.
touch fad.in
Use the input file to perform fixed anchor docking
dock6 -i fad.in
Use the following lines to answer the prompted questions as we did in rigid docking.
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 9 internal_energy_cutoff 100.0 ligand_atom_file ../001.files/2p16_wh_lig.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd yes use_rmsd_reference_mol yes rmsd_reference_filename ../001.files/2p16_wh_lig.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 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 2p16_fad write_orientations no num_scored_conformers 40 write_conformations no cluster_conformations yes cluster_rmsd_threshold 2.0 rank_ligands no
Once docking is completed an output file will be generated. (2p16_fad_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
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 2p16.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 2p16.lig.min_scored.mol2 use_database_filter no orient_ligand yes automated_matching yes receptor_site_file ../2.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 ../3.boxgrid/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.
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 005.virtual_screen.
Copy the ligand library to the same directory.
cp /gpfs/projects/AMS536/zzz.programs/VS_library_5K.mol2 ./
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
write_fragment_libraries no
user_specified_anchor no
limit_max_anchors no
min_anchor_size 5
pruning_use_clustering yes
pruning_max_orients 1000
pruning_clustering_cutoff 100
pruning_conformer_score_cutoff 100.0
pruning_conformer_score_scaling_factor 1.0
use_clash_overlap no
write_growth_tree no
use_internal_energy yes
internal_energy_rep_exp 9
internal_energy_cutoff 100.0
ligand_atom_file 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/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
You can run this job on head node for a short while, just to make sure everything is fine. Then you can kill the job by Ctrl-C. We will use the mpi version to finish this calculation. See next section.
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. (8.virtual_screen_mpi) Copy the input file and the ligand database file from 7.virtual_screen directory.
To submit the job to the seawulf cluster we are using a new file. (virtual.sh)
vi 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 2p16.virtual #PBS -V cd $PBS_O_WORKDIR mpirun -v -np 112 dock6.mpi -i virtual.in -o 2p16.virtual.mpi.out
Before submitting the job, please check the mpirun you are using
which mpirun
make sure you are using an intel mpi.
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 8.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 2nnq.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 ../1.dockprep/2nnq.rec.charged.mol2 cont_score_att_exp 6 cont_score_rep_exp 12 cont_score_rep_rad_scale 1 cont_score_use_dist_dep_dielectric yes cont_score_dielectric 4.0 cont_score_vdw_scale 1 cont_score_es_scale 1 footprint_similarity_score_secondary no pharmacophore_score_secondary no descriptor_score_secondary no gbsa_zou_score_secondary no gbsa_hawkins_score_secondary no SASA_score_secondary no amber_score_secondary no minimize_ligand yes simplex_max_iterations 1000 simplex_tors_premin_iterations 0 simplex_max_cycles 1 simplex_score_converge 0.1 simplex_cycle_converge 1.0 simplex_trans_step 1.0 simplex_rot_step 0.1 simplex_tors_step 10.0 simplex_random_seed 0 simplex_restraint_min no atom_model all vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl ligand_outfile_prefix 2nnq.virtualscreen.minimized write_orientations no num_scored_conformers 1 rank_ligands no
Create a new submission script for the minimization.
touch 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
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 9.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 ../8.cartesianmin/2nnq.virtualscreen.minimized_scored.mol2 limit_max_ligands no skip_molecule no read_mol_solvation no calculate_rmsd no use_database_filter no orient_ligand no bump_filter no score_molecules yes contact_score_primary no contact_score_secondary no grid_score_primary no grid_score_secondary no multigrid_score_primary no multigrid_score_secondary no dock3.5_score_primary no dock3.5_score_secondary no continuous_score_primary no continuous_score_secondary no footprint_similarity_score_primary no footprint_similarity_score_secondary no pharmacophore_score_primary no pharmacophore_score_secondary no descriptor_score_primary yes descriptor_score_secondary no descriptor_use_grid_score no descriptor_use_multigrid_score no descriptor_use_continuous_energy no descriptor_use_footprint_similarity yes descriptor_use_pharmacophore_score yes descriptor_use_tanimoto yes descriptor_use_hungarian yes descriptor_use_volume_overlap yes descriptor_fps_use_footprint_reference_mol2 yes descriptor_fps_footprint_reference_mol2_filename ../4.dock/2nnq.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 ../1.dockprep/2nnq.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 ../4.dock/2nnq.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 ../4.dock/2nnq.lig.min_scored.mol2 descriptor_hungarian_ref_filename ../4.dock/2nnq.lig.min_scored.mol2 descriptor_hungarian_matching_coeff -5 descriptor_hungarian_rmsd_coeff 1 descriptor_volume_reference_mol2_filename ../4.dock/2nnq.lig.min_scored.mol2 descriptor_volume_overlap_compute_method analytical descriptor_weight_fps_score 1 descriptor_weight_pharmacophore_score 1 descriptor_weight_fingerprint_tanimoto -1 descriptor_weight_hungarian 1 descriptor_weight_volume_overlap_score -1 gbsa_zou_score_secondary no gbsa_hawkins_score_secondary no SASA_score_secondary no amber_score_secondary no minimize_ligand no atom_model all vdw_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/vdw_AMBER_parm99.defn flex_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex.defn flex_drive_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/flex_drive.tbl chem_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/chem.defn pharmacophore_defn_file /gpfs/projects/AMS536/zzz.programs/dock6/parameters/ph4.defn ligand_outfile_prefix descriptor.output write_footprints yes write_hbonds yes write_orientations no num_scored_conformers 1 rank_ligands no
