This is for ligands that modify a specific residue.

(1) Make a directory:
mkdir coval_min
cd coval_min/

(2) Download pdb file from the web:
wget https://files.rcsb.org/view/5YY1.pdb

(3) Manually create 3 files using your favorate text editor: a covalent ligand, a co-facotor and a receptor file.

(3.1) for the receptor, rec.pdb:

Change the covalent cysteine to reduce name CYM.

Here is a modified cystiene residue (this is without a charge and without a hydrogen):

cat CYM.prep

0 0 2

CYSTEINE without h and without charge for covalent

CYM INT 0
CORR OMIT DU BEG
0.00000
1 DUMM DU M 0 -1 -2 0.000 0.000 0.000 0.00000
2 DUMM DU M 1 0 -1 1.449 0.000 0.000 0.00000
3 DUMM DU M 2 1 0 1.522 111.100 0.000 0.00000
4 N N M 3 2 1 1.335 116.600 180.000 -0.41570
5 H H E 4 3 2 1.010 119.800 0.000 0.27190
6 CA CX M 4 3 2 1.449 121.900 180.000 0.02130
7 HA H1 E 6 4 3 1.090 109.500 300.000 0.11240
8 CB 2C 3 6 4 3 1.525 111.100 60.000 -0.12310
9 HB2 H1 E 8 6 4 1.090 109.500 300.000 0.11120
10 HB3 H1 E 8 6 4 1.090 109.500 60.000 0.11120
11 SG SH E 8 6 4 1.810 116.000 180.000 -0.23580
12 C C M 6 4 3 1.522 111.100 180.000 0.59730
13 O O E 12 6 4 1.229 120.500 0.000 -0.56790

IMPROPER
-M CA N H
CA +M C O

DONE
STOP

Here is a frcmod file with need parameters:
cat thioe.frcmod

Nir and trent got the parms from gaff and change the names, for thio-ether
MASS

BOND
c3-SH 225.8 1.8210 SOURCE1 358 0.0075
C -SH 225.8 1.8210 SOURCE1 358 0.0075

ANGLE
2C-SH-c3 60.63 99.92 SOURCE3 14 2.0723
CT-SH-c3 60.63 99.92 SOURCE3 14 2.0723
c3-c3-SH 61.10 112.69 SOURCE3 24 2.1842
hc-c3-SH 42.51 108.76 SOURCE2 3 1.6891
SH-c3-h1 42.51 108.76 SOURCE2 3 1.6891

DIHE
X -c3-SH-X 3 1.000 0.000 3.000 JCC,7,(1986),230

IMPROPER

NONBON

These parameter values were taken from the gaff force field (/nfs/soft/amber/amber14/dat/leap/parm/gaff.dat).

(3.2) for a covalent ligand, lig.pdb:

keep the sidechain of the mofifed cystien:

cat 94F_bkup.pdb

ATOM 89 CA CYS A 12 29.215 -1.013 17.747 1.00 23.10 C
ATOM 92 CB CYS A 12 30.264 -1.172 18.849 1.00 23.68 C
ATOM 93 SG CYS A 12 29.632 -1.020 20.503 1.00 25.01 S
ATOM 1385 C31 94F A 203 16.511 -4.125 24.817 1.00 32.17 C
ATOM 1386 C30 94F A 203 15.606 -4.260 23.787 1.00 32.36 C
ATOM 1387 C32 94F A 203 17.840 -3.910 24.531 1.00 31.54 C
ATOM 1388 C29 94F A 203 16.037 -4.183 22.484 1.00 32.47 C
ATOM 1389 C19 94F A 203 21.614 -2.199 22.881 1.00 29.12 C
ATOM 1390 C5 94F A 203 23.964 -5.464 22.540 1.00 29.57 C
ATOM 1391 C18 94F A 203 22.453 -3.293 22.771 1.00 29.13 C
ATOM 1392 C23 94F A 203 18.284 -3.829 23.223 1.00 31.43 C
ATOM 1393 C22 94F A 203 19.692 -3.612 23.019 1.00 29.79 C
ATOM 1394 C24 94F A 203 17.372 -3.972 22.203 1.00 32.05 C
ATOM 1395 C3 94F A 203 21.900 -4.565 22.770 1.00 29.39 C
ATOM 1396 C2 94F A 203 20.533 -4.693 22.895 1.00 29.45 C
ATOM 1397 C20 94F A 203 20.247 -2.347 23.012 1.00 29.51 C
ATOM 1398 C7 94F A 203 23.836 -3.189 22.644 1.00 28.94 C
ATOM 1399 C14 94F A 203 29.536 0.583 20.690 1.00 25.47 C
ATOM 1400 C13 94F A 203 28.376 1.143 21.474 1.00 26.40 C
ATOM 1401 C12 94F A 203 27.093 0.850 20.795 1.00 26.10 C
ATOM 1402 C9 94F A 203 23.979 -0.953 21.731 1.00 28.62 C
ATOM 1403 C17 94F A 203 25.969 -2.149 22.595 1.00 29.17 C
ATOM 1404 C10 94F A 203 25.024 -0.390 20.794 1.00 28.12 C
ATOM 1405 C16 94F A 203 26.678 -0.809 22.586 1.00 28.22 C
ATOM 1406 C25 94F A 203 17.826 -3.897 20.794 1.00 32.91 C
ATOM 1407 N4 94F A 203 22.652 -5.665 22.653 1.00 29.29 N
ATOM 1408 N6 94F A 203 24.595 -4.289 22.534 1.00 29.70 N
ATOM 1409 N8 94F A 203 24.523 -1.948 22.642 1.00 29.33 N
ATOM 1410 N11 94F A 203 26.238 0.005 21.476 1.00 27.18 N
ATOM 1411 O15 94F A 203 26.847 1.362 19.713 1.00 25.34 O
ATOM 1412 F1 94F A 203 19.995 -5.935 22.878 1.00 29.69 F
ATOM 1413 F26 94F A 203 16.791 -3.770 19.931 1.00 34.77 F
ATOM 1414 F27 94F A 203 18.674 -2.862 20.538 1.00 31.92 F
ATOM 1415 F28 94F A 203 18.495 -5.026 20.445 1.00 33.49 F
ATOM 1416 CL 94F A 203 19.224 -0.960 23.165 1.00 29.22 CL

change and make the numbering and naming consistent:

cat 94F.pdb

ATOM 1382 CA LIG A 203 29.215 -1.013 17.747 1.00 23.10 C
ATOM 1383 CB LIG A 203 30.264 -1.172 18.849 1.00 23.68 C
ATOM 1384 SG LIG A 203 29.632 -1.020 20.503 1.00 25.01 S
ATOM 1385 C31 LIG A 203 16.511 -4.125 24.817 1.00 32.17 C
ATOM 1386 C30 LIG A 203 15.606 -4.260 23.787 1.00 32.36 C
ATOM 1387 C32 LIG A 203 17.840 -3.910 24.531 1.00 31.54 C
ATOM 1388 C29 LIG A 203 16.037 -4.183 22.484 1.00 32.47 C
ATOM 1389 C19 LIG A 203 21.614 -2.199 22.881 1.00 29.12 C
ATOM 1390 C5 LIG A 203 23.964 -5.464 22.540 1.00 29.57 C
ATOM 1391 C18 LIG A 203 22.453 -3.293 22.771 1.00 29.13 C
ATOM 1392 C23 LIG A 203 18.284 -3.829 23.223 1.00 31.43 C
ATOM 1393 C22 LIG A 203 19.692 -3.612 23.019 1.00 29.79 C
ATOM 1394 C24 LIG A 203 17.372 -3.972 22.203 1.00 32.05 C
ATOM 1395 C3 LIG A 203 21.900 -4.565 22.770 1.00 29.39 C
ATOM 1396 C2 LIG A 203 20.533 -4.693 22.895 1.00 29.45 C
ATOM 1397 C20 LIG A 203 20.247 -2.347 23.012 1.00 29.51 C
ATOM 1398 C7 LIG A 203 23.836 -3.189 22.644 1.00 28.94 C
ATOM 1399 C14 LIG A 203 29.536 0.583 20.690 1.00 25.47 C
ATOM 1400 C13 LIG A 203 28.376 1.143 21.474 1.00 26.40 C
ATOM 1401 C12 LIG A 203 27.093 0.850 20.795 1.00 26.10 C
ATOM 1402 C9 LIG A 203 23.979 -0.953 21.731 1.00 28.62 C
ATOM 1403 C17 LIG A 203 25.969 -2.149 22.595 1.00 29.17 C
ATOM 1404 C10 LIG A 203 25.024 -0.390 20.794 1.00 28.12 C
ATOM 1405 C16 LIG A 203 26.678 -0.809 22.586 1.00 28.22 C
ATOM 1406 C25 LIG A 203 17.826 -3.897 20.794 1.00 32.91 C
ATOM 1407 N4 LIG A 203 22.652 -5.665 22.653 1.00 29.29 N
ATOM 1408 N6 LIG A 203 24.595 -4.289 22.534 1.00 29.70 N
ATOM 1409 N8 LIG A 203 24.523 -1.948 22.642 1.00 29.33 N
ATOM 1410 N11 LIG A 203 26.238 0.005 21.476 1.00 27.18 N
ATOM 1411 O15 LIG A 203 26.847 1.362 19.713 1.00 25.34 O
ATOM 1412 F1 LIG A 203 19.995 -5.935 22.878 1.00 29.69 F
ATOM 1413 F26 LIG A 203 16.791 -3.770 19.931 1.00 34.77 F
ATOM 1414 F27 LIG A 203 18.674 -2.862 20.538 1.00 31.92 F
ATOM 1415 F28 LIG A 203 18.495 -5.026 20.445 1.00 33.49 F
ATOM 1416 CL LIG A 203 19.224 -0.960 23.165 1.00 29.22 CL

Uses chimera to add hydrogens:

[[File:2018_12_coval1.png|thumb|center|500px|Use chimera to add hydrogens]]

Save the protonated molecule in pdb format and rename the atoms:

cat 94F.full_mod_num.pdb

