2023.05.08.6OIM_covalent - tbalius/teb_docking_test_sets GitHub Wiki

Covalent docking tutorial using KRAS in complex with Sotorasib/AMG 510 (PDB ID 6OIM)

Here we describe performing covalent docking with the following:

  • DOCK 3 (DOCKovalent)
  • DOCK 6 attach-and-grow algorithm (not released yet)
  • DOCK 6 HDB covalent (not released yet)

Software requirements

  • Python 2 environment with the following packages: Numpy
  • Python 3 environment with the following packages: Numpy, RDKit
  • AmberTools
  • ChemAxon, Corina, and Amsol 7.1 for database generation.
  • UCSF Chimera for visualization and running Dock Prep on the command line.
  • DOCK 3.7 for preparing the receptor using blastermaster (Python 2 program) and for docking.
  • DOCK 6.12 for database generation, calculating RMSD, and for docking with attach-and-grow (not released yet).

Organizing file directories

Directories for this tutorial:

6OIM_covalent_tutorial/
|-- dockprep/
|-- build_ligand_dock3/
|-- build_ligand_dock6/
|-- run_dock3/
|-- run_dock6/
+-- run_dock6_hdb/

Scripts on GitHub for this tutorial:

0001.processing.csh
0002.run_chimera_dockprep.csh
0002p5.fix_broken_h.csh
0003.run.antechamber.csh
0004.blastermaster.csh
0004.blastermaster_manualProt.csh
0005.convert_dock3_to_dock6.csh
0006.dock6_cov_sph.csh
setup.csh
calc_rmsd_dock3.csh
calc_rmsd_dock6.csh
mol2_python3.py*
mol_covalent_Si_to_Du_solv.py*
mol2_replace_sybyl_with_ele.py*
multimol2_removeH.py*

* available in teb_script_programs

Dock Prep

Prepare the system for docking in the dockprep directory.

Scripts are available here.

Download the Dock Prep scripts from GitHub by running the following wget commands:

wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0001.processing.csh
wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0002.run_chimera_dockprep.csh
wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0002p5.fix_broken_h.csh
wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0003.run.antechamber.csh
wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0004.blastermaster.csh
wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0004.blastermaster_manualProt.csh
wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0005.convert_dock3_to_dock6.csh
wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/0006.dock6_cov_sph.csh

Copy and paste the code block into a C shell script and run the script.

>> csh download_scripts.csh

You will need to modify several of the scripts.

If there is a call to a script not available in this repository, it is most likely available in teb_scripts_programs.

Prepare the receptor and cofactors

To prepare the receptor and cofactor files, run the following scripts in order:

>> csh 0001.processing.csh
>> csh 0002.run_chimera_dockprep.csh
>> csh 0002p5.fix_broken_h.csh
>> csh 0003.run.antechamber.csh

The 0001.processing.csh script downloads the 6OIM.pdb file from the web and splits the file into receptor (with magnesium ion), ligand, and cofactor files called rec.pdb, lig.pdb, and cof.pdb, respectively. The 0002.run_chimera_dockprep.csh script calls Chimera's Dock Prep and AddH functions on the receptor and cofactor files. The 0002p5.fix_broken_h.csh script fixes a protonation error made by Chimera's AddH by replacing the name and coordinates of the incorrect hydrogen in the cofactor file with the name and coordinates of the correct hydrogen. The 0003.run.antechamber.csh script runs Amber's AnteChamber to calculate charges.

Generate spheres and grids

Use blastermaster to generate DOCK 3 spheres and grids, run the following scripts in order:

>> csh 0004.blastermaster.csh
>> csh 0004.blastermaster_manualProt.csh
  • Note: you will need to wait for the blastermaster jobs to finish before proceeding

The 0004.blastermaster.csh script runs blastermaster on the receptor without the cofactor. The 0004.blastermaster_manualProt.csh script runs blastermaster on the receptor with the cofactor. Blastermaster is called with covalent flags.

