from htmd.ui import *
from moleculekit.config import config
from moleculekit.util import maxDistance
from htmd.home import home
from os.path import join
config(viewer='ngl')
Using docking to generate starting poses for simulations#
Download the files for this tutorial from this link
Dock the protein with the ligand#
datadir = join(home(dataDir="building-protein-ligand"))
prot = Molecule(join(datadir, "trypsin.pdb"))
prot.center()
lig = Molecule(join(datadir, "BEN.cif"))
poses, scores = dock(prot, lig)
==============================
*** Open Babel Warning in PerceiveBondOrders
Failed to kekulize aromatic bonds in OBMol::PerceiveBondOrders (title is /tmp/tmpai1237b_/protein.pdb)
1 molecule converted
2024-06-12 13:30:30,555 - htmd.dock - INFO - Charges detected in ligand and will be used for docking.
1 molecule converted
################################################################# # If you used AutoDock Vina in your work, please cite: # # # # O. Trott, A. J. Olson, # # AutoDock Vina: improving the speed and accuracy of docking # # with a new scoring function, efficient optimization and # # multithreading, Journal of Computational Chemistry 31 (2010) # # 455-461 # # # # DOI 10.1002/jcc.21334 # # # # Please see http://vina.scripps.edu for more information. # ################################################################# WARNING: The search space volume > 27000 Angstrom^3 (See FAQ) Detected 20 CPUs WARNING: at low exhaustiveness, it may be impossible to utilize all CPUs Reading input ... done. Setting up the scoring function ... done. Analyzing the binding site ... done. Using random seed: -1777289594 Performing search ... 0% 10 20 30 40 50 60 70 80 90 100% |----|----|----|----|----|----|----|----|----|----| *********************************************** done. Refining results ... done. mode | affinity | dist from best mode | (kcal/mol) | rmsd l.b.| rmsd u.b. -----+------------+----------+---------- 1 -4.9 0.000 0.000 2 -4.6 24.178 24.846 3 -4.5 30.431 31.607 4 -4.4 21.896 23.222 5 -4.3 28.747 29.751 6 -4.3 24.124 24.750 7 -4.3 29.405 30.422 8 -4.3 29.932 30.983 9 -4.2 21.067 22.375 10 -4.1 18.337 18.880 11 -4.1 24.363 25.414 12 -4.1 18.396 19.036 13 -4.0 31.747 32.325 14 -3.9 19.767 20.688 15 -3.8 23.783 24.742 16 -3.6 32.220 33.387 17 -3.5 28.693 29.262 18 -3.5 29.433 30.277 19 -3.5 31.659 32.851 20 -3.5 29.840 30.736 Writing output ... done.
20 molecules converted
20 files output. The first is /tmp/tmpai1237b_/output_1.pdb
Visualize the docked poses#
mol = Molecule()
mol.append(prot)
for i, p in enumerate(poses):
mol.append(p)
mol.view(sel='protein', style='NewCartoon', hold=True)
mol.view(sel='resname MOL', style='Licorice', color=1)
NGLWidget()
Build systems from docked poses#
molbuilt = []
for i, p in enumerate(poses):
prot = Molecule(join(datadir, "trypsin.pdb"))
prot.filter('chain A and (protein or water or resname CA)')
prot.set('segid', 'P', sel='protein and noh')
prot.set('segid', 'W', sel='water')
prot.set('segid', 'CA', sel='resname CA')
prot.center()
D = maxDistance(prot, 'all')
ligand = p
ligand.set('segid','L')
ligand.set('resname','BEN')
mol = Molecule(name='combo')
mol.append(prot)
mol.append(ligand)
D = D + 15
smol = solvate(mol, minmax=[[-D, -D, -D], [D, D, D]])
topos = [join(datadir, "BEN.cif")]
params = [join(datadir, "BEN.frcmod")]
molbuilt.append(amber.build(smol, topo=topos, param=params, outdir=f'./docked/build/{i+1}/', saltconc=0.15))
if i==1: # For demonstration purposes lets only build the two first
break
2024-06-12 13:32:26,799 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.
2024-06-12 13:32:26,828 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-06-12 13:32:27,147 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
Solvating: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:01<00:00, 4.19it/s]
2024-06-12 13:32:29,547 - htmd.builder.solvate - INFO - 19297 water molecules were added to the system.
2024-06-12 13:32:31,741 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-06-12 13:32:31,748 - htmd.builder.builder - INFO - 6 disulfide bonds were added
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 26, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 42, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 117, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 184, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 174, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 198, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 149, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 163, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 8, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 138, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 110, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 211, insertion: '', segid: 'P'>
2024-06-12 13:32:32,989 - htmd.builder.amber - INFO - Starting the build.
2024-06-12 13:32:36,048 - htmd.builder.amber - INFO - Finished building.
2024-06-12 13:32:36,844 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-06-12 13:32:41,854 - htmd.builder.ionize - INFO - Adding 9 anions + 0 cations for neutralizing and 108 ions for the given salt concentration 0.15 M.
2024-06-12 13:32:53,962 - htmd.builder.amber - INFO - Starting the build.
2024-06-12 13:32:57,244 - htmd.builder.amber - INFO - Finished building.
2024-06-12 13:32:58,091 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
/home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/Bio/pairwise2.py:278: BiopythonDeprecationWarning: Bio.pairwise2 has been deprecated, and we intend to remove it in a future release of Biopython. As an alternative, please consider using Bio.Align.PairwiseAligner as a replacement, and contact the Biopython developers if you still need the Bio.pairwise2 module.
warnings.warn(
2024-06-12 13:33:05,155 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224
2024-06-12 13:33:05,300 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.
2024-06-12 13:33:05,349 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-06-12 13:33:05,599 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
Solvating: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:01<00:00, 4.34it/s]
2024-06-12 13:33:07,936 - htmd.builder.solvate - INFO - 19296 water molecules were added to the system.
2024-06-12 13:33:10,205 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-06-12 13:33:10,212 - htmd.builder.builder - INFO - 6 disulfide bonds were added
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 26, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 42, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 117, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 184, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 174, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 198, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 149, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 163, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 8, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 138, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 110, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 211, insertion: '', segid: 'P'>
2024-06-12 13:33:11,436 - htmd.builder.amber - INFO - Starting the build.
2024-06-12 13:33:14,563 - htmd.builder.amber - INFO - Finished building.
2024-06-12 13:33:15,360 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-06-12 13:33:20,275 - htmd.builder.ionize - INFO - Adding 9 anions + 0 cations for neutralizing and 108 ions for the given salt concentration 0.15 M.
2024-06-12 13:33:32,260 - htmd.builder.amber - INFO - Starting the build.
2024-06-12 13:33:35,364 - htmd.builder.amber - INFO - Finished building.
2024-06-12 13:33:36,201 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-06-12 13:33:43,212 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224
Equilibrate the build systems#
from htmd.protocols.equilibration_v3 import Equilibration
md = Equilibration()
md.runtime = 10
md.timeunits = "ns"
md.temperature = 298
builds = sorted(glob('docked/build/*/'))
for i, b in enumerate(builds):
md.write(b, f'docked/equil/{i+1}/')
mdx = LocalGPUQueue()
mdx.submit(glob('./docked/equil/*/'))
mdx.wait()
Create the production folder#
from htmd.protocols.production_v6 import Production
md = Production()
md.runtime = 50
md.timeunits = 'ns'
md.temperature = 298
equils = sorted(glob('docked/equil/*/'))
for i, b in enumerate(equils):
md.write(b, f'docked/generators/{i+1}/')
mdx = LocalGPUQueue()
mdx.submit(glob('./docked/generators/*/'))
mdx.wait()