HETATM 1 C01 LIG A 1 29.215 -1.013 17.747 1.00 0.00 C
HETATM 2 C02 LIG A 1 30.264 -1.172 18.849 1.00 0.00 C
HETATM 3 S03 LIG A 1 29.632 -1.020 20.503 1.00 0.00 S
HETATM 4 C04 LIG A 1 16.511 -4.125 24.817 1.00 0.00 C
HETATM 5 C05 LIG A 1 15.606 -4.260 23.787 1.00 0.00 C
HETATM 6 C06 LIG A 1 17.840 -3.910 24.531 1.00 0.00 C
HETATM 7 C07 LIG A 1 16.037 -4.183 22.484 1.00 0.00 C
HETATM 8 C08 LIG A 1 21.614 -2.199 22.881 1.00 0.00 C
HETATM 9 C09 LIG A 1 23.964 -5.464 22.540 1.00 0.00 C
HETATM 10 C10 LIG A 1 22.453 -3.293 22.771 1.00 0.00 C
HETATM 11 C11 LIG A 1 18.284 -3.829 23.223 1.00 0.00 C
HETATM 12 C12 LIG A 1 19.692 -3.612 23.019 1.00 0.00 C
HETATM 13 C13 LIG A 1 17.372 -3.972 22.203 1.00 0.00 C
HETATM 14 C14 LIG A 1 21.900 -4.565 22.770 1.00 0.00 C
HETATM 15 C15 LIG A 1 20.533 -4.693 22.895 1.00 0.00 C
HETATM 16 C16 LIG A 1 20.247 -2.347 23.012 1.00 0.00 C
HETATM 17 C17 LIG A 1 23.836 -3.189 22.644 1.00 0.00 C
HETATM 18 C18 LIG A 1 29.536 0.583 20.690 1.00 0.00 C
HETATM 19 C19 LIG A 1 28.376 1.143 21.474 1.00 0.00 C
HETATM 20 C20 LIG A 1 27.093 0.850 20.795 1.00 0.00 C
HETATM 21 C21 LIG A 1 23.979 -0.953 21.731 1.00 0.00 C
HETATM 22 C22 LIG A 1 25.969 -2.149 22.595 1.00 0.00 C
HETATM 23 C23 LIG A 1 25.024 -0.390 20.794 1.00 0.00 C
HETATM 24 C24 LIG A 1 26.678 -0.809 22.586 1.00 0.00 C
HETATM 25 C25 LIG A 1 17.826 -3.897 20.794 1.00 0.00 C
HETATM 26 N26 LIG A 1 22.652 -5.665 22.653 1.00 0.00 N
HETATM 27 N27 LIG A 1 24.595 -4.289 22.534 1.00 0.00 N
HETATM 28 N28 LIG A 1 24.523 -1.948 22.642 1.00 0.00 N
HETATM 29 N29 LIG A 1 26.238 0.005 21.476 1.00 0.00 N
HETATM 30 O30 LIG A 1 26.847 1.362 19.713 1.00 0.00 O
HETATM 31 F31 LIG A 1 19.995 -5.935 22.878 1.00 0.00 F
HETATM 32 F32 LIG A 1 16.791 -3.770 19.931 1.00 0.00 F
HETATM 33 F33 LIG A 1 18.674 -2.862 20.538 1.00 0.00 F
HETATM 34 F34 LIG A 1 18.495 -5.026 20.445 1.00 0.00 F
HETATM 35 CL LIG A 1 19.224 -0.960 23.165 1.00 0.00 Cl
HETATM 36 H36 LIG A 1 31.026 -0.406 18.705 1.00 0.00 H
HETATM 37 H37 LIG A 1 30.732 -2.151 18.745 1.00 0.00 H
HETATM 38 H38 LIG A 1 16.180 -4.188 25.843 1.00 0.00 H
HETATM 39 H39 LIG A 1 14.561 -4.426 24.003 1.00 0.00 H
HETATM 40 H40 LIG A 1 18.547 -3.803 25.340 1.00 0.00 H
HETATM 41 H41 LIG A 1 15.328 -4.288 21.676 1.00 0.00 H
HETATM 42 H42 LIG A 1 22.038 -1.206 22.864 1.00 0.00 H
HETATM 43 H43 LIG A 1 24.581 -6.345 22.443 1.00 0.00 H
HETATM 44 H44 LIG A 1 29.494 1.022 19.693 1.00 0.00 H
HETATM 45 H45 LIG A 1 30.457 0.921 21.166 1.00 0.00 H
HETATM 46 H46 LIG A 1 28.366 0.694 22.467 1.00 0.00 H
HETATM 47 H47 LIG A 1 28.494 2.222 21.569 1.00 0.00 H
HETATM 48 H48 LIG A 1 23.559 -0.136 22.317 1.00 0.00 H
HETATM 49 H49 LIG A 1 23.184 -1.410 21.141 1.00 0.00 H
HETATM 50 H50 LIG A 1 26.283 -2.718 23.470 1.00 0.00 H
HETATM 51 H51 LIG A 1 26.228 -2.702 21.692 1.00 0.00 H
HETATM 52 H52 LIG A 1 25.273 -1.150 20.054 1.00 0.00 H
HETATM 53 H53 LIG A 1 24.608 0.478 20.282 1.00 0.00 H
HETATM 54 H54 LIG A 1 27.752 -0.975 22.502 1.00 0.00 H
HETATM 55 H55 LIG A 1 26.468 -0.286 23.519 1.00 0.00 H
HETATM 56 H56 LIG A 1 29.692 -1.121 16.773 1.00 0.00 H
HETATM 57 H57 LIG A 1 28.758 -0.026 17.820 1.00 0.00 H
HETATM 58 H58 LIG A 1 28.448 -1.778 17.863 1.00 0.00 H
CONECT 23 21 29 52 53
CONECT 20 19 29 30
CONECT 19 18 20 46 47
CONECT 18 3 19 44 45
CONECT 24 22 29 54 55
CONECT 22 24 28 50 51
CONECT 10 8 14 17
CONECT 8 10 16 42
CONECT 15 12 14 31
CONECT 16 8 12 35
CONECT 12 11 15 16
CONECT 11 6 12 13
CONECT 13 7 11 25
CONECT 25 13 32 33 34
CONECT 7 5 13 41
CONECT 14 10 15 26
CONECT 5 7 4 39
CONECT 4 6 5 38
CONECT 6 4 11 40
CONECT 9 26 27 43
CONECT 17 10 27 28
CONECT 21 23 28 48 49
CONECT 1 2 57 56 58
CONECT 2 1 3 36 37
CONECT 35 16
CONECT 31 15
CONECT 32 25
CONECT 33 25
CONECT 34 25
CONECT 52 23
CONECT 53 23
CONECT 46 19
CONECT 47 19
CONECT 44 18
CONECT 45 18
CONECT 54 24
CONECT 55 24
CONECT 50 22
CONECT 51 22
CONECT 42 8
CONECT 41 7
CONECT 39 5
CONECT 38 4
CONECT 40 6
CONECT 43 9
CONECT 48 21
CONECT 49 21
CONECT 56 1
CONECT 57 1
CONECT 58 1
CONECT 36 2
CONECT 37 2
CONECT 29 20 23 24
CONECT 26 9 14
CONECT 27 9 17
CONECT 28 17 21 22
CONECT 30 20
CONECT 3 2 18
END

run the following script:

cat run.002.ligprep.antechamber.csh

#! /bin/tcsh

set workdir = `pwd`
cd $workdir

setenv AMBERHOME /nfs/soft/amber/amber14

rm lig; mkdir lig; cd lig

cp $workdir/94F.full_mod_num.pdb lig.pdb

$AMBERHOME/bin/antechamber -i lig.pdb -fi pdb -o lig.ante.mol2 -fo mol2

$AMBERHOME/bin/antechamber -i lig.ante.mol2 -fi mol2 -o lig.ante.charge.mol2 -fo mol2 -c bcc -at sybyl -nc 0
$AMBERHOME/bin/antechamber -i lig.ante.mol2 -fi mol2 -o lig.ante.pdb -fo pdb
$AMBERHOME/bin/antechamber -i lig.ante.charge.mol2 -fi mol2 -o lig.ante.charge.prep -fo prepi
$AMBERHOME/bin/parmchk -i lig.ante.charge.prep -f prepi -o lig.ante.charge.frcmod

modify the prep file to remove the cysteine side chain. Also add partial charge to covalent carbon to make the Cysteine (CYM) + covalent ligand an integer.

diff lig/lig.ante.charge.prep lig/lig.ante.charge.mod.prep

11,19c11
< 4 C01 c3 M 3 2 1 1.540 111.208 -180.000 -0.104100
< 5 H56 hc E 4 3 2 1.090 115.954 83.673 0.047367
< 6 H57 hc E 4 3 2 1.090 65.214 -175.770 0.047367
< 7 H58 hc E 4 3 2 1.089 45.251 -8.622 0.047367
< 8 C02 c3 M 4 3 2 1.530 133.294 -81.790 -0.002300
< 9 H36 h1 E 8 4 3 1.090 108.130 -134.600 0.074200
< 10 H37 h1 E 8 4 3 1.090 108.604 107.772 0.074200
< 11 S03 ss M 8 4 3 1.777 114.690 -13.854 -0.332200
< 12 C18 c3 M 11 8 4 1.617 102.332 -80.293 -0.008300
---
> 12 C18 c3 M 3 2 1 1.617 102.332 -80.293 -0.037200

(3.3) for a co-factory (gdp), gdp.pdb:

Add hydrogens to gdp is chimera.

build the prep and frcmod files with antechamber:

cat run.002.ligprep.antechamber_gdp.csh

#! /bin/tcsh

set workdir = `pwd`
cd $workdir

setenv AMBERHOME /nfs/soft/amber/amber14

rm lig2; mkdir lig2; cd lig2

#cp $workdir/xtal-lig.pdb lig.pdb
#cp $workdir/33443.pdb lig.pdb
cp $workdir/gdp_h.pdb lig.pdb
#sed -i 's/<0> /LIG/g' lig1.mol2

