From 8b8902599ca20940fc2161a90e1e1e24a530d6bb Mon Sep 17 00:00:00 2001 From: Stefan Doerr Date: Wed, 12 Jun 2024 13:36:07 +0300 Subject: [PATCH] updated docking tutorial --- tutorials/docking-simulation-generators.ipynb | 279 ++++++++++++++++-- 1 file changed, 248 insertions(+), 31 deletions(-) diff --git a/tutorials/docking-simulation-generators.ipynb b/tutorials/docking-simulation-generators.ipynb index 012fdfc89..90de183fa 100644 --- a/tutorials/docking-simulation-generators.ipynb +++ b/tutorials/docking-simulation-generators.ipynb @@ -2,8 +2,9 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": { + "collapsed": false, "jupyter": { "outputs_hidden": false } @@ -11,8 +12,12 @@ "outputs": [], "source": [ "from htmd.ui import *\n", - "config(viewer='ngl')\n", - "os.chdir('/webdata/nc983hu3brda/') # Don't use this command" + "from moleculekit.config import config\n", + "from moleculekit.util import maxDistance\n", + "from htmd.home import home\n", + "from os.path import join\n", + "\n", + "config(viewer='ngl')" ] }, { @@ -33,18 +38,95 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "==============================\n", + "*** Open Babel Warning in PerceiveBondOrders\n", + " Failed to kekulize aromatic bonds in OBMol::PerceiveBondOrders (title is /tmp/tmpai1237b_/protein.pdb)\n", + "\n", + "1 molecule converted\n", + "2024-06-12 13:30:30,555 - htmd.dock - INFO - Charges detected in ligand and will be used for docking.\n", + "1 molecule converted\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "#################################################################\n", + "# If you used AutoDock Vina in your work, please cite: #\n", + "# #\n", + "# O. Trott, A. J. Olson, #\n", + "# AutoDock Vina: improving the speed and accuracy of docking #\n", + "# with a new scoring function, efficient optimization and #\n", + "# multithreading, Journal of Computational Chemistry 31 (2010) #\n", + "# 455-461 #\n", + "# #\n", + "# DOI 10.1002/jcc.21334 #\n", + "# #\n", + "# Please see http://vina.scripps.edu for more information. #\n", + "#################################################################\n", + "\n", + "WARNING: The search space volume > 27000 Angstrom^3 (See FAQ)\n", + "Detected 20 CPUs\n", + "WARNING: at low exhaustiveness, it may be impossible to utilize all CPUs\n", + "Reading input ... done.\n", + "Setting up the scoring function ... done.\n", + "Analyzing the binding site ... done.\n", + "Using random seed: -1777289594\n", + "Performing search ... \n", + "0% 10 20 30 40 50 60 70 80 90 100%\n", + "|----|----|----|----|----|----|----|----|----|----|\n", + "***************************************************\n", + "done.\n", + "Refining results ... done.\n", + "\n", + "mode | affinity | dist from best mode\n", + " | (kcal/mol) | rmsd l.b.| rmsd u.b.\n", + "-----+------------+----------+----------\n", + " 1 -4.9 0.000 0.000\n", + " 2 -4.6 24.178 24.846\n", + " 3 -4.5 30.431 31.607\n", + " 4 -4.4 21.896 23.222\n", + " 5 -4.3 28.747 29.751\n", + " 6 -4.3 24.124 24.750\n", + " 7 -4.3 29.405 30.422\n", + " 8 -4.3 29.932 30.983\n", + " 9 -4.2 21.067 22.375\n", + " 10 -4.1 18.337 18.880\n", + " 11 -4.1 24.363 25.414\n", + " 12 -4.1 18.396 19.036\n", + " 13 -4.0 31.747 32.325\n", + " 14 -3.9 19.767 20.688\n", + " 15 -3.8 23.783 24.742\n", + " 16 -3.6 32.220 33.387\n", + " 17 -3.5 28.693 29.262\n", + " 18 -3.5 29.433 30.277\n", + " 19 -3.5 31.659 32.851\n", + " 20 -3.5 29.840 30.736\n", + "Writing output ... done.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "20 molecules converted\n", + "20 files output. The first is /tmp/tmpai1237b_/output_1.pdb\n" + ] } - }, - "outputs": [], + ], "source": [ - "prot = Molecule('bentryp/trypsin.pdb')\n", + "datadir = join(home(dataDir=\"building-protein-ligand\"))\n", + "\n", + "prot = Molecule(join(datadir, \"trypsin.pdb\"))\n", "prot.center()\n", - "lig = Molecule('bentryp/benzamidine.pdb')\n", + "lig = Molecule(join(datadir, \"BEN.cif\"))\n", "poses, scores = dock(prot, lig)" ] }, @@ -57,13 +139,41 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": { + "collapsed": false, "jupyter": { "outputs_hidden": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "23c06a4888f94844875cf2efa1e7ce97", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "4b26c51a2fa64d1f88c5b6832e89b546", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "NGLWidget()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "mol = Molecule()\n", "mol.append(prot)\n", @@ -82,29 +192,130 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": { + "collapsed": false, "jupyter": { "outputs_hidden": false }, "scrolled": true }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2024-06-12 13:32:26,799 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.\n", + "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\n", + "2024-06-12 13:32:27,147 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n", + "Solvating: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:01<00:00, 4.19it/s]\n", + "2024-06-12 13:32:29,547 - htmd.builder.solvate - INFO - 19297 water molecules were added to the system.\n", + "2024-06-12 13:32:31,741 - htmd.builder.amber - INFO - Detecting disulfide bonds.\n", + "2024-06-12 13:32:31,748 - htmd.builder.builder - INFO - 6 disulfide bonds were added\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2024-06-12 13:32:32,989 - htmd.builder.amber - INFO - Starting the build.\n", + "2024-06-12 13:32:36,048 - htmd.builder.amber - INFO - Finished building.\n", + "2024-06-12 13:32:36,844 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n", + "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.\n", + "2024-06-12 13:32:53,962 - htmd.builder.amber - INFO - Starting the build.\n", + "2024-06-12 13:32:57,244 - htmd.builder.amber - INFO - Finished building.\n", + "2024-06-12 13:32:58,091 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n", + "/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.\n", + " warnings.warn(\n", + "2024-06-12 13:33:05,155 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224\n", + "2024-06-12 13:33:05,300 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.\n", + "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\n", + "2024-06-12 13:33:05,599 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n", + "Solvating: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:01<00:00, 4.34it/s]\n", + "2024-06-12 13:33:07,936 - htmd.builder.solvate - INFO - 19296 water molecules were added to the system.\n", + "2024-06-12 13:33:10,205 - htmd.builder.amber - INFO - Detecting disulfide bonds.\n", + "2024-06-12 13:33:10,212 - htmd.builder.builder - INFO - 6 disulfide bonds were added\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n", + "Disulfide Bond between: UniqueResidueID\n", + " and: UniqueResidueID\n", + "\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2024-06-12 13:33:11,436 - htmd.builder.amber - INFO - Starting the build.\n", + "2024-06-12 13:33:14,563 - htmd.builder.amber - INFO - Finished building.\n", + "2024-06-12 13:33:15,360 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n", + "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.\n", + "2024-06-12 13:33:32,260 - htmd.builder.amber - INFO - Starting the build.\n", + "2024-06-12 13:33:35,364 - htmd.builder.amber - INFO - Finished building.\n", + "2024-06-12 13:33:36,201 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n", + "2024-06-12 13:33:43,212 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224\n" + ] + } + ], "source": [ "molbuilt = []\n", + "\n", "for i, p in enumerate(poses):\n", - " prot = Molecule('bentryp/trypsin.pdb')\n", + " prot = Molecule(join(datadir, \"trypsin.pdb\"))\n", " prot.filter('chain A and (protein or water or resname CA)')\n", " prot.set('segid', 'P', sel='protein and noh')\n", " prot.set('segid', 'W', sel='water')\n", " prot.set('segid', 'CA', sel='resname CA')\n", " prot.center()\n", - " from moleculekit.util import maxDistance\n", + " \n", " D = maxDistance(prot, 'all')\n", " \n", " ligand = p\n", " ligand.set('segid','L')\n", - " ligand.set('resname','MOL')\n", + " ligand.set('resname','BEN')\n", " \n", " mol = Molecule(name='combo')\n", " mol.append(prot)\n", @@ -112,11 +323,11 @@ " \n", " D = D + 15\n", " smol = solvate(mol, minmax=[[-D, -D, -D], [D, D, D]])\n", - " topos = ['top/top_all22star_prot.rtf', 'top/top_water_ions.rtf', './bentryp/benzamidine.rtf']\n", - " params = ['par/par_all22star_prot.prm', 'par/par_water_ions.prm', './bentryp/benzamidine.prm']\n", + " topos = [join(datadir, \"BEN.cif\")]\n", + " params = [join(datadir, \"BEN.frcmod\")]\n", "\n", - " molbuilt.append(charmm.build(smol, topo=topos, param=params, outdir='./docked/build/{}/'.format(i+1), saltconc=0.15))\n", - " if i==1: # For time purposes lets only build the two first\n", + " molbuilt.append(amber.build(smol, topo=topos, param=params, outdir=f'./docked/build/{i+1}/', saltconc=0.15))\n", + " if i==1: # For demonstration purposes lets only build the two first\n", " break" ] }, @@ -129,8 +340,9 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": { + "collapsed": false, "jupyter": { "outputs_hidden": false } @@ -138,19 +350,22 @@ "outputs": [], "source": [ "from htmd.protocols.equilibration_v3 import Equilibration\n", + "\n", "md = Equilibration()\n", - "md.numsteps = 1000\n", + "md.runtime = 10\n", + "md.timeunits = \"ns\"\n", "md.temperature = 298\n", "\n", - "builds = natsort(glob('docked/build/*/'))\n", + "builds = sorted(glob('docked/build/*/'))\n", "for i, b in enumerate(builds):\n", - " md.write(b, 'docked/equil/{}/'.format(i+1))" + " md.write(b, f'docked/equil/{i+1}/')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "collapsed": false, "jupyter": { "outputs_hidden": false }, @@ -166,6 +381,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "collapsed": false, "jupyter": { "outputs_hidden": false } @@ -186,6 +402,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "collapsed": false, "jupyter": { "outputs_hidden": false } @@ -198,9 +415,9 @@ "md.timeunits = 'ns'\n", "md.temperature = 298\n", "\n", - "equils = natsort(glob('docked/equil/*/'))\n", + "equils = sorted(glob('docked/equil/*/'))\n", "for i, b in enumerate(equils):\n", - " md.write(b, 'docked/generators/{}/'.format(i+1))" + " md.write(b, f'docked/generators/{i+1}/')" ] }, { @@ -235,7 +452,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -249,7 +466,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.10.13" } }, "nbformat": 4,