Molecule_Breakdown_Model_Substituents.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4f012aac",
   "metadata": {},
   "source": [
    "# **_Substituent Database Building_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85a8d437",
   "metadata": {},
   "source": [
    "## Environment Setting Up and Packages Importing:\n",
    "        !conda activate fragment-env \n",
    "        Please create this environment from the .yml file: fragment.yml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4de9af4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import time\n",
    "\n",
    "import pandas as pd\n",
    "from tqdm import tqdm \n",
    "import pickle\n",
    "import numpy as np\n",
    "from collections import Counter\n",
    "\n",
    "#RDKit:\n",
    "from rdkit.Chem import AllChem\n",
    "from rdkit.Chem import rdchem\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem.rdmolops import *\n",
    "from rdkit.Chem import PandasTools\n",
    "from rdkit.Chem.Draw import MolsToGridImage"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dabb6bf",
   "metadata": {},
   "source": [
    "### 1. Input the Data (usually for multiple-input-file cases):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b2c8f6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To check the time running the code:\n",
    "start_time = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb5d0803",
   "metadata": {},
   "source": [
    "#### (1) If you are using command line, for example .csv files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb9302b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#scr_path = sys.argv[1]\n",
    "#filename = (scr_path.split('/')[-1]).split('.')[0]\n",
    "\n",
    "#df_input = pd.read_csv(scr_path,  sep = ' ', header = None, names = ['SMILES'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06f96544",
   "metadata": {},
   "source": [
    "#### (2) If you are using Jupyter Notebook, for example .csv files:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2f32d8c",
   "metadata": {},
   "source": [
    "        Some formats that you may encounter:\n",
    "\n",
    "        (1) Read .sdf file:\n",
    "        f = 'filename.sdf'\n",
    "        df_input = PandasTools.LoadSDF(f, embedProps=True, molColName=None, smilesName='smiles')\n",
    "\n",
    "        (2) Read .csv or csv.gz file:\n",
    "        df_input = pd.read_csv('filename.csv')\n",
    "        df_input = pd.read_csv('filename.csv.gz')\n",
    "\n",
    "        (3) Read .pickle file:\n",
    "        df_input = pd.read_pickle('filename.pickle')\n",
    "\n",
    "        (4) Read .smi files:\n",
    "        df = pd.read_csv(\"filename.smi\", sep=' ', header=None, names=['SMILES'])\n",
    "\n",
    "        (5) Read .xml files:\n",
    "        import xml.etree.ElementTree as ET\n",
    "        xml_data = open('data/drugbank_full_database.xml', 'r').read()\n",
    "        root = ET.XML(xml_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e78c155c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_input = pd.read_csv('example.smi',  sep=' ', header=None, names=['SMILES'])\n",
    "df_input "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96abb6aa",
   "metadata": {},
   "source": [
    "#### For Single SMILES cases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6e735f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For example: \n",
    "# df_input = pd.DataFrame({'SMILES':['CCCC1=CC(=CC=C1B(O)O)C(O)=O']})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dc2860d",
   "metadata": {},
   "source": [
    "#### (3) Clean the Raw Database:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d42c4625",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the duplicates in the original dataset:\n",
    "df_input.drop_duplicates(subset = \"SMILES\", keep = 'last', inplace = True)\n",
    "df_input = df_input.reset_index(drop = True)\n",
    "df_input.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3038c089",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To check the size of the database:\n",
    "len(df_input)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e283d42a",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('This databse contains', len(df_input), 'compounds')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c96c5b2",
   "metadata": {},
   "source": [
    "### 2. Create a Data Output Dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b35dc7c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_output = pd.DataFrame({'Input_Smiles':df_input['SMILES']}) \n",
    "df_output['Fragment_Distribution'] = 'fragment_distribution'\n",
    "df_output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8508728",
   "metadata": {},
   "source": [
    "### 3. Break Molecules from All SMILES:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00b7af01",
   "metadata": {},
   "source": [
    "#### (1) Remove All Bonds Between a Cyclic Atom and an Acyclic Atom:\n",
    "        To obtain all the substituents on rings\n",
    "#### (2) Remove All Bonds Between Any Two Cyclic Atoms:\n",
    "        To break the ring structures\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31829937",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Neutralizing Molecules:\n",
    "def neutralize_atoms(mol):\n",
    "    pattern = Chem.MolFromSmarts(\"[+1!h0!$([*]~[-1,-2,-3,-4]),-1!$([*]~[+1,+2,+3,+4])]\")\n",
    "    at_matches = mol.GetSubstructMatches(pattern)\n",
    "    at_matches_list = [y[0] for y in at_matches]\n",
    "    if len(at_matches_list) > 0:\n",
    "        for at_idx in at_matches_list:\n",
    "            atom = mol.GetAtomWithIdx(at_idx)\n",
    "            chg = atom.GetFormalCharge()\n",
    "            hcount = atom.GetTotalNumHs()\n",
    "            atom.SetFormalCharge(0)\n",
    "            atom.SetNumExplicitHs(hcount - chg)\n",
    "            atom.UpdatePropertyCache()\n",
    "    return mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2a8a0bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "count = 0\n",
    "my_list_invaild = []\n",
    "\n",
    "for q in tqdm(range(0, len(df_output)), desc = 'Loop 1'):\n",
    "\n",
    "    try:  \n",
    "        # 1.remove all the non-bond ions:\n",
    "        qsmiles = df_output['Input_Smiles'][q]\n",
    "        qsmiles = qsmiles.replace('.[Li+]', '')\n",
    "        qsmiles = qsmiles.replace('.[Na+]', '')\n",
    "        qsmiles = qsmiles.replace('.[K+]', '')\n",
    "        qsmiles = qsmiles.replace('.[F-]', '')\n",
    "        qsmiles = qsmiles.replace('.[Cl-]', '')\n",
    "        qsmiles = qsmiles.replace('.[Br-]', '')\n",
    "        qsmiles = qsmiles.replace('.[I-]', '')\n",
    "\n",
    "        qmol = Chem.MolFromSmiles(qsmiles)\n",
    "        neutralize_atoms(qmol)\n",
    "        qmolw = Chem.RWMol(qmol)\n",
    " \n",
    "        # 2.make sure there is no more dissociation \n",
    "        if '.' not in qsmiles:\n",
    "\n",
    "            if qmol and qmolw:\n",
    "\n",
    "                # 3. Save all the cyclic atom in a list\n",
    "                cyclic_list=[]\n",
    "                for atom in qmol.GetAtoms():\n",
    "                    if atom.IsInRing() == True:\n",
    "                        i=atom.GetIdx()\n",
    "                        c = qmol.GetAtomWithIdx(i).GetSymbol()\n",
    "                        cyclic_list.append(c)\n",
    "    \n",
    "                # 4.Remove all the bonds between two cyclic atoms in the same ring:     \n",
    "                for bond in qmol.GetBonds():\n",
    "                    \n",
    "                    if bond.IsInRing() == True:\n",
    "                        qmol_begin=bond.GetBeginAtomIdx()\n",
    "                        qmol_end=bond.GetEndAtomIdx()\n",
    "                        bondatom=(qmol_begin, qmol_end)\n",
    "                        bond_a=qmol.GetAtomWithIdx(qmol_begin).GetSymbol()\n",
    "                        bond_b=qmol.GetAtomWithIdx(qmol_end).GetSymbol()\n",
    "                        if qmol.GetAtomWithIdx(qmol_begin).IsInRing() == qmol.GetAtomWithIdx(qmol_end).IsInRing() == True:\n",
    "                            qmolw.RemoveBond(qmol_begin, qmol_end)\n",
    "                                    \n",
    "                    else:\n",
    "                        qmol_begin=bond.GetBeginAtomIdx()\n",
    "                        qmol_end=bond.GetEndAtomIdx()\n",
    "                        bondatom=(qmol_begin, qmol_end)\n",
    "\n",
    "                        # 5.Remove the cyclic atom in the subtituents:\n",
    "                        if  qmol.GetAtomWithIdx(qmol_begin).IsInRing() != qmol.GetAtomWithIdx(qmol_end).IsInRing():\n",
    "                            qmolw.RemoveBond(qmol_begin, qmol_end)\n",
    "\n",
    "                        # 6.Remove the bonds between two cyclic atons in different rings:\n",
    "                        elif qmol.GetAtomWithIdx(qmol_begin).IsInRing() == qmol.GetAtomWithIdx(qmol_end).IsInRing() == True:\n",
    "                            qmolw.RemoveBond(qmol_begin, qmol_end)\n",
    "                            \n",
    "                new_mol=qmolw.GetMol()\n",
    "                new_smiles = Chem.MolToSmiles(new_mol)\n",
    "                frag_simles = new_smiles.split('.')\n",
    "                mylist = frag_simles\n",
    "\n",
    "                mylist=list(map(lambda item: item.replace('[C@@H]','C'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[C@H]','C'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[C@@]','C'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[C@]','C'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[P@@H]','P'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[P@@]','P'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[P@]','P'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[S@@H]','S'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[S@@]','S'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[S@]','S'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[N@]','N'), mylist))\n",
    "        \n",
    "                mylist=list(map(lambda item: item.replace('[nH3]','N'), mylist))\n",
    "\n",
    "                mylist=list(map(lambda item: item.replace('C-C','CC'), mylist))\n",
    "                \n",
    "                mylist=list(map(lambda item: item.replace('c','C'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('n','N'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('s','S'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('o','O'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('p','P'), mylist))\n",
    "\n",
    "                mylist=list(map(lambda item: item.replace('[ZN','[Zn'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[MO','[Mo'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[AS','[As'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[AC','[Ac'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[SN','[Sn'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[IN','[In'), mylist))\n",
    "                mylist=list(map(lambda item: item.replace('[te','[Te'), mylist))\n",
    "\n",
    "                c1 = Counter(mylist)\n",
    "                c2 = Counter(cyclic_list)\n",
    "                mylist_new = c1-c2\n",
    "\n",
    "                #repeat-count:\n",
    "                my_dict = mylist_new\n",
    "\n",
    "                # 7.Count each fragment only once (single-count)\n",
    "                my_dict = {x: 1 for x in my_dict}\n",
    "\n",
    "                key_list=list(my_dict.keys())\n",
    "                count_list=list(my_dict.values())\n",
    "                df_output['Fragment_Distribution'][q]= my_dict\n",
    "\n",
    "        else:\n",
    "            count += 1\n",
    "            my_list_invaild.append(q)\n",
    "            continue\n",
    "\n",
    "    except:\n",
    "        count += 1\n",
    "        my_list_invaild.append(q)\n",
    "        continue\n",
    "\n",
    "print(count, 'invalids')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4503e13a",
   "metadata": {},
   "source": [
    "### 4. Wrangling the Data Obtained :\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f544595",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.Remove the invalid data: \n",
    "new_df = df_output.drop(my_list_invaild)\n",
    "\n",
    "# Reset the index:\n",
    "new_df = new_df.reset_index(drop = True)\n",
    "\n",
    "print('After remove invalid data:', len(new_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "572dc2b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2.Save the pickle file:\n",
    "with open ('example-each-cpd-breakdown-substituent-without-ROMol.pickle','wb') as f:\n",
    "    pickle.dump(new_df,f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9b02421",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3.Sort the data according to the fragment frequency:\n",
    "fragment_distribution = {\n",
    "    k:[d.get(k) for d in new_df.Fragment_Distribution if k in d]\n",
    "    for k in set().union(*new_df.Fragment_Distribution)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74c40502",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_result=pd.DataFrame({'Fragment':fragment_distribution}) \n",
    "k=list(df_result['Fragment'].keys())\n",
    "arr=df_result.Fragment.values\n",
    "l1=list()\n",
    "for ls in arr:\n",
    "    while None in ls:\n",
    "        ls.remove(None)\n",
    "        next\n",
    "    else:\n",
    "        l1.append(ls)\n",
    "\n",
    "l2=list()\n",
    "for i in l1:\n",
    "    s=sum(i)\n",
    "    l2.append(s)\n",
    "\n",
    "df_result['Fragment_SMILES']=k\n",
    "df_result['Frequency']=l2\n",
    "\n",
    "df_result = df_result.sort_values(by=['Frequency'], ascending=False)\n",
    "df_result = df_result.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39dbe5f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_result = df_result.drop(columns='Fragment')\n",
    "df_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0eab3c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The number of substituents obtained is', len(df_result))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83ec09e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.Save the pickle file:\n",
    "# BE CAREFULLY EXCUTING -- WHEN SAVING AS PICKELS\n",
    "with open ('example-substituents-sorted-without-ROMol.pickle','wb') as f:\n",
    "    pickle.dump(df_result,f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a42504a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5.Finish all the code running and print the key outputs:\n",
    "end_time = time.time()\n",
    "duration = round((end_time - start_time)/60)\n",
    "print('This programme took', duration, 'minutes')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "232e5ea3",
   "metadata": {},
   "source": [
    "### 5. Display the Ring Fragments (ony for single-input-file cases, for multiple-input-file cases please use the 'merge' code):\n",
    "        \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41913d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.Add ROMol to each ring fragment:\n",
    "PandasTools.AddMoleculeColumnToFrame(df_result, smilesCol = \"Fragment_SMILES\")\n",
    "df_result\n",
    "\n",
    "# If the ROMols already have been added before, to display them again:\n",
    "# PandasTools.RenderImagesInAllDataFrames(images=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9836e67",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. Display the results (all or top 10000 frequent ring fragments) in a .html file\n",
    "df_result_top = df_result.head(10000)\n",
    "\n",
    "fmolport = open('example-top10000-substituents-with-ROMol.html','w')\n",
    "h = df_result.to_html()\n",
    "fmolport.write(h)\n",
    "fmolport.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('env-fragmentation')",
   "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.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "078c8f19e0ba7590f95e671d05f61230bb22f953ed32440aa8021a171b742402"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}