Flexible docking

Setup a directory structure for the docking

Not to end up with a folder cluttered with tons of files a folder structure like the one below this is recommended for a docking project:

~/project_folder/rec
~/project_folder/rec_pcr_00
~/project_folder/lig
~/project_folder/lig_pcr_00
~/project_folder/com
~/project_folder/multidock
~/project_folder/multidock/rec
~/project_folder/multidock/lig
~/project_folder/multidock/com
~/project_folder/multidock/hex

Starting structures

You will need a receptor and a ligand structure and optionally a complex structure. So, the first task is to retrieve receptor and ligand coordinate files of the free (unbound) form (currently only pdb format supported).

Again, just like in the normal docking example, lets say our starting structures are called:

1lig_original.pdb
1rec_original.pdb
1com_original.pdb (optional)

Place the files in the corresponding rec, lig or com folder.

Preparing the structures

Prepare your structures just like in the regular docking howto.

Now run the pdb2xplor.py script for all structures:

>>> cd rec
>>> pdb2xplor.py -i 1rec_original.pdb -exe

>>> cd ../lig
>>> pdb2xplor.py -i 1lig_original.pdb -c 1 -exe

>>> cd ../com
>>> pdb2xplor.py -i 1com_original.pdb

You have now three cleaned and minimized structures with corresponding psf-files:

rec/1REC.pdb  and  1REC.psf
lig/1LIG.pdb  and  1LIG.psf
com/1COM.pdb  and  1COM.psf

Generate a ligand ensemble

This example starts a simulation on localhost without any principal component restraints. The simulation will run 10 MD simulations in parallell (for greater ensemble diversity) and will result in 500 shapshots (separated by 1 ps, 500 x 2 fs time steps).

NOTE: This is a time consuming step! Depending on the size of your system and the hardware you are using the simulation will take anything from days to months.

The pcr force is here set to zero - i.e. this is a standard parallell MD simulation. If you would like an even greater diversity run the script with a -f argument greater than zero (hint: aim at something close to 0.2 times the number of amino acids).

Extracted from the runPCR.py help screen:

Syntax: runPCR -t |psfFolder| -h |host|
               [-f |Force| -r |resultFolder| -n |nice| -i |inpFolder| ]

Options:
       -f      force constant for PCR restraint
       -r      base folder for result (sub-folders will be created)
       -t      folder with topology (psf, pdb)
       -n      nice value
       -h      host computer (accessible via ssh)
       -i      folder with all input files, must contain restart_h2o.inp'
       -parm   folder with param19.* files

Default options:
       - f     0
       - i     ~/biskit/external/xplor/inpPcr
       - h     localhost
       - parm  ~/biskit/external/xplor/toppar
       - n     0
       - r     .
       - nstep         500

Now start the simulation. If you want to run it on another host than the localhost you need pasword free ssh access to that machine.

>>> runPcr.py -t lig -r lig_pcr_00 -h localhost

The script returns information about the simulation that is usefull if you need to stop or restart the simulation:

host: localhost
pid:  14500
nice: 0
binary:        /shared_bin/xplor/ifcxplor
result folder: /mnt/data/shared_bin/biskit/test/lig_pcr_00/pcr_00
topology:      /mnt/data/shared_bin/biskit/test/lig
force:         0

Generate a receptor ensemble

Do the same thing for the receptor as you did for the ligand under Generate a ligand ensemble.

>>> runPcr.py -t rec -r rec_pcr_00 -h localhost

Collect the ligand simulation snapshots

The following step is quite time consuming and has therefore been parallellized through PVM. The first step is thus to start the pvm deamon (if not already running).

>>> pvm

Create a model folder in the lig_pcr_00 folder

>>> cd ~/project_folder/lig_pcr_00
>>> mkdir model

pdb2model.py is a parallellized script that converts all the pdb snapshots into PDBModels. The script reads the pdb files from the lig_pcr_00 folder, strips of any water molecules and writes PDBModels into the model folder.

