{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:52:37,309 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but \"NUMEXPR_MAX_THREADS\" not set, so enforcing safe limit of 8.\n", "2024-12-05 13:52:37,310 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.\n", "/home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/pandas/core/computation/expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).\n", " from pandas.core.computation.check import NUMEXPR_INSTALLED\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. https://dx.doi.org/10.1021/acs.jctc.6b00049\n", "HTMD Documentation at: https://software.acellera.com/htmd/\n", "\n", "You are on the latest HTMD version (2.4.2+23.g4c98f4164.dirty).\n", "\n" ] } ], "source": [ "from htmd.ui import *\n", "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')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Using docking to generate starting poses for simulations\n", "\n", "Download the files for this tutorial from this [link](http://pub.htmd.org/nc983hu3brda/bentryp.tar.gz)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dock the protein with the ligand" ] }, { "cell_type": "code", "execution_count": 2, "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/tmpmhkzxhom/protein.pdb)\n", "\n", "1 molecule converted\n", "2024-12-05 13:52:42,850 - 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: 1238370264\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.9 0.042 1.620\n", " 3 -4.6 24.177 24.823\n", " 4 -4.5 30.437 31.604\n", " 5 -4.4 21.893 23.170\n", " 6 -4.3 28.748 29.753\n", " 7 -4.3 29.399 30.379\n", " 8 -4.3 2.472 3.609\n", " 9 -4.3 21.328 22.539\n", " 10 -4.3 21.329 22.516\n", " 11 -4.2 21.484 22.800\n", " 12 -4.1 29.504 30.504\n", " 13 -4.1 18.311 18.874\n", " 14 -4.1 24.407 25.335\n", " 15 -4.1 3.897 4.965\n", " 16 -4.0 24.067 25.048\n", " 17 -4.0 22.026 23.022\n", " 18 -3.8 22.494 23.397\n", " 19 -3.7 17.615 18.842\n", " 20 -3.7 30.998 32.036\n", "Writing output ... done.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "20 molecules converted\n", "20 files output. The first is /tmp/tmpmhkzxhom/output_1.pdb\n" ] } ], "source": [ "datadir = join(home(dataDir=\"building-protein-ligand\"))\n", "\n", "prot = Molecule(join(datadir, \"trypsin.pdb\"))\n", "prot.center()\n", "lig = Molecule(join(datadir, \"BEN.cif\"))\n", "poses, scores = dock(prot, lig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualize the docked poses" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "aaf9129de73f41b0aea7f3c380dc2467", "version_major": 2, "version_minor": 0 }, "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c88182225a074037be9a4e9988e3c13c", "version_major": 2, "version_minor": 0 }, "text/plain": [ "NGLWidget()" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mol = Molecule()\n", "mol.append(prot)\n", "for i, p in enumerate(poses):\n", " mol.append(p)\n", "mol.view(sel='protein', style='NewCartoon', hold=True)\n", "mol.view(sel='resname BEN', style='Licorice', color=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build systems from docked poses" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:53:44,922 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.\n", "2024-12-05 13:53:44,952 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb\n", "2024-12-05 13:53:45,252 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n", "\n", "2024-12-05 13:53:46,435 - htmd.builder.solvate - INFO - 19297 water molecules were added to the system.\n", "2024-12-05 13:53:48,964 - htmd.builder.amber - INFO - Detecting disulfide bonds.\n", "2024-12-05 13:53:48,981 - 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-12-05 13:53:50,234 - htmd.builder.amber - INFO - Starting the build.\n", "2024-12-05 13:53:53,664 - htmd.builder.amber - INFO - Finished building.\n", "2024-12-05 13:53:54,470 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n", "2024-12-05 13:53:59,698 - htmd.builder.ionize - INFO - Adding 9 anions + 0 cations for neutralizing and 108 ions for the given salt concentration 0.15 M.\n", "2024-12-05 13:54:11,829 - htmd.builder.amber - INFO - Starting the build.\n", "2024-12-05 13:54:14,825 - htmd.builder.amber - INFO - Finished building.\n", "2024-12-05 13:54:15,564 - 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-12-05 13:54:22,933 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224\n", "2024-12-05 13:54:23,071 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.\n", "2024-12-05 13:54:23,100 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb\n", "2024-12-05 13:54:23,346 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n", "\n", "2024-12-05 13:54:24,556 - htmd.builder.solvate - INFO - 19297 water molecules were added to the system.\n", "2024-12-05 13:54:27,067 - htmd.builder.amber - INFO - Detecting disulfide bonds.\n", "2024-12-05 13:54:27,075 - 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-12-05 13:54:28,367 - htmd.builder.amber - INFO - Starting the build.\n", "2024-12-05 13:54:31,500 - htmd.builder.amber - INFO - Finished building.\n", "2024-12-05 13:54:32,248 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n", "2024-12-05 13:54:37,561 - htmd.builder.ionize - INFO - Adding 9 anions + 0 cations for neutralizing and 108 ions for the given salt concentration 0.15 M.\n", "2024-12-05 13:54:49,812 - htmd.builder.amber - INFO - Starting the build.\n", "2024-12-05 13:54:52,936 - htmd.builder.amber - INFO - Finished building.\n", "2024-12-05 13:54:53,745 - moleculekit.writers - WARNING - Field \"resid\" of PDB overflows. Your data will be truncated to 4 characters.\n", "2024-12-05 13:55:01,315 - 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(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", " \n", " D = maxDistance(prot, 'all')\n", " \n", " ligand = p\n", " ligand.set('segid','L')\n", " \n", " mol = Molecule(name='combo')\n", " mol.append(prot)\n", " mol.append(ligand)\n", " \n", " D = D + 15\n", " smol = solvate(mol, minmax=[[-D, -D, -D], [D, D, D]])\n", " topos = [join(datadir, \"BEN.cif\")]\n", " params = [join(datadir, \"BEN.frcmod\")]\n", "\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Equilibrate the build systems" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:55:26,560 - acemd - INFO - Copying docked/build/1/structure.pdb to docked/equil/1/structure.pdb\n", "2024-12-05 13:55:26,564 - acemd - INFO - Copying docked/build/1/structure.prmtop to docked/equil/1/structure.prmtop\n", "2024-12-05 13:55:29,000 - acemd - INFO - Copying docked/build/2/structure.pdb to docked/equil/2/structure.pdb\n", "2024-12-05 13:55:29,003 - acemd - INFO - Copying docked/build/2/structure.prmtop to docked/equil/2/structure.prmtop\n" ] } ], "source": [ "from acemd.protocols import setup_equilibration\n", "\n", "builds = sorted(glob('docked/build/*/'))\n", "for i, bd in enumerate(builds):\n", " setup_equilibration(bd, f'docked/equil/{i+1}/', run=\"10ns\", temperature=298)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "scrolled": true }, "outputs": [], "source": [ "mdx = LocalGPUQueue()\n", "mdx.submit(glob('./docked/equil/*/'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "mdx.wait()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create the production folder" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "from acemd.protocols import setup_production\n", "\n", "equils = sorted(glob('docked/equil/*/'))\n", "for i, eq in enumerate(equils):\n", " setup_production(eq, f'docked/generators/{i+1}/', run=\"50ns\", temperature=298)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mdx = LocalGPUQueue()\n", "mdx.submit(glob('./docked/generators/*/'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mdx.wait()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.15" } }, "nbformat": 4, "nbformat_minor": 4 }