$AMBERHOME/bin/antechamber -i lig.pdb -fi pdb -o lig.ante.mol2 -fo mol2

$AMBERHOME/bin/antechamber -i lig.ante.mol2 -fi mol2 -o lig.ante.charge.mol2 -fo mol2 -c bcc -at sybyl -nc -3
$AMBERHOME/bin/antechamber -i lig.ante.mol2 -fi mol2 -o lig.ante.pdb -fo pdb
$AMBERHOME/bin/antechamber -i lig.ante.charge.mol2 -fi mol2 -o lig.ante.charge.prep -fo prepi
$AMBERHOME/bin/parmchk -i lig.ante.charge.prep -f prepi -o lig.ante.charge.frcmod

Here is the tleap input file:

cat tleap.in | awk '{print " "$0}'

set default PBradii mbondi2
# load the protein force field
source leaprc.ff12SB
# load in GAFF
source leaprc.gaff
# ions
loadamberparams /nfs/soft/amber/amber14/dat/leap/parm/frcmod.ionsjc_tip3p
loadamberparams /nfs/soft/amber/amber14/dat/leap/parm/frcmod.ionslrcm_hfe_tip3p

# load ligand and covalent parameters.
loadamberparams lig/lig.ante.charge.frcmod
loadamberparams lig2/lig.ante.charge.frcmod

loadamberparams thioe.frcmod

loadamberprep lig/lig.ante.charge.mod.prep

loadamberprep lig2/lig.ante.charge.prep

loadamberprep CYM.prep

# load pdb file
REC = loadpdb rec.pdb
LIG = loadpdb 94F.full_mod_num.pdb
COF = loadpdb gdp.pdb
#complex
COM = combine {REC COF LIG}

# draw bond between CYN and LIG
#bondByDistance COM 2.0
bond COM.172.C18 COM.12.SG
#bond COM.330.C COM.331.N
#bond COM.331.C COM.332.N
#deleteBond COM.331.SG COM.331.C

#desc COM.331
#desc COM.331.SG
#desc COM.331.C
#desc COM.331.N
#desc COM.501
desc COM.172.C18
desc COM.12.SG

saveamberparm COM com.leap.prm7 com.leap.rst7

quit

Here is the command to run leap:
$AMBERHOME/bin/tleap -s -f tleap.in > ! tleap.out

Here is the command to submit the miminization to the queue:
cat run.004.pmemd_cuda_min.csh | awk '{print " "$0}'

#setenv AMBERHOME /nfs/soft/amber/amber14

setenv AMBERHOME /nfs/soft/amber/amber14
setenv LD_LIBRARY_PATH ""
#setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"
setenv LD_LIBRARY_PATH "/nfs/soft/cuda-6.5/lib64/:\$LD_LIBRARY_PATH"

cat << EOF1 > ! 01mi.in
01mi.in: minimization with GB
&cntrl
imin = 1, maxcyc = 10000, ncyc = 500, ntmin = 1,
igb=1,
ntx = 1, ntc = 1, ntf = 1,
ntb = 0, ntp = 0,
ntwx = 1000, ntwe = 0, ntpr = 1000,
cut = 999.9,
ntr = 1,
restraintmask = '!@H=',
restraint_wt = 0.1,
/
EOF1

#$AMBERHOME/bin/pmemd.cuda -O -i 01mi.in -o 01mi.out -p com.leap.prm7 -c com.leap.rst7 -ref com.leap.rst7 -x 01mi.mdcrd -inf 01mi.info -r 01mi.rst7
#$AMBERHOME/bin/sander -O -i 01mi.in -o 01mi.out -p com.leap.prm7 -c com.leap.rst7 -ref com.leap.rst7 -x 01mi.mdcrd -inf 01mi.info -r 01mi.rst7

set pwd = `pwd`
#cd $pwd

cat << EOF > ! qsub.sander.csh
#\$ -S /bin/csh
#\$ -cwd
#\$ -q gpu.q
#\$ -o stdout
#\$ -e stderr

cd $pwd

$AMBERHOME/bin/pmemd.cuda -O -i 01mi.in -o 01mi.out -p com.leap.prm7 -c com.leap.rst7 -ref com.leap.rst7 -x 01mi.mdcrd -inf 01mi.info -r 01mi.rst7

EOF

qsub qsub.sander.csh

Visualize the minimized coordinates with Chimera or VMD, you can first convert the files to pdb format with the following commands:

$AMBERHOME/bin/ambpdb -p com.leap.prm7 < 01mi.rst7 > 01mi.pdb

$AMBERHOME/bin/ambpdb -p com.leap.prm7 < com.leap.rst7 > com.leap.pdb

File:2018 12 coval1.png

2018-12-17T21:32:57Z

TBalius: TBalius uploaded a new version of "File:2018 12 coval1.png"

File:2018 12 coval1.png

2018-12-17T21:24:12Z

TBalius: TBalius uploaded a new version of "File:2018 12 coval1.png"

File:2018 12 coval1.png

2018-12-17T21:16:19Z

TBalius: TBalius uploaded a new version of "File:2018 12 coval1.png"

2018-12-17T17:58:47Z

TBalius:

2018-06-12T00:45:30Z

TBalius: /* run blastermaster.py */

This tutorial use the 3.7.2 beta version of dock release on April 17, 2015.

This is for a Linux environment and the scripts assume that you are running on SGE queuing system.

More information and tutorials, see [[DOCK_3.7]].

== set up directories and get databases ==

Create directory called "RotationProject"

create a python file called "0000.autodude_db_download.py"

# this gets the database from the autodude webpage

import sys, os
import urllib

system = 'abl1'
url = 'http://autodude.docking.org/dude_e_db2/'

print "url = " + url

#page=requests.get(url)

webfile = urllib.urlopen(url)
page = webfile.read()
webfile.close()

splitpage=page.split('\n')

for line in splitpage:
if system in line:
file = line.replace('"',' ').split()[2]
print url+file
urllib.urlretrieve(url+file,file)

# exit()

This python script will download the dockable db2 databases from the autodude webpage.

python /mnt/nfs/home/rstein/RotationProject/autodude_db_download.py

make a subdirectory called databases:

mkdir databases

go inside.

cd databases

make directories for ligands and decoys and move the corresponding files into those directories

mkdir decoys
mv decoys*db2.gz decoys

mkdir ligands
mv ligands*db2.gz ligands

download the ligand and decoy isomeric smiles file:

wget http://autodude.docking.org/abl1/decoys_final.ism
mv decoys_final.ism decoys.ism

note that the scripts expect the name to be decoys.ism, so we changed the name.

wget http://autodude.docking.org/abl1/actives_final.ism
mv actives_final.ism ligands.ism

== run be_blasti.py==

First we need to get our protein of interest from the protein databank (pdb). We will typiclly use a receptor with a ligand bound as is the case for pdbcode 2HYY, which is the Abl kinase domain in complex with imatinib (STI571, Glivec).

Note, in the following scripts, that DOCKBASE is a environment variable that point to the DOCK3.7 code. e.g.:

setenv DOCKBASE "/path2dock3.7/DOCK"

or

setenv DOCKBASE "/nfs/home/tbalius/zzz.github/DOCK"

Note, if you get an error ImportError: No module named Bio.PDB, then install the biopython as followed:

sudo yum install python-biopython

Create the following cshell script 0001.be_balsti_py.csh by using your favorite text editor (eg vim).

Note that the "#" symbol denotes a comment to explain what the script is doing.

#!/bin/csh

# this script calls be_blasti.py which creates a receptor and ligand file from a (list of) pdbcode(s).

# msms is a molecular surface generation program needed for be_blasti.py to run
# which is put in your path
set path = ( /nfs/home/tbalius/zzz.programs/msms $path )
# you will need to have msms on you system.

set list = "2HYY" # or use `cat filename` to list your pdb codes here from a text file like pdblist_rat, to loop over each variable (pdb code) later
#set list = `cat $1`
#set list = `cat /nfs/work/users/tbalius/VDR/Enrichment/pdblist_rat `

# CHANGE THIS, according to where the magic is going to happen
#set mountdir = "/mnt/nfs/work/users/tbalius/VDR/"
set mountdir = `pwd`

# loop over pdbnames e.g. 1DB1 or list
foreach pdbname ( $list )

echo " ${pdbname} "

# for each pdb makes a directory with its name
set workdir = ${mountdir}/${pdbname}

## so you don't blow away stuff; continue means STOP here and continue with next pdb from list
if ( -s $workdir ) then
echo "$workdir exits"
continue
endif

mkdir -p ${workdir}
cd ${workdir}

# the atom type definition is needed for msms which is sym-linked into the cwd
ln -s /nfs/home/tbalius/zzz.programs/msms/atmtypenumbers .
# carbs are disregarded as ligands! if it is: carbohydrate instead of noncarbohydrate
# renumber renumbers the residue number
python $DOCKBASE/proteins/pdb_breaker/be_blasti.py --pdbcode $pdbname nocarbohydrate original_numbers | tee -a pdbinfo_using_biopython.log

# error checking looks for receptor and ligand file which should be produced by be_blasti.py
if !(-s rec.pdb) then
echo "rec.pdb is not found"
endif

mv rec.pdb temp.pdb
grep -v TER temp.pdb | grep -v END > rec.pdb

rm temp.pdb

# be_blasti.py produces peptide which may be used as a ligand if no other ligand is produced
if (-s lig.pdb) then
sed -e "s/HETATM/ATOM /g" lig.pdb > xtal-lig.pdb
else if (-s pep.pdb) then ## if no ligand and peptide
sed -e "s/HETATM/ATOM /g" pep.pdb > xtal-lig.pdb
else
echo "Warning: No ligand or peptid."
endif

end # system

running 0001.be_balsti_py.csh will run a script that comes with dock called be_blasti.

issue the following command to run the script:
csh 0001.be_balsti_py.csh

It will do the following
# Download the pdb file from the web
# Break the file into rec and ligand components

Note that you will need to have msms on you system.

[[get msms]]

For Shoichet lab members msms is already installed.

check to make sure that the right ligand was selected and the the residue is not missing anything of importance.
If this automatic procedure has not prepared these files correctly, then modify them.

