Difference between revisions of "2012 DOCK tutorial with Streptavidin"

From Rizzo_Lab
Jump to: navigation, search
(IX. Frequently Encountered Problems)
(Molecular footprints comparisons)
Line 699: Line 699:
[[File:2012_Dock_Tutorial_Footprints_2.png|thumb|375px|right|2012_DOCK_Tutorial_Footprints_2]] <BR>
==IX. Frequently Encountered Problems==
==IX. Frequently Encountered Problems==

Revision as of 01:23, 11 May 2012

For additional Rizzo Lab tutorials see DOCK Tutorials. Use this link Wiki Markup as a reference for editing the wiki. This tutorial was developed collaboratively by the AMS 536 class of 2012.

I. Introduction





Streptavidin & Biotin



Organizing Directories




II. Preparing the Receptor and Ligand

Downloading the PDB Structure (1DF8)


Preparing the Ligand and Enzyme in Chimera

If you are on "herbie" you can access Chimera directly by typing in Chimera at the prompt. Otherwise, download Chimera to your desktop and obtain your protein/nucleic acid complex of interest using Fetch (by ID).

Open Your Protein of Interest and Delete Unwanted Molecules:



  1. SelectChainBAll.
  2. ActionsAtomsDelete.


Save the file as 1DF8.dockprep.mol2.

dockprep image



Create a Receptor File:

Go to SelectResidueBTN. Then go to ActionsAtomsDelete.

Save the file as 1DF8.receptor.mol2.

Create a Receptor File with No Hydrogen Atoms:

Go to SelectChemistryElementHDelete.

Save a PDB file as 1DF8.receptor.noH.pdb.

Create a Ligand File:

Open only the 1DF8.dockprep.mol2.


Save the file as 1DF8.ligand.mol2.


III. Generating Receptor Surface and Spheres

Receptor Surface








4.0                          #maximum sphere radius in angstroms (default=4.0) 
1DF8.receptor.sph            #clustered spheres file that we want to generate

Save the INSPH file. Then use this command to generate the spheres file:

sphgen -i INSPH -o OUTSPH

-i means input; -o means output.


density type = X
reading  1DF8.receptor.dms                                                     type   R
# of atoms =    881   # of surf pts =  10771
finding spheres for   1DF8.receptor.dms                                         
dotlim =     0.000
radmax =    4.000
Minimum radius of acceptable spheres?
output to  1DF8.receptor.sph                                                    
clustering is complete     28  clusters

You can also open the spheres file that we generated in this step (1DF8_receptor.sph). This file contains detailed information about the spheres, which are divided into 28 clusters. Cluster 0 in the end of the spheres file is a combination of all the clusters.



You will be asked the following questions:

Enter name of sphere cluster file:
Enter name for output file prefix:

In this case, 28 pdb files should be generated.



Then type the command showsphere

showshpere < showsphere.in





Sphere_selector filters the output from sphgen, selecting all spheres within a user-defined radius of a target molecule. In this case, we selected the spheres within 10 angstroms of our ligand. A file called selected_spheres.sph should be created, showing the spheres that are selected. You can again visualize it by showsphere. You can also change the radius (to 8 angstroms or 6 angstroms) or manually edit the file selected_spheres.sph so that you can select the spheres you want. In this tutorial, spheres within 6 angstroms of the ligand are used for the next step (with one sphere which does not belong to the active site being manually deleted).

2012_DOCK_Tutorial_1DF8_surface_spheres_10 angstroms
2012_DOCK_Tutorial_1DF8_surface_spheres_6 angstroms

IV. Generating Box and Grid


In order to speed up docking calculations, DOCK generates a fine grid. At each point in the grid, electrostatic and the VDW probes' energies are precomputed. The energies are computed using a molecular force field. To determine the dimensions of the grid, however, we first generate a box that contains the outer boundaries for the grid calculation. The dimensions and location of the box can be determined using a program called showbox.

First create a directory where you will place the grid files.

$mkdir 03-box-grid
$cd 03-box-grid


The program asks the questions depicted in the diagram the right:

Error creating thumbnail: Unable to save thumbnail to destination

To run the program in the interactive mode, run


To feed the answers to the questions, run

$showbox < showbox.in

For example, showbox.in can contain:







 -i: read from grid.in or input_file, standard_in otherwise
     standard_out otherwise
 -t: terse program output
 -v: verbose program output

For our grid.in file, we will use the following answers:

