G-rdkit-ligand-modification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2492b2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -q https://raw.githubusercontent.com/openforcefield/2023-workshop-vignettes/master/colab_setup.ipynb\n",
    "%run colab_setup.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d826f45",
   "metadata": {},
   "outputs": [],
   "source": [
    "%env LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$CONDA_PREFIX/lib/"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "461250f5",
   "metadata": {},
   "source": [
    "# OpenFF Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf0b4910",
   "metadata": {},
   "outputs": [],
   "source": [
    "from openff.toolkit import ForceField, Molecule, Topology\n",
    "\n",
    "orig_ligand = Molecule.from_file(\"inputs/PT2385.sdf\")\n",
    "orig_top = Topology.from_pdb(\n",
    "    \"inputs/solvated_complex.pdb\", unique_molecules=[orig_ligand]\n",
    ")\n",
    "ff = ForceField(\"openff-2.0.0.offxml\", \"ff14sb_off_impropers_0.0.3.offxml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07b32834",
   "metadata": {},
   "outputs": [],
   "source": [
    "orig_ligand.visualize(backend=\"rdkit\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bebac5bf",
   "metadata": {},
   "source": [
    "# RDKit Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa9e07ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use RDKit's reaction handling to try mutating each aliphatic C-H bond to C-F\n",
    "import rdkit\n",
    "\n",
    "rdmol = orig_ligand.to_rdkit()\n",
    "rxn = rdkit.Chem.rdChemReactions.ReactionFromSmarts(\"[C:1][H:2] >> [C:1][F:2]\")\n",
    "products = rxn.RunReactants([rdmol])\n",
    "sanitized_products = []\n",
    "# mols_for_drawing = [rdmol]\n",
    "\n",
    "for product in products:\n",
    "    mol_copy = rdkit.Chem.Mol(product[0])\n",
    "    rdkit.Chem.SanitizeMol(mol_copy, rdkit.Chem.SANITIZE_ALL)\n",
    "    rdkit.Chem.AssignStereochemistryFrom3D(mol_copy)\n",
    "    rdkit.Chem.Kekulize(mol_copy, clearAromaticFlags=True)\n",
    "    rdkit.Chem.SetAromaticity(mol_copy, rdkit.Chem.AromaticityModel.AROMATICITY_MDL)\n",
    "    sanitized_products.append(mol_copy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "396bf766",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since some substitutions won't actually lead to a \"new\" molecule, filter by isomeric SMILES\n",
    "unique_offmols = [orig_ligand]\n",
    "rdmols_for_drawing = [rdmol]\n",
    "\n",
    "seen_smiles = set()\n",
    "for san_prod in sanitized_products:\n",
    "    offmol = Molecule.from_rdkit(san_prod)\n",
    "    smi = offmol.to_smiles()\n",
    "    if smi not in seen_smiles:\n",
    "        print(f\"smiles is unique: {smi}\")\n",
    "        unique_offmols.append(offmol)\n",
    "        seen_smiles.add(smi)\n",
    "        rdmols_for_drawing.append(san_prod)\n",
    "\n",
    "    else:\n",
    "        print(f\"smiles is not unique: {smi}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c39cc7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare rdmols molecules to be drawn as a grid image.\n",
    "# I don't like the 3D projections so I remove the conformers here.\n",
    "\n",
    "for rdmol_for_drawing in rdmols_for_drawing:\n",
    "    rdkit.Chem.rdchem.Mol.RemoveAllConformers(rdmol_for_drawing)\n",
    "rdkit.Chem.Draw.MolsToGridImage(rdmols_for_drawing, molsPerRow=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "457a114d",
   "metadata": {},
   "source": [
    "# OpenFF and OpenMM Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59558d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openmm\n",
    "from openff.units import Quantity, unit\n",
    "from openmm import unit as openmm_unit\n",
    "\n",
    "idx_to_sim = dict()\n",
    "n_systems = len(unique_offmols)\n",
    "\n",
    "for idx, unique_offmol in enumerate(unique_offmols):\n",
    "    print(f\"Making topology {idx+1}/{n_systems}\")\n",
    "\n",
    "    this_top = Topology.from_molecules([*orig_top.molecules][:-1] + [unique_offmol])\n",
    "    this_top.box_vectors = orig_top.box_vectors\n",
    "    pdb_filename = f\"topology{idx+1}.pdb\"\n",
    "    this_top.to_file(pdb_filename)\n",
    "\n",
    "    print(f\"Parametrizing system {idx+1}/{n_systems}\")\n",
    "\n",
    "    sys = ff.create_openmm_system(this_top)\n",
    "\n",
    "    integrator = openmm.LangevinIntegrator(\n",
    "        300 * openmm_unit.kelvin,\n",
    "        1 / openmm_unit.picosecond,\n",
    "        0.002 * openmm_unit.picoseconds,\n",
    "    )\n",
    "\n",
    "    print(f\"Minimizing system {idx+1}/{n_systems}\")\n",
    "\n",
    "    simulation = openmm.app.Simulation(this_top.to_openmm(), sys, integrator)\n",
    "    simulation.context.setPositions(this_top.get_positions().to_openmm())\n",
    "    simulation.minimizeEnergy(\n",
    "        tolerance=openmm_unit.Quantity(\n",
    "            value=50.0, unit=openmm_unit.kilojoule_per_mole / openmm_unit.nanometer\n",
    "        )\n",
    "    )\n",
    "    print(f\"Simulating system {idx+1}/{n_systems}\")\n",
    "\n",
    "    traj_filename = f\"trajectory_{idx+1}.dcd\"\n",
    "    dcd_reporter = openmm.app.DCDReporter(traj_filename, 50)\n",
    "    simulation.reporters.append(dcd_reporter)\n",
    "    simulation.context.setVelocitiesToTemperature(300 * openmm_unit.kelvin)\n",
    "    simulation.runForClockTime(15 * openmm_unit.second)\n",
    "    idx_to_sim[idx + 1] = {\n",
    "        \"topology\": this_top,\n",
    "        \"pdb_filename\": pdb_filename,\n",
    "        \"system\": sys,\n",
    "        \"simulation\": simulation,\n",
    "        \"trajectory\": traj_filename,\n",
    "    }"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "188a7ec0",
   "metadata": {},
   "source": [
    "# NGLView"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f320bae7",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import mdtraj\n",
    "import nglview\n",
    "\n",
    "SIM_TO_VISUALIZE = 1\n",
    "\n",
    "trajectory: mdtraj.Trajectory = mdtraj.load(\n",
    "    idx_to_sim[SIM_TO_VISUALIZE][\"trajectory\"],\n",
    "    top=mdtraj.load(idx_to_sim[SIM_TO_VISUALIZE][\"pdb_filename\"]).topology,\n",
    ")\n",
    "\n",
    "view = nglview.show_mdtraj(trajectory.image_molecules())\n",
    "view.add_representation(\"line\", selection=\"water\")\n",
    "view.add_representation(\n",
    "    \"hyperball\", radiusSize=1, radiusScale=0.5, selection=\"not protein and not water\"\n",
    ")\n",
    "view"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}