Visualize them with chimera or an alternive visualization program like pymol.

cd 2HYY

chimera rec.pdb lig.pdb

[[File:rec_lig_2HYY.png|thumb|center|375px|2HYY, the receptor and ligand generated from be_blasti.py.]]

== run blastermaster.py ==

Write (paste what follows) the following script using a text editor like vi. This script creates the files necessary for docking including the spheres (for orienting the ligands/decoys) and grids (for scoring the ligand/decoy poses)

WARNING: if you copy and pasted the script make sure there is no space before the "EOF"; this is because "EOF" designates the "end of file" for the cat command, if a space is there it wont stop cat'ing.

0002.blastermaster.csh

#!/bin/csh

# This script runs Ryan's blastermaster python masterscript for generating everything that dock needs, i.e. grids, spheres
# Run on sgehead as jobs are submitted to the queue

# list is same as in 001... script
set list = "2HYY"
#set list = `cat $1`
#set list = `cat /nfs/work/users/tbalius/VDR/Enrichment/pdblist_all `

set mountdir = `pwd`
#set mountdir = "/nfs/work/users/tbalius/VDR/"

# loop over all pdb(s)
foreach pdbname ( $list )

echo "${pdbname}"

set workdir = ${mountdir}/${pdbname}

# checks that 001 ran successfully and produced the directory structure as expected
# if not stops with current pdb code and continues with next one in list
if ! ( -s $workdir ) then
echo "$workdir does not exit"
continue
endif

cd $workdir

#cat xtal-lig_ori.pdb | awk '{if ($1 == "ATOM" || $1 == "HETATM"){print $0}}' | sed -e "s/HETATM/ATOM /g" > xtal-lig.pdb

# the following lines create a qsub script which submits blastermaster to the queue
cat <<EOF > qsub.csh
#!/bin/csh
#\$ -cwd
#\$ -j yes
#\$ -o stderr
#\$ -q all.q
cd $workdir
python $DOCKBASE/proteins/blastermaster/blastermaster.py --addhOptions=" -HIS -FLIPs " -v
EOF

qsub qsub.csh

end # pdbname
# going to the next pdb

# this will produce two directories:
# 1) working - contains all input and output files that are generated; not needed afterwards but as a reference
# 2) dockfiles - contains everything that is needed to run dock (copied from working)
# grids
# trim.electrostatics.phi
# vdw.vdw
# vdw.bmp
# ligand.desolv.heavy
# ligand.desolv.hydrogen
# spheres
# matching_spheres.sph

== Modifications to INDOCK ==

Before run the enrichment calculations consider modifying your INDOCK file.

For example, you might want to modify the the energy threshold maximum that controls poses (or molecules) written to a mol2 files (that is, for molecule that do not score below this thresold, poses will not written ).

Change from
mol2_score_maximum -10.0
To:
mol2_score_maximum +10.0

Another parameter to consider changing might be the heavy atom count maximum:

Change from
atom_maximum 25
to
atom_maximum 100

== Visualize the docking spheres ==

Use the showsphere program:

$DOCKBASE/proteins/showsphere/bin/showsphere

$DOCKBASE/proteins/showsphere/doshowsph.csh file.sph 1 file.pdb

for example:
$DOCKBASE/proteins/showsphere/doshowsph.csh 2HYY/dockfiles/matching_spheres.sph 1 matching_spheres.pdb

Alternatively, use the following simple cshell/awk script to convert spheres to pdb format:

cat << EOF > sphere2pdb.csh
#!/bin/csh -f
awk '$0!~/e/{ \
printf("ATOM %5d C SPH%5d%12.3f%8.3f%8.3f%6.2f%6.2f\nTER\n", \$1, \$1, \$2, \$3, \$4, 1, \$5)}' \$1
EOF

This command statement will write the commands script to a file called sphere2pdb.csh. Make sure there is not space before the second EOF (end of file) above. To run the command:

csh sphere2pdb.csh file.sph > file.pdb

cd 2HYY/working

chimera rec.pdb matching_spheres.pdb

[[File:Spheres.png|thumb|center|375px|The spheres generated from blastermaster.py in relation to the receptor.]]

The box used for scoring can also be visualized in chimera with the following command:

chimera rec.pdb matching_spheres.pdb box

[[File:Spheresbox.png|thumb|center|375px|The spheres and box generated from blastermaster.py in relation to the receptor.]]

== run enrichment calculations ==

Submit an enrichment calculation via 0003.lig-decoy_enrichment_submit.csh

We recommend using this method, as it uses the DOCK submission infrastructure.

* Write a file called 0003.lig-decoy_enrichment_submit.csh

#!/bin/csh

#This script provides a alternative way to dock a DUD-e like ligand-decoy-database for the enrichment evaluation of actives over decoys
#It assumes that ligands and decoys have been pre-prepation (see script blablabla_ToDo) which needs to be run in SF.

set filedir = "/mnt/nfs/home/jklyu/work/DOCK_tutorial" #CHANGE THIS
# this is where the work is done:
set mountdir = $filedir # Might CHANGE THIS
set dude_dir = "/mnt/nfs/home/jklyu/work/DOCK_tutorial/databases" # should contain decoy.smi and ligand.smi for ROC script 00005...csh
## TO DO - rename this outside in the dir structure and call in blbalbalbabla script
if (-s $dude_dir) then
echo " $dude_dir exist"
else
# this is something to modified in future.
# probably better to exit if it is not there.
echo "databases do not exist. "
echo "consider making a symbolic link to the database files"
endif