compute_grids                  yes
grid_spacing                   0.3
output_molecule                no
contact_score                  no
energy_score                   yes
energy_cutoff_distance         9999
atom_model                     a
attractive_exponent            6
repulsive_exponent             9
distance_dielectric            yes
dielectric_factor              4
bump_filter                    yes
bump_overlap                   0.75
receptor_file                  ../01-dockprep/1DF8.receptor.mol2
box_file                       1DF8.box.pdb
score_grid_prefix              grid

Line by line:

  1. compute scoring grids (yes)
  3. write up coordinates of the receptor into a new file
  4. compute contact grid? default is no
  11. distance dielectric factor is the coefficient of the dielectric
  14. our receptor file
  15. the box file we generated in the Box section
  16. VDW parameters file

V. Docking a Single Molecule for Pose Reproduction

Docking and Results


$touch dock.in     #create a file called “dock.in”

Run dock6.4

$dock6 -i dock.in





1st run:


ligand_outfile_prefix [output] ():   		1DF8.output


2nd run:


$vim dock.in

Parameters being changed are listed below:

orient_ligand                                               yes



The receptor_site_file should be the selected spheres file (.sph) generated in a previous step (02 surface and spheres), according to which the program orients the ligand.




Example alt text
DOCK second Run result

3rd Run:

Here we will further specify parameters to improve the result.

calculate_rmsd                                               yes
score_molecules                                              yes
num_scored_conformers                                        10

Again, further questions will be asked when you run the program. Keep most answers as the default values. Paths requiring specification are listed below:



simplex_max_iterations                                        [1000] ():20
write_conformations                                           [yes] (yes no):no
cluster_rmsd_threshold                                        [2.0] ():0.2

The simplex_max_iteration parameter specifies number of minimization cycles.



Example alt text
DOCK third Run result

4th Run:

This time we will set the ligand as flexible:

flexible_ligand                                              yes

Further questions will be asked during the run. Keep most of the answers as the default values except for:

simplex_grow_max_iterations                                   [20] ():500

The simplex_grow_max_iteration specifies the maximum number of iterations per cycle per growth step.


Example alt text
DOCK fourth Run result

5th Run:

This time we will turn the bump filter on:

bump_filter                                                  yes

The bump_filter option filters out conformations that cause clashes between atoms.

Further questions will be asked during the run. Keep answers as the default values. Specify the following path:




Example alt text
DOCK fifth Run result


Following is the final dock.in file used during the fifth run:

ligand_atom_file                                             ../01-dockprep/1DF8.ligand.mol2
limit_max_ligands                                            no
skip_molecule                                                no
read_mol_solvation                                           no
calculate_rmsd                                               yes
use_rmsd_reference_mol                                       no
use_database_filter                                          no
orient_ligand                                                yes
automated_matching                                           yes
receptor_site_file                                           ../02-surface-spheres/selected_spheres_06A.sph
max_orientations                                             500
critical_points                                              no
chemical_matching                                            no
use_ligand_spheres                                           no
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
flexible_ligand                                              yes
min_anchor_size                                              40
pruning_use_clustering                                       yes
pruning_max_orients                                          100
pruning_clustering_cutoff                                    100
pruning_conformer_score_cutoff                               25
use_clash_overlap                                            no
write_growth_tree                                            no
bump_filter                                                  no
score_molecules                                              yes
contact_score_primary                                        no
contact_score_secondary                                      no
grid_score_primary                                           yes
grid_score_secondary                                         no
grid_score_rep_rad_scale                                     1
grid_score_vdw_scale                                         1
grid_score_es_scale                                          1
grid_score_grid_prefix                                       ../03-grid/1DF8.grid
dock3.5_score_secondary                                      no
continuous_score_secondary                                   no
gbsa_zou_score_secondary                                     no
gbsa_hawkins_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                                                ../03-grid/vdw_AMBER_parm99.defn
flex_defn_file                                               ../03-grid/flex.defn
flex_drive_file                                              ../03-grid/flex_drive.tbl
ligand_outfile_prefix                                        1DF8.output
write_orientations                                           no
num_scored_conformers                                        10
write_conformations                                          no
cluster_conformations                                        yes
cluster_rmsd_threshold                                       0.2
rank_ligands                                                 no


$vim sub.dock.csh

wherein you write:

#! /bin/tcsh 
#PBS -l nodes=1:ppn=1 
#PBS -l walltime=01:00:00 
#PBS -o zzz.qsub.out 
#PBS -j oe 

cd ~/AMS536/DOCK_tutorial/05-dock_qsub 