Check the matching_spheres.sph file created by blastermaster in the blastermaster_cof/dockfiles directory, the last three lines should be formatted as follows:

 9001  -4.51400  -4.02200   2.39800   1.000    1 0  0
 9002  -5.43200  -4.54500   1.28800   1.000    1 0  0
 9003  -6.34400  -3.26000   0.40900   1.000    1 0  0

Blastermaster generates grids for DOCK 3. To convert the grids to DOCK 6 readable grids, run the 0005.convert_dock3_to_dock6.csh script.

>> csh 0005.convert_dock3_to_dock6.csh

To generate covalent spheres for covalent docking with DOCK 6, run the 0005.dock6_cov_sph.csh script.

>> csh 0006.dock6_cov_sph.csh

This script greps for the atom lines of the gamma sulfur, beta carbon, and alpha carbon of cysteine 12 from rec.pdb, saves them to a PDB file called cov_sph.pdb, and then calls pdbtosph generate the cov_sph.sph sphere file.

Check the cov_sph.pdb file inside the created dock6_cov_spheres directory, it should be formatted as follows:

ATOM     93  SG  CYS A  12      -6.344  -3.260   0.409  1.00 26.89           S
ATOM     92  CB  CYS A  12      -5.432  -4.545   1.288  1.00 24.31           C
ATOM     89  CA  CYS A  12      -4.514  -4.022   2.398  1.00 23.24           C

The DOCK 6 attach-and-grow covalent docking algorithm converts three atoms from the covalent cysteine sidechain into spheres for orienting: the alpha carbon, beta carbon, and gamma sulfur. The spheres must be in this specific order: 1) gamma 2) beta 3) alpha. This order is different from the covalent spheres file created by blastermaster used for covalent docking with DOCK 3.

To manually prepare covalent spheres for docking with DOCK 6, determine the gamma, beta, and alpha atoms of the covalent residue of interest, save the relevant atom lines to a pdb or mol2 file in the correct order, and convert the file to a sphere file using pdbtosph (distributed with DOCK 3) or mol2toSPH.py (available in teb_scripts_programs). This can also be applied to covalent residues other than cysteine (e.g. for serine use the gamma oxygen, beta carbon, and alpha carbon).

DOCK 3

Build ligand from SMILES

Database generation requires DOCK 6.12 and the following programs: ChemAxon, Corina, and Amsol 7.1.

Prepare the ligand database for covalent docking with DOCK 3 inside the build_ligand_dock3 directory.

Create a SMILES file called mov.smi containing the SMILES string of sotorasib.

O=C(N1CCN(C2=NC(N(c3c(C(C)C)nccc3C)c4nc(c5c(F)cccc5O)c(F)cc42)=O)[C@@H](C)C1)CC[SiH3] MOV

The SMILES has been modified for covalent docking by attaching a silicon atom where the molecule forms a covalent bond. This atom will be used to attach the ligand onto the receptor during docking.

RDKit can be used to process multiple SMILES by applying a reaction that appends the silicon atoms to your covalent warhead(s).

Sample DB2 generation scripts are distributed with DOCK 6.12. The sample DB2 generation pipeline requires the following programs: ChemAxon, Corina, and Amsol 7.1. ChemAxon's cxcalc is used for protonation, Corina is used to generate a 3d conformation of the ligand from SMILES, and Amsol 7.1 is used to calculate solvation free energies. DOCK 6 is used to generate multiple conformations of the ligand which are saved in a DB2 file. You will need to modify the scripts to run using the paths to your installed programs or to make program substitutions.

Run the ligand building script, the arguments passed are a SMILES file with your covalent ligand SMILES and a covalent flag.

>> csh ${dock6_path}/template_pipeline/hdb_lig_gen/generate/build_ligand_simple_with_dock6.csh mov.smi --covalent

  • ${dock6_path} is a variable (placeholder) for the DOCK 6 installation path (e.g. /home/baliuste/zzz.github/dock6)