set list = "2HYY" # CHANGE THIS (pdbname)
foreach pdbname ( $list )
# creates "ligands" and "decoys" and has the aim to dock all of the subsets for those two
foreach db_type ( "ligands" "decoys" )
set workdir1 = "${mountdir}/${pdbname}/${db_type}"
set workdir2 = "${mountdir}/${pdbname}"
#
echo $mountdir
echo $workdir1
echo $workdir2
#
mkdir -p ${workdir1}
cd ${workdir1}
#creat dirlist for *.db2.gz files prepared for docking
ls ${dude_dir}/${db_type}/*.db2.gz > ${db_type}_files.txt
#copy the files needed for dock
cp ${workdir2}/INDOCK ${workdir1}
ln -s ${workdir2}/dockfiles/ ${workdir1}
#use dirlist to creat chunks for job submission
python /nfs/home/tbalius/zzz.github/DOCK/docking/setup/setup_db2_zinc15_file_number.py ./ chunk ./${db_type}_files.txt 500 count
#
csh $DOCKBASE/docking/submit/submit.csh

end # db_type
end # pdbname

* Run the above script

csh 0003.lig-decoy_enrichment_submit.csh

== combine scores and poses ==

Write this file as 0004.combineScoresAndPoses.csh. This script combines your docking runs (for ligands and decoys) and stores them in the extract_all.txt file. It also creates a .mol2 file containing all top scoring poses for each ligand/decoy.

Note that if you used the alternative (which is preferred) submission in the above step, then you will need to modify the directory structure in the below script.

From:
${mountdir}/${pdbname}/ligands-decoys/${db_type}/allChunksCombined
To:
${mountdir}/${pdbname}/${db_type}/allChunksCombined

#!/bin/csh

# This script combines the results from the ligand-decoy run 0003 (all chunks) into a combine file containing dock scores from OUTDOCK files
# Three files are produced (one for lig, decoy and both)
# and: a file which has top poses as specified (e.g. top 1000 molecules with 2 poses each); two files (for lig and for decoys)

# to remove dir
# rm -fr pdbs/3O1D/ligands-decoys/ligands/allChunksCombined/ pdbs/3O1D/ligands-decoys/decoys/allChunksCombined/ pdbs/3O1D/ligands-decoys/dockedLigDecoyCombined/

set filedir = "/mnt/nfs/home/rstein/RotationProject"
set mountdir = "/mnt/nfs/home/rstein/RotationProject"
set d37 = $DOCKBASE/analysis/

cd $mountdir

set list = "2HYY"
#set list = `cat filename`
#set list = `cat $1`

foreach pdbname ( $list )

foreach db_type ( "ligands" "decoys" )

set workdir = ${mountdir}/${pdbname}/ligands-decoys/${db_type}/allChunksCombined

echo $pdbname

#ls -l ${mountdir}/${pdbname}/${db_type}/

mkdir -p ${workdir}
cd ${workdir}

# creates a file called dirlist that contains the full path of all directories with docked runs (chunks)
ls -ld ${mountdir}/${pdbname}/ligands-decoys/${db_type}/* | awk '/chunk/{print $9}' > dirlist

#ls -ld ${mountdir}/${pdbname}/ligands-decoys/${db_type}/*

# for debuging
#echo "print $db_type dirlist:"
#cat dirlist

# script extracts scores from all docking runs specified in dirlist
$d37/extract_all.py
# script gets poses for top scoring molecules and produces poses.mol2 (default name)
$d37/getposes.py -d ${mountdir}/${pdbname}/ligands-decoys/${db_type}

end # db_type

## combine decoyes and actives
set workdir = ${mountdir}/${pdbname}/ligands-decoys/dockedLigDecoyCombined

rm -rf ${workdir}
mkdir -p ${workdir}
cd ${workdir}

cat ${mountdir}/${pdbname}/ligands-decoys/ligands/allChunksCombined/dirlist ${mountdir}/${pdbname}/ligands-decoys/decoys/allChunksCombined/dirlist > dirlist

# for debuging
#echo "print ALL dirlist"
#cat dirlist

$d37/extract_all.py
#$d37/getposes.py -d ${mountdir}/${pdbname} # doesn't work yet; not really needed
#getposes.py -z -l 1000 -x 2 -f extract_all.sort.uniq.txt -o ligands.1000.mol2 -d /mnt/nfs/work/users/fischer/VDR/27Jan2014_learningDOCKrgc/Enrichment/1DB1/DOCKING/ligands

end # pdbname

[[File:liganddock.png|thumb|center|375px|The pose generated from a ligand (purple) docked to Abl1 compared to imatinib (blue).]]

[[File:decoydock.png|thumb|center|375px|The pose generated from a decoy (green) docked to Abl1 compared to imatinib (blue).]]

== create AUC plot of ligands and decoys ==

Make sure that you are using a version of python that has matplotlib/numpy/scipy modules:

For Shoichet user source the following:
source /nfs/soft/python/envs/complete/latest/env.sh
or
source /nfs/soft/python/envs/complete/latest/env.csh
Also, if you are running the processing script remotely make sure to use X11 forwarding by including a -X in the ssh command.
eg:
ssh gimel.ucsf.bkslab.org -X

Write a file called 0005.AUCplot_of-lig-decoys.csh.

#!/bin/csh

# This script creates a log adjusted AUC (ROC) plot with ligand vs decoy results
# need X11 forwarding enabled when running remotely (ssh sgehead -X).

set filedir = "/mnt/nfs/home/rstein/RotationProject" #CHANGE THIS
set mountdir = "/mnt/nfs/home/rstein/RotationProject" #CHANGE THIS
set d37 = $DOCKBASE/analysis
set dude_dir = "/mnt/nfs/home/rstein/RotationProject/databases" # should contain decoy.smi and ligands.smi

# ln -s /mnt/nfs/work/users/fischer/VDR/lig-decoy-db/ligands.mod.smi /mnt/nfs/work/users/fischer/VDR/lig-decoy-db/ligands.smi

# CHANGE THIS
set list = "2HYY"
#set list = `cat filename`
#set list = `cat $1`

foreach pdbname ( $list )

set workdir = ${mountdir}/${pdbname}/ROC_ligdecoy/

# This script will not work without the following line:
echo "HERE is the HAWK"

# checks that previous script 0003 has produced mol2 files
if (! ( -s $mountdir/${pdbname}/ligands-decoys/decoys/allChunksCombined/poses.mol2) && ! (-s $mountdir/${pdbname}/ligands-decoys/ligands/allChunksCombined/poses.mol2 )) then
ls -l $mountdir/${pdbname}/ligands-decoys/decoys/allChunksCombined/poses.mol2
ls -l $mountdir/${pdbname}/ligands-decoys/ligands/allChunksCombined/poses.mol2
echo "skipping ${pdbname}. cannot generate ROC"
continue
endif

rm -rf $workdir
mkdir -p $workdir
cd $workdir

#wget http://dude.docking.org/targets/aa2ar/actives_final.ism

# reads ZINC ids (ligand or decoy molecule names)
# everything
awk '{print $2}' $dude_dir/decoys.ism > decoys.name # note that you may have to change the column ($2) based on where the SMILES codes are
awk '{printf "%9s\n", $3}' $dude_dir/ligands.ism > ligands.name # note that you may have to change the column ($3) based on where the SMILES codes are
#things that finished docking
awk '{print $3}' $mountdir/${pdbname}/ligands-decoys/decoys/allChunksCombined/extract_all.sort.uniq.txt > decoys.finished.name
awk '{print $3}' $mountdir/${pdbname}/ligands-decoys/ligands/allChunksCombined/extract_all.sort.uniq.txt > ligands.finished.name

cat ${mountdir}/${pdbname}/ligands-decoys/ligands/allChunksCombined/dirlist ${mountdir}/${pdbname}/ligands-decoys/decoys/allChunksCombined/dirlist > dirlist

#which enrich.py
set enrich_py = $d37/enrich.py
set plots_py = $d37/plots.py

pwd
# calculates AUCs, stores in txt file which is then plotted for finished ligands and decoys
python ${enrich_py} -i . -o . --ligand-file=ligands.finished.name --decoy-file=decoys.finished.name
python ${plots_py} -i . -o . --ligand-file=ligands.finished.name --decoy-file=decoys.finished.name -l $pdbname

mv roc.txt roc.finished.txt
mv roc_own.txt roc_own.finished.txt
mv roc_own.png roc_own.finished.png

#
# calculates AUCs, stores in txt file which is then plotted for all ligands and decoys
# - i is the flag for the input directory, this dir should contain the extract_all.sort.uniq.txt.
# the scripts enrich_py and plots_py will go through the extract file and look for the ligand and decoy names.
# when it finds them it will populate the ROC cruve. these values are devied by the total number of ligand or decoys.
# note that often not all ligands and not all decoys finish so the point (1,1) is always included and interpolations is performed . . .
#
#python ${enrich_py} -i $mountdir/${pdbname}/ligands-decoys/dockedLigDecoyCombined/ -o . --ligand-file=ligands.name --decoy-file=decoys.name
#python ${plots_py} -i $mountdir/${pdbname}/ligands-decoys/dockedLigDecoyCombined/ -o . --ligand-file=ligands.name --decoy-file=decoys.name -l $pdbname
python ${enrich_py} -i . -o . --ligand-file=ligands.name --decoy-file=decoys.name
python ${plots_py} -i . -o . --ligand-file=ligands.name --decoy-file=decoys.name -l $pdbname

end #pbdname

The ROC plot for ligands and decoys docking to Abl1 looks like this:

[[File:AblROCplot.png|thumb|center|375px|The ROC plot generated for the ligands and decoys.]]

This shows poor enrichment, though it is better than random (indicated by the dotted line). This is to be expected as docking to kinases is notoriously difficult. Two changes to our procedure would be to delete some of the spheres to focus docking on a particular region of the protein we are most interested in, as well as increasing the polarity of the hinge region of Abl1 to allow for more hydrogen bonding.

== improving docking ==
===Modify matching spheres===
The docking could be potentially improved by modifying the spheres and orienting the ligands/decoys to only those spheres within the hinge region of Abl1. The following image show modified spheres:

[[File:modspheres.png|thumb|center|375px|Modified spheres to improve docking.]]

To modify spheres: first, convert them to pdb format; then, visualize them in your favorite program (pymol, chimera, etc.); then, delete, move, or add atoms to this file; and finally convert it back to the sph format.

As discussed above [[http://wiki.bkslab.org/index.php/DOCK_3.7_2015/04/15_abl1_Tutorial#Visualize_the_docking_spheres]], you may use doshowsph.csh to convert the spheres to pdb format.

Here is the program that will convert a pdb file into a sphere file (this program take 2 inputs: name of pdbfile to convert and the name of the sphere file that you want to create.).
$DOCKBASE/proteins/pdbtosph/bin/pdbtosph

for example:
$DOCKBASE/proteins/pdbtosph/bin/pdbtosph matching_spheres_mod.pdb matching_spheres_mod.sph

The 0003.lig-decoy_enrichment.csh would be modified to include this line and renamed to 0003.lig-decoy_enrichment_mod_sph.csh:

!/bin/csh

#This script docks a DUD-e like ligand-decoy-database to evaluate the enrichment performance of actives over decoys
#It assumes that ligands and decoys have been pre-prepation (see script blablabla_ToDo) which needs to be run in SF.

# filedir is where your rec.pdb and xtal-lig.pdb and dockfiles directory live
set filedir = "/mnt/nfs/home/rstein/RotationProject" #CHANGE THIS
# this is where the work is done:
set mountdir = $filedir # Might CHANGE THIS
set dude_dir = "/mnt/nfs/home/rstein/RotationProject/databases" # should contain decoy.smi and ligand.smi for ROC script 00005...csh
## TO DO - rename this outside in the dir structure and call in blbalbalbabla script
if (-s $dude_dir) then
echo " $dude_dir exist"
else
# this is something to modified in future.
# probably better to exit if it is not there.
echo "databases do not exist. "
echo "consider making a symbolic link to the database files"
#echo "making a symbolic link:"
#echo "ln -s /mnt/nfs/work/users/fischer/VDR/27Jan2014_learningDOCKrgc/databases_all_xtal-ligand_decoy $dude_dir"
#ln -s /mnt/nfs/work/users/fischer/VDR/27Jan2014_learningDOCKrgc/databases_all_xtal-ligand_decoy $dude_dir
endif

# change if you want to use a different or consistent dock version
set dock = ${DOCKBASE}/docking/DOCK/bin/dock64

set list = "2HYY"
#set list = `cat $1`
#set list = `cat file`
# CHANGE THIS (pdbname)
foreach pdbname ( $list )

# creates "ligands" and "decoys" and has the aim to dock all of the subsets for those two
foreach db_type ( "ligands" "decoys" )
set workdir1 = "${mountdir}/${pdbname}/ligands-decoys_sphmod/${db_type}"

mkdir -p ${workdir1}
cd ${workdir1}
# puts dockfiles in the right relative-path that INDOCK file expects
ln -s $filedir/${pdbname}/dockfiles .

set count = '1'

# loop over database files to put each into a seperate chunk
foreach dbfile (`ls $dude_dir/${db_type}/${db_type}*.db2.gz`)

echo $dbfile

set chunk = "chunk$count"

set workdir2 = ${workdir1}/$chunk
## so you don't blow away stuff
if ( -s $workdir2 ) then
echo "$workdir2 exits"
continue
endif

#rm -rf ${workdir}
mkdir -p ${workdir2}
cd ${workdir2}

# copy INDOCK file of choice in right location
#cp $filedir/zzz.dock3_input/INDOCK .
#cp $filedir/INDOCK_match20K INDOCK
#cp $filedir/INDOCK_5k_TolerantClash INDOCK # CHANGE THIS
cp $filedir/${pdbname}/INDOCK .
# modified the dock file using sed. here we change some key sampling parameters; sed -i changes input file internally (overwrites), -e changes file externally (pipes it to screen or into file if redirected)
#sed -i "s/bump_maximum 50.0/bump_maximum 500.0/g" INDOCK
#sed -i "s/bump_rigid 50.0/bump_rigid 500.0/g" INDOCK
#sed -i "s/check_clashes yes/check_clashes no/g" INDOCK
sed -i "s/receptor_sphere_file ..\/dockfiles\/matching_spheres.sph/receptor_sphere_file ..\/..\/..\/working\/matching_spheres_mod.sph/g" INDOCK

ln -s $dbfile .

set dbf = `ls *.gz`

echo "./$dbf"

# says what to dock and where it sits
echo "./$dbf" > split_database_index\
# writes submission script that runs dock on the sgehead queue
cat <<EOF > DOCKING_${db_type}.csh
#\$ -S /bin/csh
#\$ -cwd
#\$ -q all.q
#\$ -o stdout
#\$ -e stderr

cd ${workdir2}
echo "starting . . ."
date
echo $dock
$dock
date
echo "finished . . ."

EOF
qsub DOCKING_${db_type}.csh
# alternatively if you don't want to run it on the queue but locally comment in this instead:
#csh DOCKING_${lig_type}.csh &

@ count = ${count} + 1
# counter is chuch dir

end # dbfile
end # db_type
end # pdbname

=== make the hing region more polar ===
see the following page:
[[DOCK_3.7_tart]]

==Virtual Screening==

===database setup===

This part of the tutorial is tailored for shoichet lab use. An outside user of dock might need to deviate from what is described.

Go to zinc and select your compounds of interested:
'[http://zinc15.docking.org/tranches/home http://zinc15.docking.org/tranches/home]'

This is the tranches page which allows users to select the region of chemical space of interest.

lets select the fragment preset.

on lets download the index file. This file contains the location of each database on our cluster. outside users will need to download the databases themselves.

now lets setup the directorys for docking by running the following script:

python /nfs/home/tbalius/zzz.github/DOCK/docking/setup/setup_db2_zinc15_file_number.py ./ vs_frag /nfs/work/tbalius/database_ph4/frags.txt 500 count

The file /nfs/work/tbalius/database_ph4/frags.txt should be changed to that you downloaded from ZINC.

The above script has 5 parameters:
:(1) path where directories will be located (present directory);
:(2) prefix name of the directories;
:(3) the file that contains the db2 files locations;
:(4) the number of directories to be created; and
:(5) the type of run: count (evenly distributes the db2 file among the dirs, this is much faster than the other options), size (It will try and make the directory of equal size), or both (will try and satisfy both criteria).

Note that this script is avable in later beta versons of DOCK3.7.

===submitting the docking calculations===

This script will submit a job to the queue for each of the docking directorys created by the setup script.

$DOCKBASE/docking/submit/submit.csh

DOCK3.7 is a serial program and is parallelized by submiting many serial jobs to the queue.

===combining the results===

After your docking jobs have all completed, This script will combine all your results into an extract_all file.
$DOCKBASE/analysis/extract_all.py

This script will create a mol2 file with the top scoring molecules:
$DOCKBASE/analysis/getposes.py

=== curating and hit-picking ===

Typically we will visualize the ligands in UCSF Chimera using the veiwdock tool (using the DOCK4, 5, or 6 format).

DOCK 3.7 2018/06/05 abl1 Tutorial

2018-06-12T00:44:36Z

TBalius: /* run blastermaster.py */