$qsub sub.dock.csh


$chmod +x sub.dock.csh


VI. Virtual Screening

Virtual Screening Preparation



tar -cvf DOCK-Tutorial.tar DOCK-Tutorial/
tar -xvf DOCK-Tutorial.tar DOCK-Tutorial

Virtual Screening Protocol


< ligand_atom_file                                             ligands.mol2 
> ligand_atom_file                                             ../01-dockprep/1DF8.ligand.mol2 
< ligand_outfile_prefix                                        1DF8.vs.output 
> flex_defn_file                                               ../03-box-grid/flex.defn 
> flex_drive_file                                              ../03-box-grid/flex_drive.tbl 
> ligand_outfile_prefix                                        1DF8.output 
< num_scored_conformers                                        1 
> num_scored_conformers                                        10 
< cluster_conformations                                        no 
> cluster_conformations                                        yes 
> cluster_rmsd_threshold                                       0.2 



min_anchor_size: The minimum number of atoms for an anchor.
pruning_use_clustering: Pruning conformations during the clustering process.
use_internal_energy: Using repulsive VDM to avoid internal clashes.


#! /bin/tcsh
#PBS -l nodes=4:ppn=2
#PBS -l walltime=01:00:00
#PBS -o zzz.qsub.out
#PBS -j oe
echo "Running on ${nprocs} processors"
echo ""
echo "processor list are:"
cd ~/AMS536/DOCK_Tutorial/06-virtualscreen

Finally, we can submit the job and perform virtual screening.

qsub sub.virtual_screen.csh

Virtual Screening Results



Molecule: ZINC33171556

Anchors:               1
Orientations:          500
Conformations:         116

   Grid Score:          -52.373383
     Grid_vdw:          -51.643253
      Grid_es:           -0.730129
   Int_energy:            6.125062


Molecule: ZINC20605433

Elapsed time:  0 seconds

ERROR:  Could not complete growth.
        and try increasing max_orientations.

If you download the multi-mol2 file from seawulf using:



Example of Virtual Screening Result

VII. Running DOCK in Serial and in Parallel on Seawulf



Running DOCK in Serial on a Single Processor


#PBS -l nodes=1:ppn=1        
#PBS -l walltime=01:00:00   
#PBS -N dock6           
#PBS -M user@ic.sunysb.edu 
#PBS -j oe                   
#PBS -o pbs.out            

cd /nfs/user03/sudipto/DOCK_Tutorial                     

Here is an explanation of the code and format:

#!/bin/tcsh                                                 #Execute script with tcsh
#PBS -N dock6                                               #Name of your job
#PBS -o pbs.out                                             #Name of your output file


Running DOCK in Parallel using MPI


#! /bin/tcsh
#PBS -l nodes=4:ppn=2
#PBS -l walltime=24:00:00
#PBS -o zzz.qsub.out
#PBS -j oe


echo "Running on ${nprocs} processors"

cd /nfs/user03/sudipto/DOCK_Tutorial


For more information, see PBS Queue

Serial Calculation for Pose Reproduction

Parallel Virtual Screen

VIII. Molecular Footprints

Preparation (minimization) of Reference Molecule

The footprint of the reference molecule guides the selection of molecules/poses in the footprint rescoring method, therefore careful preparation of this reference molecule is necessary. Typically the reference molecule will be an endogenous substrate or known inhibitor from an x-ray crystallographic structure.