After the script finishes running, check that DB2 files were generated inside the db_build_working directory. The resulting DB2 files are called output_anchorX_scored.db2.gz.

>> ls db_build_working/MOV/MOV_*_db2/*.db2.gz

Run covalent docking with DOCK 3

Run DOCK in the run_dock3 directory.

Download the setup script to run covalent docking with DOCK 3.

wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/setup.csh

This script creates two files: split_database_index and INDOCK. The split_database_index file contains the paths of the ligand database (.db2.gz) files prepared in the previous steps. The INDOCK file is an input file containing the parameters to run covalent docking with DOCK 3, and is created by copying and modifying the INDOCK file created by blastermaster. The script also creates a symbolic link to the dockfiles directory.

The created INDOCK file should be formatted as follows:

DOCK 3.7 parameter
#####################################################
# NOTE: split_database_index is reserved to specify a list of files
# defults for large scale docking.
ligand_atom_file               split_database_index
#####################################################
#                             OUTPUT
output_file_prefix            test.
#####################################################
#                             MATCHING
match_method                  2
distance_tolerance            0.05
match_goal                    1000
distance_step                 0.05
distance_maximum              0.5
timeout                       10.0
nodes_maximum                 4
nodes_minimum                 4
bump_maximum                  1000.0
bump_rigid                    1000.0
mol2_score_maximum            1000.0
#####################################################
#                             COLORING
chemical_matching             no
case_sensitive                no
#####################################################
#                             SEARCH MODE
atom_minimum                  4
atom_maximum                  25
number_save                   1
number_write                  1
flush_int                     100
#molecules_maximum            100000
check_clashes                 no
do_premax                     no
do_clusters                   no
#####################################################
#                             SCORING
ligand_desolvation            volume
#vdw_maximum                   1.0e10
ligand_desolv_scale           1.0
electrostatic_scale           1.0
vdw_scale                     1.0
internal_scale                0.0
per_atom_scores               no
#####################################################
#                             DOCKovalent
dockovalent                   yes
bond_len                      1.8
bond_ang1                     109.5
bond_ang2                     109.5
len_range                     0.0
len_step                      0.1
ang1_range                    10.0
ang2_range                    10.0
ang1_step                     2.5
ang2_step                     2.5
#####################################################
#                    MINIMIZATION
minimize                      no
sim_itmax                     500
sim_trnstep                   0.2
sim_rotstep                   5.0
sim_need_to_restart           1.0
sim_cnvrge                    0.1
min_cut                       1.0e15
iseed                         777
#####################################################
# INPUT FILES / THINGS THAT CHANGE
receptor_sphere_file          ./dockfiles/matching_spheres.sph
vdw_parameter_file            ./dockfiles/vdw.parms.amb.mindock
delphi_nsize                  99
flexible_receptor             no
total_receptors               1
############## grids/data for one receptor
rec_number                    1
rec_group                     1
rec_group_option              1
solvmap_file                  ./dockfiles/ligand.desolv.heavy
hydrogen_solvmap_file         ./dockfiles/ligand.desolv.hydrogen
delphi_file                   ./dockfiles/trim.electrostatics.phi
chemgrid_file                 ./dockfiles/vdw.vdw
bumpmap_file                  ./dockfiles/vdw.bmp
############## end of INDOCK

Run DOCK 3.

>> ${DOCKBASE}/bin/dock64

  • ${DOCKBASE} is a variable (placeholder) for the DOCK 3 installation path (e.g. /home/baliuste/zzz.github/DOCK)

Once DOCK 3 finishes running, unzip the output mol2 file.

>> gunzip test.mol2.gz

The unzipped mol2 file can be visualized in Chimera using ViewDock.

Calculate RMSD

Download the following scripts from the teb_scripts_programs repository:

wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/mol2_python3.py
wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/mol2_replace_sybyl_with_ele.py
wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/multimol2_removeH.py

The mol2_replace_sybyl_with_ele.py script takes a mol2 file as input and outputs a mol2 file with the sybyl atom types replaced with the element symbol. The multimol2_removeH.py removes the hydrogens from a mol2 file. The mol2_python.py script is required to run the other two scripts.

Download the calc_rmsd_dock3.csh script from this repository.

wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/calc_rmsd_dock3.csh

This script creates the files needed for calculating RMSD before calling DOCK 6 to calculate RMSD. First, the lig.pdb file created during Dock Prep is converted into a Mol2 file called ref.mol2 using Chimera. Then the mol2_replace_sybyl_with_ele.py script is run on the reference Mol2 file, ref.mol2, and the output file created by DOCK, test.mol2. The multimol2_removeH.py script is run to remove hydrogens from test.mol2. Finally, a DOCK 6 input file called rmsd.in is created and DOCK 6 is run to calculate RMSD.

Modify the script by changing the chimera path, the parameter file paths in the rmsd.in file, and the path to call DOCK 6 to your Chimera and DOCK 6 installation path.

Run the script.

>> csh calc_rmsd_dock3.csh

The rmsd_output_scored.mol2 file contains the poses with the corresponding RMSD values in the header. The Hungarian (symmetry-corrected) heavy-atom RMSD is labeled as HA_RMSDh.

##########                            HA_RMSDh:              10.216
##########                            HA_RMSDh:               5.045

DOCK 6

Build ligand from SMILES

Prepare the ligand database for covalent docking with DOCK 6 attach-and-grow inside the build_ligand_dock6 directory.

Create a SMILES file called mov.smi containing the SMILES string of sotorasib.

O=C(N1CCN(C2=NC(N(c3c(C(C)C)nccc3C)c4nc(c5c(F)cccc5O)c(F)cc42)=O)[C@@H](C)C1)CC[SiH2][SiH3] MOV

The SMILES has been modified for covalent docking by attaching two silicon atoms where the molecule forms a covalent bond. For covalent docking with DOCK 6 attach-and-grow, the ligand is prepared with two dummy atoms instead of one. The two dummy atoms are aligned to the gamma sulfur and beta carbon atoms of the covalent cysteine residue during docking.

RDKit can be used to process multiple SMILES by applying a reaction that appends the silicon atoms to your covalent warhead(s).

Run the ligand building script, the arguments passed are a SMILES file with your covalent ligand SMILES and a covalent flag.

>> csh ${dock6_path}/template_pipeline/hdb_lig_gen/generate/build_ligand_simple_with_dock6.csh mov.smi --covalent

Download a script from the teb_scripts_programs repository:

wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/mol_covalent_Si_to_Du_solv.py

This script converts silicon atoms to dummy atoms and appends the solvation parameters.

Paste the following code into a C shell script and run the script:

foreach mol2 (`ls db_build_working/*/*_output.mol2`)
  python mol_covalent_Si_to_Du_solv.py $mol2 ${mol2:r}.solv temp_Du 103 >> log
  cat temp_Du.mol2 >> ligand_cov.mol2
  rm temp_Du.mol2
end

>> csh mk_mol2_cov.csh

This will loop over each Mol2 file created by the ligand building script and run the mol_covalent_Si_to_Du_solv.py script, the output for each Mol2 file is appended to a file called ligand_cov.mol2. The ligand_cov.mol2 file will be used as input for covalent docking with DOCK 6.

Run covalent docking with DOCK 6 attach-and-grow

Create an input file for DOCK 6 called covalent.in and format it as follows:

conformer_search_type                                        covalent
pruning_use_clustering                                       yes
pruning_max_orients                                          1000
pruning_clustering_cutoff                                    100
covalent_bondlength                                          1.7:0.1:1.91
covalent_bondlength2                                         1.7:0.1:1.91
covalent_angle                                               103.0:1.0:105.0
covalent_dihedral_step                                       10.0
pruning_conformer_score_cutoff                               1000.0
use_clash_overlap                                            no
write_growth_tree                                            no
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
internal_energy_cutoff                                       100.0
ligand_atom_file                                             ../build_ligand_dock6/ligand_cov.mol2
limit_max_ligands                                            no
skip_molecule                                                no
read_mol_solvation                                           yes
calculate_rmsd                                               yes
use_rmsd_reference_mol                                       no
use_database_filter                                          no
orient_ligand                                                yes
automated_matching                                           yes
receptor_site_file                                           ../dockprep/dock6_cov_spheres/cov_sph.sph
max_orientations                                             1000
critical_points                                              no
chemical_matching                                            no
use_ligand_spheres                                           no
bump_filter                                                  no
score_molecules                                              yes
contact_score_primary                                        no
grid_score_primary                                           no
gist_score_primary                                           no
multigrid_score_primary                                      no
dock3.5_score_primary                                        yes
dock3.5_vdw_score                                            yes
dock3.5_grd_prefix                                           ../dockprep/blastermaster_cof/dock6files/chem52
dock3.5_electrostatic_score                                  yes
dock3.5_ligand_desolvation_score                             volume
dock3.5_solvent_occlusion_file                               ../dockprep/blastermaster_cof/dock6files/ligand.desolv.heavy
dock3.5_redistribute_positive_desolvation                    no
dock3.5_hydrogen_desolvation_grid                            yes
dock3.5_hydrogen_solvent_occlusion_file                      ../dockprep/blastermaster_cof/dock6files/ligand.desolv.hydrogen
dock3.5_receptor_desolvation_score                           no
dock3.5_write_atomic_energy_contrib                          no
dock3.5_score_vdw_scale                                      1
dock3.5_score_es_scale                                       1
minimize_ligand                                              yes
minimize_anchor                                              no
minimize_flexible_growth                                     yes
use_advanced_simplex_parameters                              no
minimize_flexible_growth_ramp                                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_grow_max_iterations                                  0
simplex_grow_tors_premin_iterations                          500
simplex_final_min                                            no
simplex_random_seed                                          0
simplex_restraint_min                                        no
atom_model                                                   all
vdw_defn_file                                                ${dock6_path}/parameters/vdw_AMBER_parm94.dock3_7.defn
flex_defn_file                                               ${dock6_path}/parameters/flex.defn
flex_drive_file                                              ${dock6_path}/parameters/flex_drive.tbl
ligand_outfile_prefix                                        covalent.out
write_mol_solvation                                          no
write_orientations                                           no
num_scored_conformers                                        10
write_conformations                                          no
cluster_conformations                                        yes
cluster_rmsd_threshold                                       0.5
score_threshold                                              100.0
rank_ligands                                                 no

Replace ${dock6_path} in the parameter file lines with the DOCK 6 installation path.

The conformer_search_type input parameter is set to covalent, and four parameters to sample the covalent bond environment are introduced: covalent_bondlength, covalent_bondlength2, covalent_angle, and covalent_dihedral_step.

Here is a diagram of the covalent bond environment and values being sampled:

The covalent_bondlength, covalent_bondlength2, and covalent_angle parameters can given a single value or multiple values to iterate over using the form start:step:stop. In the example parameter file given above, the covalent_bondlength and covalent_bondlength2 parameter is set to 1.7:0.1:1.91, so bond lengths of 1.7, 1.8, and 1.9 will be sampled during docking for the two bonds. The covalent_angle parameter is set to 103.0:1.0:105.0, so covalent angles of 103.0, 104.0, and 105.0 will be sampled during docking. In total, the molecule will be docked in 27 different covalent bond environments. At the time this tutorial was written (November 2023), this docking took ~14 minutes (~833 seconds); when docking with a single bond environment, the docking took ~27 seconds.

Run DOCK 6.

>> ${dock6_path}/bin/dock6 -i covalent.in -o covalent.out -v

The covalent.out_scored.mol2 file can be visualized in Chimera using ViewDock.

The output Mol2 file is written with the dummy atoms still attached to the ligand. The two dummy atoms should overlap with the gamma sulfur and beta carbon of cysteine 12.

Calculate RMSD

Download the following scripts from the teb_scripts_programs repository:

wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/mol2_python3.py
wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/mol2_replace_sybyl_with_ele.py
wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/multimol2_removeH.py

The mol2_replace_sybyl_with_ele.py script takes a mol2 file as input and outputs a mol2 file with the sybyl atom types replaced with the element symbol. The multimol2_removeH.py removes the hydrogens from a mol2 file. The mol2_python.py script is required to run the other two scripts. These are the same scripts used to prepare Mol2 files for RMSD calculation for covalent docking with DOCK 3.

Download the calc_rmsd_dock6.csh script from this repository.

wget https://raw.githubusercontent.com/tbalius/teb_docking_test_sets/main/scripts_for_tutorial/scripts_for_2023.05.08.6OIM_covalent/calc_rmsd_dock6.csh

This script creates the files needed for calculating RMSD before calling DOCK 6 to calculate RMSD. First, the lig.pdb file created during Dock Prep is converted into a Mol2 file called ref.mol2 using Chimera. The covalent.out_scored.mol2 is opened in Chimera and the dummy atoms are deleted. Then the mol2_replace_sybyl_with_ele.py script is run on the reference Mol2 file, ref.mol2, and the output file created by DOCK, covalent.out_scored.mol2. The multimol2_removeH.py script is run to remove hydrogens. Finally, a DOCK 6 input file called rmsd.in is created and DOCK 6 is run to calculate RMSD.

Modify the script by changing the chimera path, the parameter file paths in the rmsd.in file, and the path to call DOCK 6 to your Chimera and DOCK 6 installation paths.

Run the script.

>> csh calc_rmsd_dock6.csh

The rmsd_output_scored.mol2 file contains the poses with the corresponding RMSD values in the header. The Hungarian (symmetry-corrected) heavy-atom RMSD is labeled as HA_RMSDh.

##########                            HA_RMSDh:               0.749
##########                            HA_RMSDh:               1.275
##########                            HA_RMSDh:               0.696

The values for the bond environment from the crystal structure (PDB ID 6OIM) are d1 = 1.804331 Å, d2 = 1.805427 Å, and θ = 116.120121°. Even though we do not sample the exact bond angle, DOCK can still reproduce the crystallographic ligand pose.

Process SMILES with RDKit

We use RDKit's rdkit.Chem.rdChemReactions module to process multiple SMILES.

Download a script from the teb_scripts_programs repository:

wget https://raw.githubusercontent.com/tbalius/teb_scripts_programs/tree/master/zzz.scripts/simple_reaction_file.py

This script will apply a SMARTS reaction to each of your SMILES.

Sotorasib is used as an example:

O=C(N1CCN(C2=NC(N(c3c(C(C)C)nccc3C)c4nc(c5c(F)cccc5O)c(F)cc42)=O)[C@@H](C)C1)C=C MOV

Sotorasib has an acrylamide covalent warhead.

To append one silicon atom:

>> python simple_reaction_file.py "[SiH3:1].[C:2]=[C:3][C:4]=[O:5]>>[SiH3:1][C:2][C:3][C:4]=[O:5]" "[SiH3]" mov.smi mov_si.smi

To append two silicon atoms:

>> python simple_reaction_file.py "[SiH3:1][SiH2:2].[C:3]=[C:4][C:5]=[O:6]>>[SiH3:1][SiH2:2][C:3][C:4][C:5]=[O:6]" "[SiH2][SiH3]" mov.smi mov_si.smi

Change the reaction to target different covalent warheads. For vinyl sulfonamides:

>> python simple_reaction_file.py "[SiH3:1][SiH2:2].[C:3]=[C:4][S:5](=[O:6])=[O:7]>>[SiH3:1][SiH2:2][C:3][C:4][S:5](=[O:6])=[O:7]" "[SiH2][SiH3]" input.smi output.smi

Run covalent docking with DOCK 6 HDB covalent

We can also dock covalent DB2 files using DOCK 6 HDB covalent.

To prepare the ligand database, see the DOCK 3 section.

Copy or regenerate the split_database_index file used for DOCK 3.

../build_ligand_dock3/db_build_working/MOV/MOV_0_db2/output_anchor1_scored.db2.gz
../build_ligand_dock3/db_build_working/MOV/MOV_1_db2/output_anchor1_scored.db2.gz

Create an input file for DOCK 6 called covalent.in and format it as follows:

conformer_search_type                                        HDB_covalent
num_per_search                                               10
skip_broken                                                  no
hdb_db2_input_file                                           split_database_index
hdb_db2_search_score_threshold                               1000.0
hdb_db2_search_covalent_bondlength1                          1.8
hdb_db2_search_covalent_bondlength2                          1.8
hdb_db2_search_covalent_bondlength3                          1.8
hdb_db2_search_covalent_angle1                               100.0
hdb_db2_search_covalent_angle2                               100.0
hdb_db2_search_covalent_angle3                               100.0
hdb_db2_search_covalent_dihederal_step1                      10.0
hdb_db2_search_covalent_dihederal_step2                      10.0
hdb_write_dummy_mol2                                         yes
use_internal_energy                                          yes
internal_energy_rep_exp                                      12
internal_energy_cutoff                                       100.0
use_database_filter                                          no
orient_ligand                                                yes
automated_matching                                           yes
receptor_site_file                                           ../dockprep/dock6_cov_spheres/cov_sph.sph
max_orientations                                             1000
critical_points                                              no
chemical_matching                                            no
use_ligand_spheres                                           no
bump_filter                                                  no
score_molecules                                              yes
contact_score_primary                                        no
grid_score_primary                                           no
gist_score_primary                                           no
DistJoin_score_primary                                       no
multigrid_score_primary                                      no
dock3.5_score_primary                                        yes
dock3.5_vdw_score                                            yes
dock3.5_grd_prefix                                           ../dockprep/blastermaster_cof/dock6files/chem52
dock3.5_electrostatic_score                                  yes
dock3.5_ligand_desolvation_score                             volume
dock3.5_solvent_occlusion_file                               ../dockprep/blastermaster_cof/dock6files/ligand.desolv.heavy
dock3.5_redistribute_positive_desolvation                    no
dock3.5_hydrogen_desolvation_grid                            yes
dock3.5_hydrogen_solvent_occlusion_file                      ../dockprep/blastermaster_cof/dock6files/ligand.desolv.hydrogen
dock3.5_receptor_desolvation_score                           no
dock3.5_write_atomic_energy_contrib                          no
dock3.5_score_vdw_scale                                      1
dock3.5_score_es_scale                                       1
minimize_ligand                                              no
simplex_final_min                                            no
atom_model                                                   all
vdw_defn_file                                                ${dock6_path}/parameters/vdw_AMBER_parm94.dock3_7.defn
flex_defn_file                                               ${dock6_path}/parameters/flex.defn
flex_drive_file                                              ${dock6_path}/parameters/flex_drive.tbl
ligand_outfile_prefix                                        covalent.out
write_mol_solvation                                          no
write_orientations                                           no
num_scored_conformers                                        10
write_conformations                                          no
cluster_conformations                                        yes
cluster_rmsd_threshold                                       0.5
score_threshold                                              100.0
rank_ligands                                                 no

Run DOCK 6.

>> ${dock6_path}/bin/dock6 -i covalent.in -o covalent.out -v

The output Mol2 is written with one dummy atom attached to the ligand. Setting hdb_write_dummy_mol2 to no will tell DOCK to not write out the dummy atom in the output, which can be helpful for calculating RMSD.

The calc_rmsd_dock6.csh script used to calculate RMSD for attach-and-grow can be also used for HDB covalent.