>>> pdb2model.py -i lig_pcr_00/*.pdb -a -wat -o model/

pdb2traj.py collects the PDBModels created by pdb2model.py and fits them all to a given reference structure, deletes the waters and writes a trajectory output file: traj.dat

>>> pdb2traj.py -i model/*.model -r ../lig/1LIG.pdb -wat -f

If you don't need all the trajectory frames can remove some of them from the trajectort. This will speed up the subsequent calculations in the example below we keep every 5th frame of the trajectory.

>>> thinTraj.py -i traj.dat -step 5

Collect the receptor simulation snapshots

Do the same thing for the receptor as you did for the ligand under Collect the ligand simulation snapshots.

>>> cd ~/project_folder/rec_pcr_00
>>> mkdir model
>>> pdb2model.py -i rec_pcr_00/*.pdb -a -wat -o model/
>>> pdb2traj.py -i model/*.model -r ../lig/1REC.pdb -wat -f

Cluster the trajectories

We don't want to use all 500 snapshots for the docking, rather we want to select representative structures from the ensemble. In this example we select 10 ligand conformers and we do this by means of fuzzy clustering, where each picked structure is a cluster center.

selectModels.py will cluster the trajectory snapshots in ten clusters. The structure closest to the center of each cluster will be written to disc as a PDBModel and a model dictionary.

The -ref option tells the script to add the reference structure of the trajectory being clustered to the model dictionary, thus we will end up with a model dictionary containing 11 structures (10 cluster centers and the minimized structure from Preparing the structures.

>>> cd ~/project_folder/multidock/lig
>>> selectModels.py -i ../../lig_pcr_00/traj.dat -psf ../../lig/1LIG.psf -n 10 -ref

and the same for the receptor:

>>> cd ~/project_folder/multidock/rec
>>> selectModels.py -i ../../rec_pcr_00/traj.dat -psf ../../rec/1REC.psf -n 10 -ref

Prepare the complex structures for docking

This step has already been explained thoroughly in the docking howto.

Use PCR2hex.py to convert the complex pdb files that can be used by HEX.

>>> cd ~/project_folder/multidock/com
>>> PCR2hex.py  -psf ../../com/1COM.psf -pdb ../../com/1COM.pdb

Add additional info to the models

Now well add structure related information to all models that will be used in the docking (acessible surface area, surface curvature, residue conservation etc.)

>>> cd ~/project_folder/multidock/lig
>>> dope.py -s ../../lig/1LIG.pdb  -so ../../lig/1LIG.model -i *.model  -dic 1LIG_model.dic

>>> cd ~/project_folder/multidock/rec
>>> dope.py -s ../../rec/1REC.pdb  -so ../../rec/1REC.model -i *.model  -dic 1REC_model.dic

>>> cd ~/project_folder/multidock/com
>>> dope.py -s ../../com/1COM.pdb  -so ../../com/1COM.model -i *.model  -dic 1COM_model.dic

Run the flexible docking

In this example we are docking 11 receptor structures against 11 ligand structures (i.e. 121 single docking runs) and we are using the known complex for benchmarking. The HEX macro files and the stdout from HEX will be written to disc in the output folder during docking. When the docking is done the macro files and the logfiles are compressed.

This starts the multidocking across the accessible hosts (defined in ~/.biskit/hosts.dat). The output will be a complex list with 11x11x512 docking solutions.

>>> multidock.py -rdic rec/1A2P_model.dic -ldic lig/1A19_models.dic -com com/1BGS_hex.pdb

Calculate complex list data

This step has already been explained thoroughly in the docking howto.

>>> contacter.py -i complexes.cl -a -ref ../../com/ref.complex

Check complex list data completeness

A simple script that will tell you if all the data has been calculated for all the complexes. If not you can use contacter.py script to recalculate the missing data (see the contacter.py help info).

>>> inspectComplexList.py complexes_cont.cl

Analyze docking result

This is just an example of a simple analysis script that prints various information about the quality of the docking results.

>>> a_rmsd_vs_dock_performance.py -cl complexes_cont.cl -ref ../../com/ref.complex -key fnac_10

and another scripts that visualizes the result as contour plots

>>> a_multidock_contour.py -cl complexes_cont.cl