Preferably, the molecule should be prepared with hydrogens and charges using AMBER and subsequently converted to mol2 format using the 'top2mol2' command. This will provide you with a 'monotonic' mol2 file -- meaning that all the hydrogens will be listed with the corresponding protein residue, rather than separately at the end of the file (like Chimera's 'DockPrep' tool does).

Once the molecule is prepared, it should be minimized in the receptor using a tethered minimization. This is especially important if the crystal structure was refined without hydrogens (typically hydrogens are not used in refinement unless the resolution is ~1 Angstrom or better). Minimization is needed to help alleviate any steric clashes involving the newly added hydrogens (which would otherwise disrupt your footprint results).

An example input file for use on seawulf (using dock6.5):

ligand_atom_file                                             ligand.mol2
limit_max_ligands                                            no
skip_molecule                                                no
read_mol_solvation                                           no
calculate_rmsd                                               yes
use_rmsd_reference_mol                                       no
use_database_filter                                          no
orient_ligand                                                no
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
flexible_ligand                                              no
bump_filter                                                  no
score_molecules                                              yes
contact_score_primary                                        no
contact_score_secondary                                      no
grid_score_primary                                           no
grid_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                                      receptor.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
descriptor_score_secondary                                   no
gbsa_zou_score_secondary                                     no
gbsa_hawkins_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                                                /nfs/user03/tbalius/zzz.programs/dock6.5/parameters/vdw_AMBER_parm99.defn
flex_defn_file                                               /nfs/user03/tbalius/zzz.programs/dock6.5/parameters/flex.defn
flex_drive_file                                              /nfs/user03/tbalius/zzz.programs/dock6.5/parameters/flex_drive.tbl
ligand_outfile_prefix                                        xtalmin
write_orientations                                           no
num_scored_conformers                                        1
rank_ligands                                                 no

Molecular footprints comparisons

After rescoring of the molecules by footprint-based score, typically three output files will be generated: a mol2 file showing the ligands information, a footprint.txt file showing the footprints infomation, and a hbond.txt file showing the hydrogen bond information. The footprint.txt file contain per residue interaction of the receptor with the ligand, for both reference molecule and the docked molecule. Shown below is a VDW footprints comparison.


IX. Frequently Encountered Problems


How to Open Files

It is important to remember that if you save files to your desktop as mol2 files or pdb files when preparing your ligand and receptor for docking, they cannot be opened using Windows because Windows only opens pdf files using Adobe Reader. You can only access these files through Chimera.

If you want to view a file in chimera, first transfer the file from "herbie" or "seawulf" to your laptop using WinSCP. Then open Chimera on your laptop to open the files.

Woo Suk

Distinguishing overlapped spheres





Installation of DOCK

If you want to install DOCK on your own desktop or laptop, and you are a beginner of Unix/Linux, you might encounter some problems during installation. A good reference for the installation is the DOCK User Manual(http://dock.compbio.ucsf.edu/DOCK_6/dock6_manual.htm). The problem I encountered was: after I used gnu as configuration file to configure the Makefile, I successfully made the config.h file, but when I was trying to build the DOCK executables, an error message " g77 command not found " appeared. And it is quite inconvenient to install the g77 command. A way to solve this problem is to open the config.h file, and manually change the g77 in that file to gfortran.

Use the Correct Path of Your File

One frequently encountered problem is using the wrong path of your file. For example, when you use sphere_selector to generate spheres around the ligand, it requires two files--one is the spheres file .sph and another is your ligand file, but your ligand file is probably not in your working directory! So you need either copy your ligand file to the current directory or specify the correct path of your ligand file like below:


This is also important when you run DOCK. Remember to use the correct path and name of the your file.

Notice the Shell That You are Using

This tutorial is based on the shell tcsh, and maybe you are using a different shell. This will not cause any trouble in most of the time, but you may want to notice this in some circumstances, for example when you want to execute some .csh file. One way to solve this problem is that simply type "tcsh" in the command line, and you will be running tcsh in another shell.


One of the problems one may encounter with grid creation is the sampling location for ligand binding. It is important to create a grid and dock the ligand to the correct location in the protein. If the grid is created in the incorrect location, the ligand binding sites will be sampled within the grid, and the best binding site will not be sampled. In our case, we are docking it to an area where a ligand is known to bind to streptavidin. However, this may not always be the case.

One solution is to dock to the center of mass of the spheres clusters created with spheregen and assess the center of mass as viable or not based on your chemical intuition.

Another way is to assess the protein for buried sites where there may also be a number of spheres. Both are options in the showbox program.

Finally, one can attempt to find proteins with similar sequence and determine where their ligand binding sites are, and then pick the binding site in the protein you are docking to using that information.



Keep the previous dock6 output .mol2 files



$mv 1DF8.output_scored.mol2 first_1DF8.output_scored.mol2


Some recommended parameters we can change when running dock

max_orientations                                             500
min_anchor_size                                              40
pruning_conformer_score_cutoff                               25
pruning_clustering_cutoff                                    100

Explanations of these parameters can be found on the website provided by Long Fei. Increase the number of these parameters can lead to increased running time. But it enables more thorough cycling and may give you better results.


Always Remember Your Path

One situation often appears is that the program will return an error saying "Could not open 1DF8.ligands.mol2 for reading." or something like that indicating that the program couldn't find the input file you specified in you dock.in or vs.in or some other files. In this case, the most often reason is that the path you specified in your input file is wrong. So the first thing you might want to do is to check those paths, pay special attention to those paths where you use "." to indicate the same folder and ".." to indicate the folder up one level because this is the place where you are most likely to make mistakes. If you are not sure of the usage of "." and "..", it will be safer for you to use the absolute paths of the files. And you can use the command pwd to ask the system where you are if you don't know that.


Delete Spheres Manually


1. Transfer sph file to pbd file by using showsphere.

DOCK spheres within 6.0 ang of ligands
cluster     1   number of spheres in cluster    22
  60  28.86000  11.09173   4.97943   2.337  586 0  0
  67  28.67840   9.74620  12.13940   1.400   60 0  0
 161  27.11509  13.42709  13.48051   1.401  385 0  0
 174  27.79940  12.82468  12.04078   2.475  480 0  0
 183  34.13903  14.01393   7.31591   3.590  691 0  0
 187  36.78641  17.02434   9.54436   3.481  691 0  0
 385  30.01753  14.79655  14.46633   1.402  174 0  0
 463  25.74631  14.83304   9.98248   1.400  589 0  0


Molecule preparation for molecular footprinting

A sample script:

#check programs
foreach prog (antechamber parmchk tleap)
        set prog_version = `which $prog`
        if ( $prog_version == /nfs/user03/tbalius/amber10_ompi/bin/$prog ) then
                echo $prog OK, current path is: `which $prog` 
                echo please double-check $prog path, current path is: `which $prog`

#set your home directory on seawulf
set mydir = "$HOME"

# set up your ligands
foreach lig (ligand_name1 ligand_name2 ligand_name3)

        #set working dir and go there
        echo setting up directory for: $lig
        mkdir  "$mydir/DOCK_project2/10-ligandprep-top2mol2/$lig"
        set workdir = "$mydir/DOCK_project2/10-ligandprep-top2mol2/$lig"
        cd $workdir
        cp ../$lig.mol2 .
        echo 'Done!'

        #convert mol2 to pdb and prep files
        echo converting mol2 to pdb and prep files
        antechamber -i $lig.mol2 -fi mol2  -o $lig.ante.pdb  -fo pdb
        antechamber -i $lig.mol2 -fi mol2  -o $lig.ante.prep -fo prepi
        echo 'Done!'

        #check for missing parameters
        echo checking parameters for ligand: $lig
        parmchk -i $lig.ante.prep -f  prepi -o $lig.ante.frcmod
        echo 'Done!' 

        #kludgey way to make input files
        echo making input file for ligand: $lig
        echo "set default PBradii mbondi2" > tleap.$lig.in                                          
        echo "source leaprc.ff99SB" >> tleap.$lig.in        
        echo "loadoff ../rizzo_amber7.ionparms/ions.lib" >> tleap.$lig.in
        echo "loadamberparams ../rizzo_amber7.ionparms/ions.frcmod" >> tleap.$lig.in             
        echo "source leaprc.gaff" >> tleap.$lig.in
        echo "loadamberparams ../rizzo_amber7.ionparms/gaff.dat.rizzo" >> tleap.$lig.in
        echo "loadamberparams $lig.ante.frcmod" >> tleap.$lig.in
        echo "loadamberprep $lig.ante.prep" >> tleap.$lig.in     
        echo "lig = loadpdb $lig.ante.pdb" >> tleap.$lig.in      
        echo "saveamberparm lig $lig.gas.leap.prm7 $lig.gas.leap.rst7" >> tleap.$lig.in  
        echo "solvateBox lig TIP3PBOX 10.0" >> tleap.$lig.in                           
        echo "saveamberparm lig $lig.wat.leap.prm7 $lig.wat.leap.rst7" >> tleap.$lig.in  
        echo "charge lig" >> tleap.$lig.in                                                      
        echo "quit" >> tleap.$lig.in
        echo 'Done!' 

        #run tleap
        echo running tleap for ligand: $lig
        tleap -s -f tleap.$lig.in > tleap.$lig.out
        echo 'Done!' 

        #make a new mol2 file
        echo creating final mol2 file for ligand: $lig
        top2mol2 -p $lig.gas.leap.prm7 -c $lig.gas.leap.rst7 -o $lig.gas.leap.mol2
        echo 'Done!'

        #make a new pdb file
        echo creating final pdb file for ligand: $lig
        ambpdb -p $lig.gas.leap.prm7 -tit "pdb" < $lig.gas.leap.rst7 > $lig.gas.leap.pdb
        echo 'Done!'

Keeping your data in sync






rsync -arv source target

Deleting jobs


qstat -u username


qdel jobid