This tutorial use the 3.7.2 beta version of dock release on April 17, 2015.

This is for a Linux environment and the scripts assume that you are running on SGE queuing system.

More information and tutorials, see [[DOCK_3.7]].

== set up directories and get databases ==

Create directory called "RotationProject"

create a python file called "0000.autodude_db_download.py"

# this gets the database from the autodude webpage

import sys, os
import urllib

system = 'abl1'
url = 'http://autodude.docking.org/dude_e_db2/'

print "url = " + url

#page=requests.get(url)

webfile = urllib.urlopen(url)
page = webfile.read()
webfile.close()

splitpage=page.split('\n')

for line in splitpage:
if system in line:
file = line.replace('"',' ').split()[2]
print url+file
urllib.urlretrieve(url+file,file)

# exit()

This python script will download the dockable db2 databases from the autodude webpage.

python /mnt/nfs/home/rstein/RotationProject/autodude_db_download.py

make a subdirectory called databases:

mkdir databases

go inside.

cd databases

make directories for ligands and decoys and move the corresponding files into those directories

mkdir decoys
mv decoys*db2.gz decoys

mkdir ligands
mv ligands*db2.gz ligands

download the ligand and decoy isomeric smiles file:

wget http://autodude.docking.org/abl1/decoys_final.ism
mv decoys_final.ism decoys.ism

note that the scripts expect the name to be decoys.ism, so we changed the name.

wget http://autodude.docking.org/abl1/actives_final.ism
mv actives_final.ism ligands.ism

== run be_blasti.py==

First we need to get our protein of interest from the protein databank (pdb). We will typiclly use a receptor with a ligand bound as is the case for pdbcode 2HYY, which is the Abl kinase domain in complex with imatinib (STI571, Glivec).

Note, in the following scripts, that DOCKBASE is a environment variable that point to the DOCK3.7 code. e.g.:

setenv DOCKBASE "/path2dock3.7/DOCK"

or

setenv DOCKBASE "/nfs/home/tbalius/zzz.github/DOCK"

Note, if you get an error ImportError: No module named Bio.PDB, then install the biopython as followed:

sudo yum install python-biopython

Create the following cshell script 0001.be_balsti_py.csh by using your favorite text editor (eg vim).

Note that the "#" symbol denotes a comment to explain what the script is doing.

#!/bin/csh

# this script calls be_blasti.py which creates a receptor and ligand file from a (list of) pdbcode(s).

# msms is a molecular surface generation program needed for be_blasti.py to run
# which is put in your path
set path = ( /nfs/home/tbalius/zzz.programs/msms $path )
# you will need to have msms on you system.

set list = "2HYY" # or use `cat filename` to list your pdb codes here from a text file like pdblist_rat, to loop over each variable (pdb code) later
#set list = `cat $1`
#set list = `cat /nfs/work/users/tbalius/VDR/Enrichment/pdblist_rat `

# CHANGE THIS, according to where the magic is going to happen
#set mountdir = "/mnt/nfs/work/users/tbalius/VDR/"
set mountdir = `pwd`

# loop over pdbnames e.g. 1DB1 or list
foreach pdbname ( $list )

echo " ${pdbname} "

# for each pdb makes a directory with its name
set workdir = ${mountdir}/${pdbname}

## so you don't blow away stuff; continue means STOP here and continue with next pdb from list
if ( -s $workdir ) then
echo "$workdir exits"
continue
endif

mkdir -p ${workdir}
cd ${workdir}

# the atom type definition is needed for msms which is sym-linked into the cwd
ln -s /nfs/home/tbalius/zzz.programs/msms/atmtypenumbers .
# carbs are disregarded as ligands! if it is: carbohydrate instead of noncarbohydrate
# renumber renumbers the residue number
python $DOCKBASE/proteins/pdb_breaker/be_blasti.py --pdbcode $pdbname nocarbohydrate original_numbers | tee -a pdbinfo_using_biopython.log

# error checking looks for receptor and ligand file which should be produced by be_blasti.py
if !(-s rec.pdb) then
echo "rec.pdb is not found"
endif

mv rec.pdb temp.pdb
grep -v TER temp.pdb | grep -v END > rec.pdb

rm temp.pdb

# be_blasti.py produces peptide which may be used as a ligand if no other ligand is produced
if (-s lig.pdb) then
sed -e "s/HETATM/ATOM /g" lig.pdb > xtal-lig.pdb
else if (-s pep.pdb) then ## if no ligand and peptide
sed -e "s/HETATM/ATOM /g" pep.pdb > xtal-lig.pdb
else
echo "Warning: No ligand or peptid."
endif

end # system

running 0001.be_balsti_py.csh will run a script that comes with dock called be_blasti.

issue the following command to run the script:
csh 0001.be_balsti_py.csh

It will do the following
# Download the pdb file from the web
# Break the file into rec and ligand components

Note that you will need to have msms on you system.

[[get msms]]

For Shoichet lab members msms is already installed.

check to make sure that the right ligand was selected and the the residue is not missing anything of importance.
If this automatic procedure has not prepared these files correctly, then modify them.

Visualize them with chimera or an alternive visualization program like pymol.

cd 2HYY

chimera rec.pdb lig.pdb

[[File:rec_lig_2HYY.png|thumb|center|375px|2HYY, the receptor and ligand generated from be_blasti.py.]]

== run blastermaster.py ==

Write (paste what follows) the following script using a text editor like vi. This script creates the files necessary for docking including the spheres (for orienting the ligands/decoys) and grids (for scoring the ligand/decoy poses)

WARNING: if you copy and pasted the script make sure there is no space before the "EOF"; this is because "EOF" designates the "end of file" for the cat command, if a space is there it wont stop cat'ing.

0002.blastermaster.csh

#!/bin/csh

# This script runs Ryan's blastermaster python masterscript for generating everything that dock needs, i.e. grids, spheres
# Run on sgehead as jobs are submitted to the queue

# list is same as in 001... script
set list = "2HYY"
#set list = `cat $1`
#set list = `cat /nfs/work/users/tbalius/VDR/Enrichment/pdblist_all `

set mountdir = `pwd`
#set mountdir = "/nfs/work/users/tbalius/VDR/"

# loop over all pdb(s)
foreach pdbname ( $list )

echo "${pdbname}"

set workdir = ${mountdir}/${pdbname}

# checks that 001 ran successfully and produced the directory structure as expected
# if not stops with current pdb code and continues with next one in list
if ! ( -s $workdir ) then
echo "$workdir does not exit"
continue
endif

cd $workdir

#cat xtal-lig_ori.pdb | awk '{if ($1 == "ATOM" || $1 == "HETATM"){print $0}}' | sed -e "s/HETATM/ATOM /g" > xtal-lig.pdb

# the following lines create a qsub script which submits blastermaster to the queue
cat <<EOF > qsub.csh
#!/bin/csh
#\$ -cwd
#\$ -j yes
#\$ -o stderr
#\$ -q all.q
cd $workdir
python $DOCKBASE/proteins/blastermaster/blastermaster.py --addhOptions=" -HIS -FLIPs " -v
EOF

qsub qsub.csh

end # pdbname
# going to the next pdb

# this will produce two directories:
# 1) working - contains all input and output files that are generated; not needed afterwards but as a reference
# 2) dockfiles - contains everything that is needed to run dock (copied from working)
# grids
# trim.electrostatics.phi
# vdw.vdw
# vdw.bmp
# ligand.desolv.heavy
# ligand.desolv.hydrogen
# spheres
# matching_spheres.sph

