BLOX_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "useful-indie",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "RDKit ERROR: [20:24:22] SMILES Parse Error: syntax error while parsing: NNC(=O)CNC(=O)\\C=N\\#N\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: Failed parsing SMILES 'NNC(=O)CNC(=O)\\C=N\\#N' for input: 'NNC(=O)CNC(=O)\\C=N\\#N'\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: syntax error while parsing: O=C1NC(=O)\\C(=N/#N)\\C=N1\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: Failed parsing SMILES 'O=C1NC(=O)\\C(=N/#N)\\C=N1' for input: 'O=C1NC(=O)\\C(=N/#N)\\C=N1'\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: syntax error while parsing: NC(=O)CNC(=O)\\C=N\\#N\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: Failed parsing SMILES 'NC(=O)CNC(=O)\\C=N\\#N' for input: 'NC(=O)CNC(=O)\\C=N\\#N'\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: syntax error while parsing: CCCCN(CC(O)C1=C\\C(=N/#N)\\C(=O)C=C1)N=O\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: Failed parsing SMILES 'CCCCN(CC(O)C1=C\\C(=N/#N)\\C(=O)C=C1)N=O' for input: 'CCCCN(CC(O)C1=C\\C(=N/#N)\\C(=O)C=C1)N=O'\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: syntax error while parsing: NC(COC(=O)\\C=N/#N)C(=O)O\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: Failed parsing SMILES 'NC(COC(=O)\\C=N/#N)C(=O)O' for input: 'NC(COC(=O)\\C=N/#N)C(=O)O'\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: syntax error while parsing: CCN(CC(O)C1=CC(=O)\\C(=N\\#N)\\C=C1)N=O\n",
      "RDKit ERROR: [20:24:22] SMILES Parse Error: Failed parsing SMILES 'CCN(CC(O)C1=CC(=O)\\C(=N\\#N)\\C=C1)N=O' for input: 'CCN(CC(O)C1=CC(=O)\\C(=N\\#N)\\C=C1)N=O'\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: syntax error while parsing: NNC(=O)CNC(=O)\\C=N\\#N\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: Failed parsing SMILES 'NNC(=O)CNC(=O)\\C=N\\#N' for input: 'NNC(=O)CNC(=O)\\C=N\\#N'\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: syntax error while parsing: O=C1NC(=O)\\C(=N/#N)\\C=N1\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: Failed parsing SMILES 'O=C1NC(=O)\\C(=N/#N)\\C=N1' for input: 'O=C1NC(=O)\\C(=N/#N)\\C=N1'\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: syntax error while parsing: NC(=O)CNC(=O)\\C=N\\#N\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: Failed parsing SMILES 'NC(=O)CNC(=O)\\C=N\\#N' for input: 'NC(=O)CNC(=O)\\C=N\\#N'\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: syntax error while parsing: CCCCN(CC(O)C1=C\\C(=N/#N)\\C(=O)C=C1)N=O\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: Failed parsing SMILES 'CCCCN(CC(O)C1=C\\C(=N/#N)\\C(=O)C=C1)N=O' for input: 'CCCCN(CC(O)C1=C\\C(=N/#N)\\C(=O)C=C1)N=O'\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: syntax error while parsing: NC(COC(=O)\\C=N/#N)C(=O)O\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: Failed parsing SMILES 'NC(COC(=O)\\C=N/#N)C(=O)O' for input: 'NC(COC(=O)\\C=N/#N)C(=O)O'\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: syntax error while parsing: CCN(CC(O)C1=CC(=O)\\C(=N\\#N)\\C=C1)N=O\n",
      "RDKit ERROR: [20:24:23] SMILES Parse Error: Failed parsing SMILES 'CCN(CC(O)C1=CC(=O)\\C(=N\\#N)\\C=C1)N=O' for input: 'CCN(CC(O)C1=CC(=O)\\C(=N\\#N)\\C=C1)N=O'\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import urllib.request \n",
    "from rdkit import Chem, rdBase\n",
    "from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools\n",
    "from rdkit.ML.Descriptors import MoleculeDescriptors\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "url = 'https://raw.githubusercontent.com/onecoinbuybus/Database_chemoinformatics/master/smiles_cas_N6512.smi'\n",
    "urllib.request.urlretrieve(url, 'ames.txt') \n",
    "df = pd.read_csv('ames.txt',header=None, sep='\\t') \n",
    "df.columns = ['smiles', 'CAS_NO', 'activity']\n",
    "PandasTools.AddMoleculeColumnToFrame(frame=df, smilesCol='smiles')\n",
    "\n",
    "none_list=[]\n",
    "for i in range(df.shape[0]):\n",
    "    if Chem.MolFromSmiles(df['smiles'][i]) is None:\n",
    "        none_list.append(i)\n",
    "        \n",
    "df=df.drop(none_list)\n",
    "mols=[Chem.MolFromSmiles(smile) for smile in df['smiles']]\n",
    "\n",
    "maccskeys = []\n",
    "for m in mols:\n",
    "    maccskey = [x for x in AllChem.GetMACCSKeysFingerprint(m)]\n",
    "    maccskeys.append(maccskey)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ultimate-folder",
   "metadata": {},
   "outputs": [],
   "source": [
    "descriptor_names = ['qed', 'MolLogP']\n",
    "descriptor_calculator = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)\n",
    "descriptors = pd.DataFrame(\n",
    "    [descriptor_calculator.CalcDescriptors(mol) for mol in mols[:110]],\n",
    "    columns=descriptor_names\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "defined-plaintiff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.18458152,  6.3494    ],\n",
       "       [ 0.64928404,  3.9456    ],\n",
       "       [ 0.39401407, -1.3818    ],\n",
       "       [ 0.41063321,  5.17742   ],\n",
       "       [ 0.88511218,  3.2809    ],\n",
       "       [ 0.45584159,  2.21164   ],\n",
       "       [ 0.30569748,  6.0138    ],\n",
       "       [ 0.3683784 ,  0.6731    ],\n",
       "       [ 0.61264255,  1.8412    ],\n",
       "       [ 0.60254391,  4.1766    ]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "properties=np.array(descriptors)\n",
    "properties[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "sporting-poison",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=np.array(maccskeys[:110])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "demonstrated-paste",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "from multiprocessing import Pool\n",
    "import copy as cp\n",
    "import random\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "import sys, csv, time, argparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "silent-philadelphia",
   "metadata": {},
   "outputs": [],
   "source": [
    "# x\n",
    "features_observed = data[:10]\n",
    "features_unchecked = data[10:]\n",
    "\n",
    "#y\n",
    "properties_observed = properties[:10]\n",
    "properties_unchecked = properties[10:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "decimal-dietary",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc = StandardScaler()\n",
    "sc.fit(features_observed)\n",
    "sc_features_observed = sc.transform(features_observed)\n",
    "sc_features_unchecked = sc.transform(features_unchecked)\n",
    "sc_property = StandardScaler() \n",
    "sc_property.fit(properties_observed)\n",
    "sc_properties_observed = sc_property.transform(properties_observed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "diverse-astronomy",
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_model(prediction_model, x_train, y_train):\n",
    "    if prediction_model == 'RF': \n",
    "        params = {'n_estimators':[10, 50, 100]}\n",
    "        gridsearch = GridSearchCV(RandomForestRegressor(), param_grid=params, cv = 3, scoring=\"r2\", n_jobs=parallel, verbose = 1)\n",
    "        gridsearch.fit(x_train,y_train)\n",
    "        model =  RandomForestRegressor(n_estimators = gridsearch.best_params_['n_estimators'])\n",
    "        model.fit(x_train, y_train)\n",
    "        return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "understood-cleaner",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n"
     ]
    }
   ],
   "source": [
    "parallel = 1\n",
    "\n",
    "#根据不同物性建造不同的model\n",
    "model_list = []\n",
    "for d in range(2):\n",
    "    model = build_model('RF', sc_features_observed, properties_observed[:,d])\n",
    "    model_list.append(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "outstanding-genesis",
   "metadata": {},
   "outputs": [],
   "source": [
    "#将备选的东西里面的值预测一遍\n",
    "\n",
    "predicted_properties_list = []\n",
    "for d in range(2):\n",
    "    predicted_properties_list.append(model_list[d].predict(sc_features_unchecked))\n",
    "predicted_properties_list = np.array(predicted_properties_list).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "solar-swing",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stein Novelty Calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "actual-parallel",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hesgau(x, y, sigma):\n",
    "    dim = len(x)\n",
    "    dist = np.sum(np.power(x-y, 2))\n",
    "    return (dim/sigma - dist/sigma**2)*np.exp(-dist/(2*sigma))\n",
    "\n",
    "def stein_novelty(point, data_list, sigma):\n",
    "    n = len(data_list)\n",
    "    score = 0\n",
    "    score = np.sum([hesgau(point, data_list[k,:], sigma) for k in range(n)])\n",
    "    score = score/(n*(n+1)/2)\n",
    "    return -score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "advisory-elder",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc_predicted_properties_list = sc_property.transform(predicted_properties_list) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "copyrighted-router",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.05654122932088566"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stein_novelty(sc_predicted_properties_list[0], sc_properties_observed, sigma=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "restricted-juice",
   "metadata": {},
   "outputs": [],
   "source": [
    "# finding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "foster-boating",
   "metadata": {},
   "outputs": [],
   "source": [
    "def recommend_next(prediction_model, features_observed, features_unchecked, properties_observed):\n",
    "\n",
    "    sc = StandardScaler()\n",
    "    sc.fit(features_observed)\n",
    "    sc_features_observed = sc.transform(features_observed)\n",
    "    sc_features_unchecked = sc.transform(features_unchecked)\n",
    "    sc_property = StandardScaler() \n",
    "    sc_property.fit(properties_observed)\n",
    "    sc_properties_observed = sc_property.transform(properties_observed)\n",
    "    \n",
    "    #根据不同物性建造不同的model\n",
    "    model_list = []\n",
    "    for d in range(2):\n",
    "        model = build_model(prediction_model, sc_features_observed, properties_observed[:,d])\n",
    "        model_list.append(model)\n",
    "\n",
    "    predicted_properties_list = []\n",
    "    for d in range(2):\n",
    "        predicted_properties_list.append(model_list[d].predict(sc_features_unchecked))\n",
    "    predicted_properties_list = np.array(predicted_properties_list).T\n",
    "    \n",
    "    #Calc. Stein Novelty\n",
    "    sc_predicted_properties_list = sc_property.transform(predicted_properties_list) \n",
    "    sn_data = [stein_novelty(point, sc_properties_observed, sigma=1) for point in sc_predicted_properties_list]\n",
    "    \n",
    "    #Select and save next candidate\n",
    "    maximum_index = np.argmax(sn_data)\n",
    "    \n",
    "    return maximum_index, predicted_properties_list[maximum_index], sn_data[maximum_index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "closing-strap",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exploration: 0\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 54 predicted_properties [0.59782408 0.903816  ] Stein novelty -0.03669726271314275\n",
      "Exploration: 1\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 12 predicted_properties [ 0.44592185 -0.031954  ] Stein novelty -0.034621463301366825\n",
      "Exploration: 2\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 11 predicted_properties [0.45004822 0.6821412 ] Stein novelty -0.03998524519295293\n",
      "Exploration: 3\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 10 predicted_properties [0.41952621 5.3883508 ] Stein novelty -0.0369188367749818\n",
      "Exploration: 4\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 28 predicted_properties [0.40760474 0.841738  ] Stein novelty -0.03912689300512455\n",
      "Exploration: 5\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 2 predicted_properties [0.38500347 2.4778688 ] Stein novelty -0.03932647186115183\n",
      "Exploration: 6\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 69 predicted_properties [0.38655677 1.5327426 ] Stein novelty -0.037895522058455734\n",
      "Exploration: 7\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 34 predicted_properties [0.61650409 1.2718916 ] Stein novelty -0.02718343501167143\n",
      "Exploration: 8\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 31 predicted_properties [0.61024363 1.013306  ] Stein novelty -0.023085220010071485\n",
      "Exploration: 9\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Fitting 3 folds for each of 3 candidates, totalling 9 fits\n",
      "Recommended_index 21 predicted_properties [0.35536024 2.8514252 ] Stein novelty -0.015680248493443138\n"
     ]
    }
   ],
   "source": [
    "num_loop=10\n",
    "\n",
    "for l in range(num_loop):\n",
    "        print('Exploration:', l)\n",
    "        recommended_index, predicted_properties, SN = recommend_next('RF', features_observed, features_unchecked, properties_observed)\n",
    "        print('Recommended_index', recommended_index, 'predicted_properties', predicted_properties, 'Stein novelty', SN)\n",
    "\n",
    "        #Add the experimental or simulation result of the recommended data\n",
    "        features_observed = np.append(features_observed, [features_unchecked[recommended_index]], axis = 0)\n",
    "        properties_observed = np.append(properties_observed, [properties_unchecked[recommended_index]], axis = 0)         \n",
    "\n",
    "        #Removed the recommend data\n",
    "        features_unchecked = np.delete(features_unchecked, recommended_index, axis = 0)\n",
    "        properties_unchecked = np.delete(properties_unchecked, recommended_index, axis = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cutting-repeat",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}