AMES_interpretability.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import os\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"]= '2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "exec(open('AMES_aux.py').read()) #load functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#parameters for data generation\n",
    "pad_len = 1200 #maximal length of data\n",
    "asym = ['Li', 'Na', 'K', 'Be', 'Mg', 'Ca', 'Sc', 'Cr', 'Fe', 'Pt', 'Cu', 'Hg', 'B', 'Al', 'C0', 'C1', 'C2', 'C3', \n",
    "        'Si', 'Sn', 'N0', 'N1', 'N2', 'P0', 'As', 'Sb', 'O0', 'O1', 'S0', 'S1', 'Se', 'F', 'Cl', 'Br', 'I']\n",
    "arom = ['c', 'n', 'o', 's']\n",
    "other = ['X']\n",
    "one_hot = asym + arom + other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#model parameters\n",
    "nr_atoms = int(pad_len/8)\n",
    "input_shape = (150,len(one_hot)+4)\n",
    "output_shape = 1\n",
    "filters = [1024, 1024, 1024]\n",
    "dense_layers = [512]\n",
    "n = 100\n",
    "batch_size = 64\n",
    "opt = optimizers.Adam(lr=0.001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### load precalculated data or ... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pickle.load(open('/publicdata/graphconv/publication/data/ames_pad1200_2410.p','rb'))\n",
    "xTest, pTest, yTest, test_mols, xVal, pVal, yVal, val_mols, xTrain, pTrain, yTrain, tr_mols = data.values()\n",
    "pTest = pTest.toarray().reshape(-1,pad_len,int(pad_len/8+1))\n",
    "pVal = pVal.toarray().reshape(-1,pad_len,int(pad_len/8+1))\n",
    "pTrain = pTrain.toarray().reshape(-1,pad_len,int(pad_len/8+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ... generate data from molecules (not needed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#data generation from scratch\n",
    "xVal, pVal = Generator(val_mols, one_hot, pad_len)\n",
    "xTest, pTest = Generator(test_mols, one_hot, pad_len)\n",
    "xTrain, pTrain = Generator(tr_mols, one_hot, pad_len)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### load model..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_file = 'ames.h5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = load_model(model_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ... or train from scratch (not needed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = create_model(input_shape, output_shape, filters, dense_layers, opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_1 (InputLayer)            (None, 150, 44)      0                                            \n",
      "__________________________________________________________________________________________________\n",
      "input_2 (InputLayer)            (None, 1200, 151)    0                                            \n",
      "__________________________________________________________________________________________________\n",
      "zero_padding1d_1 (ZeroPadding1D (None, 151, 44)      0           input_1[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "dot_1 (Dot)                     (None, 1200, 44)     0           input_2[0][0]                    \n",
      "                                                                 zero_padding1d_1[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_1 (Conv1D)               (None, 600, 1024)    90112       dot_1[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling1d_1 (AveragePoo (None, 150, 1024)    0           conv1d_1[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "zero_padding1d_2 (ZeroPadding1D (None, 151, 1024)    0           average_pooling1d_1[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "dot_2 (Dot)                     (None, 1200, 1024)   0           input_2[0][0]                    \n",
      "                                                                 zero_padding1d_2[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_2 (Conv1D)               (None, 600, 1024)    2097152     dot_2[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling1d_2 (AveragePoo (None, 150, 1024)    0           conv1d_2[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "zero_padding1d_3 (ZeroPadding1D (None, 151, 1024)    0           average_pooling1d_2[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "dot_3 (Dot)                     (None, 1200, 1024)   0           input_2[0][0]                    \n",
      "                                                                 zero_padding1d_3[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_3 (Conv1D)               (None, 600, 1024)    2097152     dot_3[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling1d_3 (AveragePoo (None, 150, 1024)    0           conv1d_3[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling1d_4 (AveragePoo (None, 1, 1024)      0           average_pooling1d_3[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "reshape_1 (Reshape)             (None, 1024)         0           average_pooling1d_4[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "dense_1 (Dense)                 (None, 512)          524800      reshape_1[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "dropout_1 (Dropout)             (None, 512)          0           dense_1[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "dense_2 (Dense)                 (None, 1)            513         dropout_1[0][0]                  \n",
      "==================================================================================================\n",
      "Total params: 4,809,729\n",
      "Trainable params: 4,809,729\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model, tr, val, tst = run_model(model, [xTrain,pTrain], yTrain, [xVal, pVal], yVal, [xTest, pTest], yTest, \n",
    "                                n=n, batch_size=batch_size, model_file='ames_from_scratch')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### evaluate model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pval = model.predict([xVal, pVal])\n",
    "ptst = model.predict([xTest, pTest])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.9188677817301143, 0.8006412380518888)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roc_auc_score(yVal, pval), roc_auc_score(yTest, ptst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### interpret model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = xVal\n",
    "p = pVal\n",
    "y = yVal\n",
    "mols = val_mols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "#layers are a list of indices of the atom representation layers\n",
    "#returns an array of atom representations and molecule representations\n",
    "atom_act, mol_act = get_activations(model, x, p, layers=[5,9,13,15], filters=filters, bs=8, pad_len=pad_len) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "#returns unqiue atom representations, the indices of the molecules from which the atoms originate and the atom indicies\n",
    "#in these molecules\n",
    "unique, idx, a_idx = reduce_act(atom_act, mols, nr_atoms=nr_atoms) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [],
   "source": [
    "#function for atom/substructure scores\n",
    "atom_scoring = K.function([model.layers[16].input],[model.layers[-1].output]) #adjust layers index to your architecture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [],
   "source": [
    "#returns the scores for each atom/substructure 3rd layer representation \n",
    "a_scores = atom_scoring([unique[:,-filters[-1]:]/nr_atoms, False])[0][:,0] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms_sorted = np.argsort(a_scores)[::-1][:50] #take 50 most important substructures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "substr = [get_substruct(np.array(mols)[idx][e],int(a_idx[e]),3) for e in atoms_sorted] #we used the 3 convolutional layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAAEsCAIAAAD2HxkiAAAEhklEQVR4nO3c227aUBBAUVz1/3/ZfSBFqImaw8Xex/FajygCR/JmBptkWdf1AnR+1QcAZydCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKcxbJ8/wg/kggnorpzEuFE1lWHZyRCiIlwLobhCYkwsyxf96bDs/ldH8AZXRtb1/o4mIMIdzWYn2F4KtbRPVw3z2W5rOvoAPz8Y/Nn6Vbnc0zCbY1vnp9/5p9HbuNx5j32+kbDQ0S4lS2CuT7bJCf6l7/g9Z1ihsM7EBG+2fKxga3bnYj5iT7/QD4WEb7NNb8N47uTrKa3D3j/f9H8PeJwRPgGe+Z3s+dq+mjwOnyICJ/3d/PcO797m57ug6OPF4nwGXuOvm8b22I1ff0J3eocJ8LH7L95jsy6d62m7x19Ohy0hKvUsSQf/O5e/XIZaOPpDl8ZfS6WvkiERzLS2KNJyC8nwoMZnHWbXpx0wea9RHg84fwx+rbgwszx7P/lNaNvUyI8qn1uiBt9OxDhgQ3eIfzc6rf1Gn17EuGxDa6m4zPT6NufCH+CwT/VH5l+8tufCPkgv4p/b3EWvkQ2LRGeiA7nJEKIifBcDMMJifB0XICZje+OQswkhJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoiJEGIihJgIISZCiIkQYiKEmAghJkKIiRBiIoSYCCEmQoj9AVnU00W3lc2AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<PIL.Image.Image image mode=RGB size=300x300 at 0x7F9E58E65828>"
      ]
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Draw.MolToImage(Chem.MolFromSmarts(Chem.MolToSmarts(substr[0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}