Before run the enrichment calculations consider modifying your INDOCK file.

For example, you might want to modify the the energy threshold maximum that controls poses (or molecules) written to a mol2 files (that is, for molecule that do not score below this thresold, poses will not written ).

Change from
mol2_score_maximum -10.0
To:
mol2_score_maximum +10.0

Another parameter to consider changing might be the heavy atom count maximum:

Change from
atom_maximum 25
to
atom_maximum 100

== Visualize the docking spheres ==

Use the showsphere program:

$DOCKBASE/proteins/showsphere/bin/showsphere

$DOCKBASE/proteins/showsphere/doshowsph.csh file.sph 1 file.pdb

for example:
$DOCKBASE/proteins/showsphere/doshowsph.csh 2HYY/dockfiles/matching_spheres.sph 1 matching_spheres.pdb

Alternatively, use the following simple cshell/awk script to convert spheres to pdb format:

cat << EOF > sphere2pdb.csh
#!/bin/csh -f
awk '$0!~/e/{ \
printf("ATOM %5d C SPH%5d%12.3f%8.3f%8.3f%6.2f%6.2f\nTER\n", \$1, \$1, \$2, \$3, \$4, 1, \$5)}' \$1
EOF

This command statement will write the commands script to a file called sphere2pdb.csh. Make sure there is not space before the second EOF (end of file) above. To run the command:

csh sphere2pdb.csh file.sph > file.pdb

cd 2HYY/working

chimera rec.pdb matching_spheres.pdb

[[File:Spheres.png|thumb|center|375px|The spheres generated from blastermaster.py in relation to the receptor.]]

The box used for scoring can also be visualized in chimera with the following command:

chimera rec.pdb matching_spheres.pdb box

[[File:Spheresbox.png|thumb|center|375px|The spheres and box generated from blastermaster.py in relation to the receptor.]]

== run enrichment calculations ==

Submit an enrichment calculation via 0003.lig-decoy_enrichment_submit.csh

We recommend using this method, as it uses the DOCK submission infrastructure.

* Write a file called 0003.lig-decoy_enrichment_submit.csh

#!/bin/csh

#This script provides a alternative way to dock a DUD-e like ligand-decoy-database for the enrichment evaluation of actives over decoys
#It assumes that ligands and decoys have been pre-prepation (see script blablabla_ToDo) which needs to be run in SF.

set filedir = "/mnt/nfs/home/jklyu/work/DOCK_tutorial" #CHANGE THIS
# this is where the work is done:
set mountdir = $filedir # Might CHANGE THIS
set dude_dir = "/mnt/nfs/home/jklyu/work/DOCK_tutorial/databases" # should contain decoy.smi and ligand.smi for ROC script 00005...csh
## TO DO - rename this outside in the dir structure and call in blbalbalbabla script
if (-s $dude_dir) then
echo " $dude_dir exist"
else
# this is something to modified in future.
# probably better to exit if it is not there.
echo "databases do not exist. "
echo "consider making a symbolic link to the database files"
endif

