rdkit/Clustering_conformers_wth_RdKit.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "NMaTRukbUzH2",
        "outputId": "7f577a5a-4378-42c9-fc9d-a92146aef613"
      },
      "outputs": [],
      "source": [
        "# If RdKit is not install uncomment the command below to install\n",
        "!pip install rdkit-pypi"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "IND5U2zqVHfo",
        "outputId": "439a8fc3-8803-43d9-dba2-1f7b430a6c43"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "2022.09.5\n"
          ]
        }
      ],
      "source": [
        "# Import Necessary Libraries\n",
        "import rdkit\n",
        "from rdkit import Chem\n",
        "from rdkit.Chem import Draw\n",
        "from rdkit.Chem import rdDistGeom\n",
        "from rdkit.Chem import rdMolAlign\n",
        "from rdkit.Chem.Draw import IPythonConsole\n",
        "from rdkit.ML.Cluster import Butina\n",
        "print(rdkit.__version__)\n",
        "\n",
        "\"\"\"\n",
        "We will use the newer version of Distance Geomentry (DG) module called ETKDG.\n",
        "The conformer generator ETKDG is a stochastic search method that utilizes distance \n",
        "geometry together with knowledge derived from experimental crystal structures.\n",
        "It has been shown to generate good conformers for acyclic, flexible molecules.\n",
        "This work builds on ETKDG to improve conformer generation of molecules containing\n",
        "small or large aliphatic (i.e., non-aromatic) rings.\n",
        "Reference:\n",
        "Improving Conformer Generation for Small Rings and Macrocycles Based on Distance Geometry and Experimental Torsional-Angle Preferences\n",
        "https://pubs.acs.org/doi/10.1021/acs.jcim.0c00025\n",
        "\n",
        "\"\"\"\n",
        "\n",
        "def generate_conformers(molecule, num_conformers=300, random_seed=0xd06f00d, num_threads=10):\n",
        "    \"\"\"\n",
        "    Generate conformers for a given molecule using the ETKDG method.\n",
        "    \n",
        "    Args:\n",
        "    - molecule (rdkit.Chem.rdchem.Mol): RDKit molecule object.\n",
        "    - num_conformers (int): Number of conformers to generate.\n",
        "    - random_seed (int): Seed for random number generation.\n",
        "    - num_threads (int): Number of threads to use.\n",
        "    \n",
        "    Returns:\n",
        "    - list: List of conformer IDs.\n",
        "    \"\"\"\n",
        "    ps = rdDistGeom.ETKDGv3()\n",
        "    ps.randomSeed = random_seed\n",
        "    ps.numThreads = num_threads\n",
        "    cids = rdDistGeom.EmbedMultipleConfs(molecule, num_conformers, ps)\n",
        "    return cids\n",
        "\n",
        "def align_and_cluster_conformers(molecule, cids, core_smiles, threshold=1.5):\n",
        "    \"\"\"\n",
        "    Align and cluster the conformers of a molecule.\n",
        "    \n",
        "    Args:\n",
        "    - molecule (rdkit.Chem.rdchem.Mol): RDKit molecule object with conformers.\n",
        "    - cids (list): List of conformer IDs.\n",
        "    - core_smiles (str): SMILES string for the core substructure used for alignment.\n",
        "    - threshold (float): RMSD threshold for clustering.\n",
        "    \n",
        "    Returns:\n",
        "    - list: List of clusters.\n",
        "    \"\"\"\n",
        "    dists = []\n",
        "    for i in range(len(cids)):\n",
        "        for j in range(i):\n",
        "            dists.append(rdMolAlign.GetBestRMS(molecule, molecule, i, j))\n",
        "\n",
        "    core = Chem.MolFromSmiles(core_smiles)\n",
        "    rdMolAlign.AlignMolConformers(molecule, atomIds=molecule.GetSubstructMatch(core))\n",
        "    clusters = Butina.ClusterData(dists, len(cids), threshold, isDistData=True, reordering=True)\n",
        "    \n",
        "    return clusters\n",
        "\n",
        "# We are using \"-Methoxycoumarin-4-acetic Acid N-Succinimidyl Ester\" as an example\n",
        "m_noH = Chem.MolFromSmiles('COC1=CC2=C(C=C1)C(=CC(=O)O2)CC(=O)ON3C(=O)CCC3=O')\n",
        "m_H = Chem.AddHs(m_noH)\n",
        "\n",
        "# Generate conformers\n",
        "cids = generate_conformers(m_H)\n",
        "\n",
        "# Remove hydrogens for alignment and clustering\n",
        "m_3d = Chem.RemoveHs(m_H)\n",
        "\n",
        "# Align and cluster conformers\n",
        "core_smiles = 'C1=CC2=C(C=C1)C(=CC(=O)O2)C'\n",
        "clusters = align_and_cluster_conformers(m_3d, cids, core_smiles)\n",
        "\n",
        "# Visualization and further analysis can be added as needed\n"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.9.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}