set list = "2HYY" # CHANGE THIS (pdbname)
foreach pdbname ( $list )
# creates "ligands" and "decoys" and has the aim to dock all of the subsets for those two
foreach db_type ( "ligands" "decoys" )
set workdir1 = "${mountdir}/${pdbname}/${db_type}"
set workdir2 = "${mountdir}/${pdbname}"
#
echo $mountdir
echo $workdir1
echo $workdir2
#
mkdir -p ${workdir1}
cd ${workdir1}
#creat dirlist for *.db2.gz files prepared for docking
ls ${dude_dir}/${db_type}/*.db2.gz > ${db_type}_files.txt
#copy the files needed for dock
cp ${workdir2}/INDOCK ${workdir1}
ln -s ${workdir2}/dockfiles/ ${workdir1}
#use dirlist to creat chunks for job submission
python /nfs/home/tbalius/zzz.github/DOCK/docking/setup/setup_db2_zinc15_file_number.py ./ chunk ./${db_type}_files.txt 500 count
#
csh $DOCKBASE/docking/submit/submit.csh

end # db_type
end # pdbname

* Run the above script

csh 0003.lig-decoy_enrichment_submit.csh

== combine scores and poses ==

Write this file as 0004.combineScoresAndPoses.csh. This script combines your docking runs (for ligands and decoys) and stores them in the extract_all.txt file. It also creates a .mol2 file containing all top scoring poses for each ligand/decoy.

Note that if you used the alternative (which is preferred) submission in the above step, then you will need to modify the directory structure in the below script.

From:
${mountdir}/${pdbname}/ligands-decoys/${db_type}/allChunksCombined
To:
${mountdir}/${pdbname}/${db_type}/allChunksCombined

#!/bin/csh

# This script combines the results from the ligand-decoy run 0003 (all chunks) into a combine file containing dock scores from OUTDOCK files
# Three files are produced (one for lig, decoy and both)
# and: a file which has top poses as specified (e.g. top 1000 molecules with 2 poses each); two files (for lig and for decoys)

# to remove dir
# rm -fr pdbs/3O1D/ligands-decoys/ligands/allChunksCombined/ pdbs/3O1D/ligands-decoys/decoys/allChunksCombined/ pdbs/3O1D/ligands-decoys/dockedLigDecoyCombined/

set filedir = "/mnt/nfs/home/rstein/RotationProject"
set mountdir = "/mnt/nfs/home/rstein/RotationProject"
set d37 = $DOCKBASE/analysis/

cd $mountdir

set list = "2HYY"
#set list = `cat filename`
#set list = `cat $1`

foreach pdbname ( $list )

foreach db_type ( "ligands" "decoys" )

set workdir = ${mountdir}/${pdbname}/ligands-decoys/${db_type}/allChunksCombined

echo $pdbname

#ls -l ${mountdir}/${pdbname}/${db_type}/

mkdir -p ${workdir}
cd ${workdir}

# creates a file called dirlist that contains the full path of all directories with docked runs (chunks)
ls -ld ${mountdir}/${pdbname}/ligands-decoys/${db_type}/* | awk '/chunk/{print $9}' > dirlist

#ls -ld ${mountdir}/${pdbname}/ligands-decoys/${db_type}/*

# for debuging
#echo "print $db_type dirlist:"
#cat dirlist

# script extracts scores from all docking runs specified in dirlist
$d37/extract_all.py
# script gets poses for top scoring molecules and produces poses.mol2 (default name)
$d37/getposes.py -d ${mountdir}/${pdbname}/ligands-decoys/${db_type}

end # db_type

## combine decoyes and actives
set workdir = ${mountdir}/${pdbname}/ligands-decoys/dockedLigDecoyCombined

rm -rf ${workdir}
mkdir -p ${workdir}
cd ${workdir}

cat ${mountdir}/${pdbname}/ligands-decoys/ligands/allChunksCombined/dirlist ${mountdir}/${pdbname}/ligands-decoys/decoys/allChunksCombined/dirlist > dirlist

# for debuging
#echo "print ALL dirlist"
#cat dirlist

$d37/extract_all.py
#$d37/getposes.py -d ${mountdir}/${pdbname} # doesn't work yet; not really needed
#getposes.py -z -l 1000 -x 2 -f extract_all.sort.uniq.txt -o ligands.1000.mol2 -d /mnt/nfs/work/users/fischer/VDR/27Jan2014_learningDOCKrgc/Enrichment/1DB1/DOCKING/ligands

end # pdbname

[[File:liganddock.png|thumb|center|375px|The pose generated from a ligand (purple) docked to Abl1 compared to imatinib (blue).]]

[[File:decoydock.png|thumb|center|375px|The pose generated from a decoy (green) docked to Abl1 compared to imatinib (blue).]]

== create AUC plot of ligands and decoys ==

Make sure that you are using a version of python that has matplotlib/numpy/scipy modules:

For Shoichet user source the following:
source /nfs/soft/python/envs/complete/latest/env.sh
or
source /nfs/soft/python/envs/complete/latest/env.csh
Also, if you are running the processing script remotely make sure to use X11 forwarding by including a -X in the ssh command.
eg:
ssh gimel.ucsf.bkslab.org -X

Write a file called 0005.AUCplot_of-lig-decoys.csh.

#!/bin/csh

# This script creates a log adjusted AUC (ROC) plot with ligand vs decoy results
# need X11 forwarding enabled when running remotely (ssh sgehead -X).

set filedir = "/mnt/nfs/home/rstein/RotationProject" #CHANGE THIS
set mountdir = "/mnt/nfs/home/rstein/RotationProject" #CHANGE THIS
set d37 = $DOCKBASE/analysis
set dude_dir = "/mnt/nfs/home/rstein/RotationProject/databases" # should contain decoy.smi and ligands.smi

# ln -s /mnt/nfs/work/users/fischer/VDR/lig-decoy-db/ligands.mod.smi /mnt/nfs/work/users/fischer/VDR/lig-decoy-db/ligands.smi

# CHANGE THIS
set list = "2HYY"
#set list = `cat filename`
#set list = `cat $1`

foreach pdbname ( $list )

set workdir = ${mountdir}/${pdbname}/ROC_ligdecoy/

# This script will not work without the following line:
echo "HERE is the HAWK"

# checks that previous script 0003 has produced mol2 files
if (! ( -s $mountdir/${pdbname}/ligands-decoys/decoys/allChunksCombined/poses.mol2) && ! (-s $mountdir/${pdbname}/ligands-decoys/ligands/allChunksCombined/poses.mol2 )) then
ls -l $mountdir/${pdbname}/ligands-decoys/decoys/allChunksCombined/poses.mol2
ls -l $mountdir/${pdbname}/ligands-decoys/ligands/allChunksCombined/poses.mol2
echo "skipping ${pdbname}. cannot generate ROC"
continue
endif

rm -rf $workdir
mkdir -p $workdir
cd $workdir

#wget http://dude.docking.org/targets/aa2ar/actives_final.ism

# reads ZINC ids (ligand or decoy molecule names)
# everything
awk '{print $2}' $dude_dir/decoys.ism > decoys.name # note that you may have to change the column ($2) based on where the SMILES codes are
awk '{printf "%9s\n", $3}' $dude_dir/ligands.ism > ligands.name # note that you may have to change the column ($3) based on where the SMILES codes are
#things that finished docking
awk '{print $3}' $mountdir/${pdbname}/ligands-decoys/decoys/allChunksCombined/extract_all.sort.uniq.txt > decoys.finished.name
awk '{print $3}' $mountdir/${pdbname}/ligands-decoys/ligands/allChunksCombined/extract_all.sort.uniq.txt > ligands.finished.name

cat ${mountdir}/${pdbname}/ligands-decoys/ligands/allChunksCombined/dirlist ${mountdir}/${pdbname}/ligands-decoys/decoys/allChunksCombined/dirlist > dirlist

#which enrich.py
set enrich_py = $d37/enrich.py
set plots_py = $d37/plots.py

pwd
# calculates AUCs, stores in txt file which is then plotted for finished ligands and decoys
python ${enrich_py} -i . -o . --ligand-file=ligands.finished.name --decoy-file=decoys.finished.name
python ${plots_py} -i . -o . --ligand-file=ligands.finished.name --decoy-file=decoys.finished.name -l $pdbname

mv roc.txt roc.finished.txt
mv roc_own.txt roc_own.finished.txt
mv roc_own.png roc_own.finished.png

#
# calculates AUCs, stores in txt file which is then plotted for all ligands and decoys
# - i is the flag for the input directory, this dir should contain the extract_all.sort.uniq.txt.
# the scripts enrich_py and plots_py will go through the extract file and look for the ligand and decoy names.
# when it finds them it will populate the ROC cruve. these values are devied by the total number of ligand or decoys.
# note that often not all ligands and not all decoys finish so the point (1,1) is always included and interpolations is performed . . .
#
#python ${enrich_py} -i $mountdir/${pdbname}/ligands-decoys/dockedLigDecoyCombined/ -o . --ligand-file=ligands.name --decoy-file=decoys.name
#python ${plots_py} -i $mountdir/${pdbname}/ligands-decoys/dockedLigDecoyCombined/ -o . --ligand-file=ligands.name --decoy-file=decoys.name -l $pdbname
python ${enrich_py} -i . -o . --ligand-file=ligands.name --decoy-file=decoys.name
python ${plots_py} -i . -o . --ligand-file=ligands.name --decoy-file=decoys.name -l $pdbname

end #pbdname

The ROC plot for ligands and decoys docking to Abl1 looks like this:

[[File:AblROCplot.png|thumb|center|375px|The ROC plot generated for the ligands and decoys.]]

This shows poor enrichment, though it is better than random (indicated by the dotted line). This is to be expected as docking to kinases is notoriously difficult. Two changes to our procedure would be to delete some of the spheres to focus docking on a particular region of the protein we are most interested in, as well as increasing the polarity of the hinge region of Abl1 to allow for more hydrogen bonding.

== improving docking ==
===Modify matching spheres===
The docking could be potentially improved by modifying the spheres and orienting the ligands/decoys to only those spheres within the hinge region of Abl1. The following image show modified spheres:

[[File:modspheres.png|thumb|center|375px|Modified spheres to improve docking.]]

To modify spheres: first, convert them to pdb format; then, visualize them in your favorite program (pymol, chimera, etc.); then, delete, move, or add atoms to this file; and finally convert it back to the sph format.

As discussed above [[http://wiki.bkslab.org/index.php/DOCK_3.7_2015/04/15_abl1_Tutorial#Visualize_the_docking_spheres]], you may use doshowsph.csh to convert the spheres to pdb format.

Here is the program that will convert a pdb file into a sphere file (this program take 2 inputs: name of pdbfile to convert and the name of the sphere file that you want to create.).
$DOCKBASE/proteins/pdbtosph/bin/pdbtosph

for example:
$DOCKBASE/proteins/pdbtosph/bin/pdbtosph matching_spheres_mod.pdb matching_spheres_mod.sph

The 0003.lig-decoy_enrichment.csh would be modified to include this line and renamed to 0003.lig-decoy_enrichment_mod_sph.csh:

!/bin/csh

#This script docks a DUD-e like ligand-decoy-database to evaluate the enrichment performance of actives over decoys
#It assumes that ligands and decoys have been pre-prepation (see script blablabla_ToDo) which needs to be run in SF.

# filedir is where your rec.pdb and xtal-lig.pdb and dockfiles directory live
set filedir = "/mnt/nfs/home/rstein/RotationProject" #CHANGE THIS
# this is where the work is done:
set mountdir = $filedir # Might CHANGE THIS
set dude_dir = "/mnt/nfs/home/rstein/RotationProject/databases" # should contain decoy.smi and ligand.smi for ROC script 00005...csh
## TO DO - rename this outside in the dir structure and call in blbalbalbabla script
if (-s $dude_dir) then
echo " $dude_dir exist"
else
# this is something to modified in future.
# probably better to exit if it is not there.
echo "databases do not exist. "
echo "consider making a symbolic link to the database files"
#echo "making a symbolic link:"
#echo "ln -s /mnt/nfs/work/users/fischer/VDR/27Jan2014_learningDOCKrgc/databases_all_xtal-ligand_decoy $dude_dir"
#ln -s /mnt/nfs/work/users/fischer/VDR/27Jan2014_learningDOCKrgc/databases_all_xtal-ligand_decoy $dude_dir
endif

# change if you want to use a different or consistent dock version
set dock = ${DOCKBASE}/docking/DOCK/bin/dock64

set list = "2HYY"
#set list = `cat $1`
#set list = `cat file`
# CHANGE THIS (pdbname)
foreach pdbname ( $list )

# creates "ligands" and "decoys" and has the aim to dock all of the subsets for those two
foreach db_type ( "ligands" "decoys" )
set workdir1 = "${mountdir}/${pdbname}/ligands-decoys_sphmod/${db_type}"

mkdir -p ${workdir1}
cd ${workdir1}
# puts dockfiles in the right relative-path that INDOCK file expects
ln -s $filedir/${pdbname}/dockfiles .

set count = '1'

# loop over database files to put each into a seperate chunk
foreach dbfile (`ls $dude_dir/${db_type}/${db_type}*.db2.gz`)

echo $dbfile

set chunk = "chunk$count"

set workdir2 = ${workdir1}/$chunk
## so you don't blow away stuff
if ( -s $workdir2 ) then
echo "$workdir2 exits"
continue
endif

#rm -rf ${workdir}
mkdir -p ${workdir2}
cd ${workdir2}

# copy INDOCK file of choice in right location
#cp $filedir/zzz.dock3_input/INDOCK .
#cp $filedir/INDOCK_match20K INDOCK
#cp $filedir/INDOCK_5k_TolerantClash INDOCK # CHANGE THIS
cp $filedir/${pdbname}/INDOCK .
# modified the dock file using sed. here we change some key sampling parameters; sed -i changes input file internally (overwrites), -e changes file externally (pipes it to screen or into file if redirected)
#sed -i "s/bump_maximum 50.0/bump_maximum 500.0/g" INDOCK
#sed -i "s/bump_rigid 50.0/bump_rigid 500.0/g" INDOCK
#sed -i "s/check_clashes yes/check_clashes no/g" INDOCK
sed -i "s/receptor_sphere_file ..\/dockfiles\/matching_spheres.sph/receptor_sphere_file ..\/..\/..\/working\/matching_spheres_mod.sph/g" INDOCK

ln -s $dbfile .

set dbf = `ls *.gz`

echo "./$dbf"

# says what to dock and where it sits
echo "./$dbf" > split_database_index\
# writes submission script that runs dock on the sgehead queue
cat <<EOF > DOCKING_${db_type}.csh
#\$ -S /bin/csh
#\$ -cwd
#\$ -q all.q
#\$ -o stdout
#\$ -e stderr

cd ${workdir2}
echo "starting . . ."
date
echo $dock
$dock
date
echo "finished . . ."

EOF
qsub DOCKING_${db_type}.csh
# alternatively if you don't want to run it on the queue but locally comment in this instead:
#csh DOCKING_${lig_type}.csh &

@ count = ${count} + 1
# counter is chuch dir

end # dbfile
end # db_type
end # pdbname

=== make the hing region more polar ===
see the following page:
[[DOCK_3.7_tart]]

==Virtual Screening==

===database setup===

This part of the tutorial is tailored for shoichet lab use. An outside user of dock might need to deviate from what is described.

Go to zinc and select your compounds of interested:
'[http://zinc15.docking.org/tranches/home http://zinc15.docking.org/tranches/home]'

This is the tranches page which allows users to select the region of chemical space of interest.

lets select the fragment preset.

on lets download the index file. This file contains the location of each database on our cluster. outside users will need to download the databases themselves.

now lets setup the directorys for docking by running the following script:

python /nfs/home/tbalius/zzz.github/DOCK/docking/setup/setup_db2_zinc15_file_number.py ./ vs_frag /nfs/work/tbalius/database_ph4/frags.txt 500 count

The file /nfs/work/tbalius/database_ph4/frags.txt should be changed to that you downloaded from ZINC.

The above script has 5 parameters:
:(1) path where directories will be located (present directory);
:(2) prefix name of the directories;
:(3) the file that contains the db2 files locations;
:(4) the number of directories to be created; and
:(5) the type of run: count (evenly distributes the db2 file among the dirs, this is much faster than the other options), size (It will try and make the directory of equal size), or both (will try and satisfy both criteria).

Note that this script is avable in later beta versons of DOCK3.7.

===submitting the docking calculations===

This script will submit a job to the queue for each of the docking directorys created by the setup script.

$DOCKBASE/docking/submit/submit.csh

DOCK3.7 is a serial program and is parallelized by submiting many serial jobs to the queue.

===combining the results===

After your docking jobs have all completed, This script will combine all your results into an extract_all file.
$DOCKBASE/analysis/extract_all.py

This script will create a mol2 file with the top scoring molecules:
$DOCKBASE/analysis/getposes.py

=== curating and hit-picking ===

Typically we will visualize the ligands in UCSF Chimera using the veiwdock tool (using the DOCK4, 5, or 6 format).