Lab_3_Multi-modal and Multi-task.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab Three Multi-modal and Multi-task\n",
    "student name: bochuan wu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# requirement\n",
    "In this lab, you will implement a multi-task network (possibly multi-modal) that classifies interactions between compounds and ligands in the ChEMBL database (https://www.ebi.ac.uk/chembl/ ). The objective is to classify which ligands bind to which targets. Each target will be a separate task working from a  shared ligand representation. \n",
    "\n",
    "The following grading rubric will be used:\n",
    "\n",
    "[10 points] Download the ChEMBL database or setup queries to download a subset go the database.(done)\n",
    "\n",
    "[20 points] Filter the database to the top 100 targets in the database. You will need a definition of \"top\" such as the targets with the most assays. From these top 100 targets, save the ligands that have an assay result for each of the targets. (done)\n",
    "\n",
    "[5 points] Binarize the binding affinity for each ligand in the assay. That is, convert the continuous measure of binding to binary. You should use the column 'IC50' for this calculation. Anything below 300 nM should be considered as an active binding. Anything above 10 uM should be considered non-binding (inactive). (done)\n",
    "\n",
    "[10 points] Featurize each ligand using RDKit (https://www.rdkit.org). This will convert the ligand representation into a binary vector of features. Mention any hyper parameters you use.(done)\n",
    "\n",
    "[10 points] Train a multi-task model (one model for each target). Use an 80/20 split for each target.(done)\n",
    "\n",
    "[20 points] Report the results using AUC, BEDROC, and Enrichment factor. These metrics are easily calculated using the RDKit scoring library. Discuss the results (you will need to look up each evaluation metric to interpret the result). You may be interested in page 6 of the following document: https://www.dropbox.com/s/6je37ml475vg3ep/Srinivas2018ImplictDescriptorUnderReview.pdf?dl=0 (done)\n",
    "\n",
    "[10 points] Finally, you have free reign to perform any other analysis. A suggested analysis is to add an additional mode of input data such as another fingerprint for each ligand. (done)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ChemBL\n",
    "ChEMBL is an actively maintained and used chemical database of bioactive molecules with drug-like properties. The European Bioinformatics Institute (EBI) is in charge of keeping up to date with it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratory Data Analysis(EDA)\n",
    "*  Data Collection\n",
    "*  Visualization\n",
    "*  Data Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Data Collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "f = open('data.csv', 'r', encoding=\"ISO-8859-1\")\n",
    "lines = f.readlines()\n",
    "f.close()\n",
    "flis = [] # will have [id, ic50_val, unit, target, smiles] \n",
    "l=0\n",
    "\n",
    "for line in lines[1:]: # first line contains the headers (therefore skiped)\n",
    "    # clean up the dataset \n",
    "    line = line.replace('\\x00','')\n",
    "    line = line.replace('\",','\"!')\n",
    "    line = line.replace('\"','')\n",
    "    lis = line.strip().split('!')\n",
    "    l=1\n",
    "    if len(lis) < 36:\n",
    "        continue\n",
    "    flis.append([lis[0], lis[3], lis[4], lis[8], lis[36]])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1065319"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(flis)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some of the data have smiles or IC50 missing\n",
    "# removed them \n",
    "lists = []\n",
    "for i in flis:\n",
    "    if len(i[1]) != 0 and len(i[-1]) != 0:\n",
    "        #print (i)\n",
    "        lists.append(i)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import DataFrame\n",
    "df = DataFrame(flis, columns = ['ids', 'ic50_val', 'unit', 'target', 'smiles'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 1065319 entries, 0 to 1065318\n",
      "Data columns (total 5 columns):\n",
      "ids         1065319 non-null object\n",
      "ic50_val    1065319 non-null object\n",
      "unit        1065319 non-null object\n",
      "target      1065319 non-null object\n",
      "smiles      1065319 non-null object\n",
      "dtypes: object(5)\n",
      "memory usage: 40.6+ MB\n"
     ]
    }
   ],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ids</th>\n",
       "      <th>ic50_val</th>\n",
       "      <th>unit</th>\n",
       "      <th>target</th>\n",
       "      <th>smiles</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CHEMBL17430</td>\n",
       "      <td>11600</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL664913</td>\n",
       "      <td>Nc1nc(N)c2c(CSc3ccc(Cl)c(Cl)c3)c[nH]c2n1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CHEMBL119370</td>\n",
       "      <td>9000</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL812106</td>\n",
       "      <td>NC(=S)\\N=N\\c1c(O)[nH]c2ccc(I)cc12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CHEMBL239370</td>\n",
       "      <td>17900</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL811877</td>\n",
       "      <td>[Cl-].C[N+](C)(CCCN1c2ccccc2Sc3ccc(Cl)cc13)Cc4...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CHEMBL61415</td>\n",
       "      <td>4200</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL664920</td>\n",
       "      <td>COc1cc(OC)cc(c1)N(C)Cc2c[nH]c3nc(N)nc(O)c23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CHEMBL22</td>\n",
       "      <td>2700</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL664913</td>\n",
       "      <td>COc1cc(Cc2cnc(N)nc2N)cc(OC)c1OC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            ids ic50_val unit        target  \\\n",
       "0   CHEMBL17430    11600   nM  CHEMBL664913   \n",
       "1  CHEMBL119370     9000   nM  CHEMBL812106   \n",
       "2  CHEMBL239370    17900   nM  CHEMBL811877   \n",
       "3   CHEMBL61415     4200   nM  CHEMBL664920   \n",
       "4      CHEMBL22     2700   nM  CHEMBL664913   \n",
       "\n",
       "                                              smiles  \n",
       "0           Nc1nc(N)c2c(CSc3ccc(Cl)c(Cl)c3)c[nH]c2n1  \n",
       "1                  NC(=S)\\N=N\\c1c(O)[nH]c2ccc(I)cc12  \n",
       "2  [Cl-].C[N+](C)(CCCN1c2ccccc2Sc3ccc(Cl)cc13)Cc4...  \n",
       "3        COc1cc(OC)cc(c1)N(C)Cc2c[nH]c3nc(N)nc(O)c23  \n",
       "4                    COc1cc(Cc2cnc(N)nc2N)cc(OC)c1OC  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#count frequency of target\n",
    "def count(df):\n",
    "    count_dic=dict(zip(*np.unique(df['target'],return_counts=True)))\n",
    "    return df['target'].apply(lambda x:count_dic[x] ) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['count'] = count(df) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ids</th>\n",
       "      <th>ic50_val</th>\n",
       "      <th>unit</th>\n",
       "      <th>target</th>\n",
       "      <th>smiles</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CHEMBL17430</td>\n",
       "      <td>11600</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL664913</td>\n",
       "      <td>Nc1nc(N)c2c(CSc3ccc(Cl)c(Cl)c3)c[nH]c2n1</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CHEMBL119370</td>\n",
       "      <td>9000</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL812106</td>\n",
       "      <td>NC(=S)\\N=N\\c1c(O)[nH]c2ccc(I)cc12</td>\n",
       "      <td>33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CHEMBL239370</td>\n",
       "      <td>17900</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL811877</td>\n",
       "      <td>[Cl-].C[N+](C)(CCCN1c2ccccc2Sc3ccc(Cl)cc13)Cc4...</td>\n",
       "      <td>19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CHEMBL61415</td>\n",
       "      <td>4200</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL664920</td>\n",
       "      <td>COc1cc(OC)cc(c1)N(C)Cc2c[nH]c3nc(N)nc(O)c23</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CHEMBL22</td>\n",
       "      <td>2700</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL664913</td>\n",
       "      <td>COc1cc(Cc2cnc(N)nc2N)cc(OC)c1OC</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            ids ic50_val unit        target  \\\n",
       "0   CHEMBL17430    11600   nM  CHEMBL664913   \n",
       "1  CHEMBL119370     9000   nM  CHEMBL812106   \n",
       "2  CHEMBL239370    17900   nM  CHEMBL811877   \n",
       "3   CHEMBL61415     4200   nM  CHEMBL664920   \n",
       "4      CHEMBL22     2700   nM  CHEMBL664913   \n",
       "\n",
       "                                              smiles  count  \n",
       "0           Nc1nc(N)c2c(CSc3ccc(Cl)c(Cl)c3)c[nH]c2n1     16  \n",
       "1                  NC(=S)\\N=N\\c1c(O)[nH]c2ccc(I)cc12     33  \n",
       "2  [Cl-].C[N+](C)(CCCN1c2ccccc2Sc3ccc(Cl)cc13)Cc4...     19  \n",
       "3        COc1cc(OC)cc(c1)N(C)Cc2c[nH]c3nc(N)nc(O)c23     17  \n",
       "4                    COc1cc(Cc2cnc(N)nc2N)cc(OC)c1OC     16  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "871\n"
     ]
    }
   ],
   "source": [
    "#find the frequency count for 100th target\n",
    "rank = df['target'].value_counts()\n",
    "c=0\n",
    "for i in rank:\n",
    "    if c == 100:\n",
    "        print(i)\n",
    "        break\n",
    "    c += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[df['count']>=871]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert ic50 to float type\n",
    "import pandas as pd\n",
    "\n",
    "df['ic50_val']= pd.to_numeric(df['ic50_val'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 111077 entries, 16863 to 1064540\n",
      "Data columns (total 6 columns):\n",
      "ids         111077 non-null object\n",
      "ic50_val    24776 non-null float64\n",
      "unit        111077 non-null object\n",
      "target      111077 non-null object\n",
      "smiles      111077 non-null object\n",
      "count       111077 non-null int64\n",
      "dtypes: float64(1), int64(1), object(4)\n",
      "memory usage: 5.9+ MB\n"
     ]
    }
   ],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter(x):\n",
    "    if x > 10000:\n",
    "        return int(1)\n",
    "    elif x < 300:\n",
    "        return int(0)\n",
    "    else:\n",
    "        return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['output']=df['ic50_val'].apply(filter)\n",
    "df = df.dropna(axis=0,how='any')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 15259 entries, 17840 to 1064540\n",
      "Data columns (total 7 columns):\n",
      "ids         15259 non-null object\n",
      "ic50_val    15259 non-null float64\n",
      "unit        15259 non-null object\n",
      "target      15259 non-null object\n",
      "smiles      15259 non-null object\n",
      "count       15259 non-null int64\n",
      "output      15259 non-null float64\n",
      "dtypes: float64(2), int64(1), object(4)\n",
      "memory usage: 953.7+ KB\n"
     ]
    }
   ],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ids</th>\n",
       "      <th>ic50_val</th>\n",
       "      <th>unit</th>\n",
       "      <th>target</th>\n",
       "      <th>smiles</th>\n",
       "      <th>count</th>\n",
       "      <th>output</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>17840</th>\n",
       "      <td>CHEMBL1479490</td>\n",
       "      <td>50000.00</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL1614432</td>\n",
       "      <td>Cc1ccc(cc1)N(CC(=O)NC2CCCC2)C(=O)c3ccc(C)s3</td>\n",
       "      <td>1128</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17843</th>\n",
       "      <td>CHEMBL1502002</td>\n",
       "      <td>18963.56</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL1614432</td>\n",
       "      <td>CCOC(=O)C1CCCN(C1)c2ccc(c3nonc23)[N+](=O)[O-]</td>\n",
       "      <td>1128</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17845</th>\n",
       "      <td>CHEMBL1462559</td>\n",
       "      <td>14819.11</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL1614432</td>\n",
       "      <td>OCc1ccc(CO[C@H]2C[C@H](C=C(O2)C(=O)N3CCCCCCC3)...</td>\n",
       "      <td>1128</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18068</th>\n",
       "      <td>CHEMBL1313596</td>\n",
       "      <td>13997.64</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL1614432</td>\n",
       "      <td>COc1ccc(COC(=O)Nc2cc3c(nn(C)c3s2)C(F)(F)F)cc1</td>\n",
       "      <td>1128</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18072</th>\n",
       "      <td>CHEMBL3191794</td>\n",
       "      <td>15131.44</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL1614432</td>\n",
       "      <td>Cc1[nH]c2ccccc2c1\\C=N\\NC(=S)NCc3ccccc3</td>\n",
       "      <td>1128</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 ids  ic50_val unit         target  \\\n",
       "17840  CHEMBL1479490  50000.00   nM  CHEMBL1614432   \n",
       "17843  CHEMBL1502002  18963.56   nM  CHEMBL1614432   \n",
       "17845  CHEMBL1462559  14819.11   nM  CHEMBL1614432   \n",
       "18068  CHEMBL1313596  13997.64   nM  CHEMBL1614432   \n",
       "18072  CHEMBL3191794  15131.44   nM  CHEMBL1614432   \n",
       "\n",
       "                                                  smiles  count  output  \n",
       "17840        Cc1ccc(cc1)N(CC(=O)NC2CCCC2)C(=O)c3ccc(C)s3   1128     1.0  \n",
       "17843      CCOC(=O)C1CCCN(C1)c2ccc(c3nonc23)[N+](=O)[O-]   1128     1.0  \n",
       "17845  OCc1ccc(CO[C@H]2C[C@H](C=C(O2)C(=O)N3CCCCCCC3)...   1128     1.0  \n",
       "18068      COc1ccc(COC(=O)Nc2cc3c(nn(C)c3s2)C(F)(F)F)cc1   1128     1.0  \n",
       "18072             Cc1[nH]c2ccccc2c1\\C=N\\NC(=S)NCc3ccccc3   1128     1.0  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2b5033bb1d0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAECCAYAAAAVYxsVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAELhJREFUeJzt3X+s3XV9x/Hna+1w/tgE5EqwrSubnQ7NnKzDbibLIgsUNZY/JKlZRuOadNlw023JhO2PJiqJZMtwZKJrpLMYQyXMhU5R1qDGLJPCRQgKiL0BpVeYXNPC3Ixi8b0/7qfj0M/pr3uuPRfP85HcnO/3/fl8v/d9kravfn+dk6pCkqRBPzPuBiRJS4/hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpM7ycTewUGeccUatXr163G1I0nPKXXfd9d2qmjrWvOdsOKxevZrp6elxtyFJzylJvnU88zytJEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqPGefkH6uWH35Z8bdwk+Nb37gzeNuQZoYHjlIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjrHDIck25M8nuRrA7W/TfL1JPcm+dckpw6MXZFkJsmDSS4cqK9vtZkklw/Uz06yJ8neJJ9McspivkFJ0ok7niOHjwHrD6vtBl5TVb8GfAO4AiDJOcBG4NVtm2uTLEuyDPgQcBFwDvD2NhfgKuDqqloDHAA2j/SOJEkjO2Y4VNWXgP2H1f69qg621duBlW15A7Czqn5YVQ8DM8B57Wemqh6qqqeAncCGJAHeCNzUtt8BXDzie5IkjWgxrjn8IfDZtrwC2DcwNttqR6q/BHhiIGgO1YdKsiXJdJLpubm5RWhdkjTMSOGQ5G+Ag8AnDpWGTKsF1Ieqqm1Vtbaq1k5NTZ1ou5Kk47TgD95Lsgl4C3B+VR36B30WWDUwbSXwaFseVv8ucGqS5e3oYXC+JGlMFnTkkGQ98B7grVX1/YGhXcDGJM9LcjawBrgDuBNY0+5MOoX5i9a7Wqh8AXhb234TcPPC3ookabEcz62sNwBfBl6ZZDbJZuAfgZ8Hdie5J8lHAKrqPuBG4H7gc8BlVfV0Oyp4J3Ar8ABwY5sL8yHzF0lmmL8Gcd2ivkNJ0gk75mmlqnr7kPIR/wGvqiuBK4fUbwFuGVJ/iPm7mSRJS4RPSEuSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOscMhyTbkzye5GsDtdOT7E6yt72e1upJck2SmST3Jjl3YJtNbf7eJJsG6r+R5Kttm2uSZLHfpCTpxBzPkcPHgPWH1S4HbquqNcBtbR3gImBN+9kCfBjmwwTYCrweOA/YeihQ2pwtA9sd/rskSSfZMcOhqr4E7D+svAHY0ZZ3ABcP1K+vebcDpyY5C7gQ2F1V+6vqALAbWN/GfqGqvlxVBVw/sC9J0pgs9JrDmVX1GEB7fWmrrwD2DcybbbWj1WeH1CVJY7TYF6SHXS+oBdSH7zzZkmQ6yfTc3NwCW5QkHctCw+E77ZQQ7fXxVp8FVg3MWwk8eoz6yiH1oapqW1Wtraq1U1NTC2xdknQsCw2HXcChO442ATcP1C9tdy2tA55sp51uBS5Iclq7EH0BcGsb+16Sde0upUsH9iVJGpPlx5qQ5Abgd4Ezkswyf9fRB4Abk2wGHgEuadNvAd4EzADfB94BUFX7k7wPuLPNe29VHbrI/cfM3xH1fOCz7UeSNEbHDIeqevsRhs4fMreAy46wn+3A9iH1aeA1x+pDknTy+IS0JKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOiOFQ5I/T3Jfkq8luSHJzyU5O8meJHuTfDLJKW3u89r6TBtfPbCfK1r9wSQXjvaWJEmjWnA4JFkB/BmwtqpeAywDNgJXAVdX1RrgALC5bbIZOFBVrwCubvNIck7b7tXAeuDaJMsW2pckaXSjnlZaDjw/yXLgBcBjwBuBm9r4DuDitryhrdPGz0+SVt9ZVT+sqoeBGeC8EfuSJI1gweFQVd8G/g54hPlQeBK4C3iiqg62abPAira8AtjXtj3Y5r9ksD5kG0nSGIxyWuk05v/XfzbwMuCFwEVDptahTY4wdqT6sN+5Jcl0kum5ubkTb1qSdFxGOa30e8DDVTVXVT8CPgX8NnBqO80EsBJ4tC3PAqsA2viLgf2D9SHbPEtVbauqtVW1dmpqaoTWJUlHM0o4PAKsS/KCdu3gfOB+4AvA29qcTcDNbXlXW6eNf76qqtU3truZzgbWAHeM0JckaUTLjz1luKrak+Qm4CvAQeBuYBvwGWBnkve32nVtk+uAjyeZYf6IYWPbz31JbmQ+WA4Cl1XV0wvtS5I0ugWHA0BVbQW2HlZ+iCF3G1XVD4BLjrCfK4ErR+lFkrR4fEJaktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQZKRySnJrkpiRfT/JAkt9KcnqS3Un2ttfT2twkuSbJTJJ7k5w7sJ9Nbf7eJJtGfVOSpNGMeuTwD8DnqupVwGuBB4DLgduqag1wW1sHuAhY0362AB8GSHI6sBV4PXAesPVQoEiSxmPB4ZDkF4DfAa4DqKqnquoJYAOwo03bAVzcljcA19e824FTk5wFXAjsrqr9VXUA2A2sX2hfkqTRjXLk8EvAHPDPSe5O8tEkLwTOrKrHANrrS9v8FcC+ge1nW+1IdUnSmIwSDsuBc4EPV9XrgP/lmVNIw2RIrY5S73eQbEkynWR6bm7uRPuVJB2nUcJhFpitqj1t/Sbmw+I77XQR7fXxgfmrBrZfCTx6lHqnqrZV1dqqWjs1NTVC65Kko1lwOFTVfwH7kryylc4H7gd2AYfuONoE3NyWdwGXtruW1gFPttNOtwIXJDmtXYi+oNUkSWOyfMTt/xT4RJJTgIeAdzAfODcm2Qw8AlzS5t4CvAmYAb7f5lJV+5O8D7izzXtvVe0fsS9J0ghGCoequgdYO2To/CFzC7jsCPvZDmwfpRdJ0uLxCWlJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1Rg6HJMuS3J3k02397CR7kuxN8skkp7T689r6TBtfPbCPK1r9wSQXjtqTJGk0i3Hk8C7ggYH1q4Crq2oNcADY3OqbgQNV9Qrg6jaPJOcAG4FXA+uBa5MsW4S+JEkLNFI4JFkJvBn4aFsP8EbgpjZlB3BxW97Q1mnj57f5G4CdVfXDqnoYmAHOG6UvSdJoRj1y+CDwV8CP2/pLgCeq6mBbnwVWtOUVwD6ANv5km///9SHbSJLGYMHhkOQtwONVdddgecjUOsbY0bY5/HduSTKdZHpubu6E+pUkHb9RjhzeALw1yTeBncyfTvogcGqS5W3OSuDRtjwLrAJo4y8G9g/Wh2zzLFW1rarWVtXaqampEVqXJB3NgsOhqq6oqpVVtZr5C8qfr6rfB74AvK1N2wTc3JZ3tXXa+Oerqlp9Y7ub6WxgDXDHQvuSJI1u+bGnnLD3ADuTvB+4G7iu1a8DPp5khvkjho0AVXVfkhuB+4GDwGVV9fRPoC9J0nFalHCoqi8CX2zLDzHkbqOq+gFwyRG2vxK4cjF6kSSNziekJUkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1PlJfIe0pOeA1Zd/Ztwt/FT55gfePO4WFpVHDpKkjuEgSeoYDpKkjuEgSeoYDpKkzoLDIcmqJF9I8kCS+5K8q9VPT7I7yd72elqrJ8k1SWaS3Jvk3IF9bWrz9ybZNPrbkiSNYpQjh4PAX1bVrwLrgMuSnANcDtxWVWuA29o6wEXAmvazBfgwzIcJsBV4PXAesPVQoEiSxmPB4VBVj1XVV9ry94AHgBXABmBHm7YDuLgtbwCur3m3A6cmOQu4ENhdVfur6gCwG1i/0L4kSaNblGsOSVYDrwP2AGdW1WMwHyDAS9u0FcC+gc1mW+1I9WG/Z0uS6STTc3Nzi9G6JGmIkcMhyYuAfwHeXVX/fbSpQ2p1lHpfrNpWVWurau3U1NSJNytJOi4jhUOSn2U+GD5RVZ9q5e+000W018dbfRZYNbD5SuDRo9QlSWMyyt1KAa4DHqiqvx8Y2gUcuuNoE3DzQP3SdtfSOuDJdtrpVuCCJKe1C9EXtJokaUxG+eC9NwB/AHw1yT2t9tfAB4Abk2wGHgEuaWO3AG8CZoDvA+8AqKr9Sd4H3Nnmvbeq9o/QlyRpRAsOh6r6D4ZfLwA4f8j8Ai47wr62A9sX2oskaXH5hLQkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqbNkwiHJ+iQPJplJcvm4+5GkSbYkwiHJMuBDwEXAOcDbk5wz3q4kaXItiXAAzgNmquqhqnoK2AlsGHNPkjSxlko4rAD2DazPtpokaQyWj7uBJkNq1U1KtgBb2ur/JHnwJ9rV5DgD+O64mziWXDXuDjQm/vlcXL94PJOWSjjMAqsG1lcCjx4+qaq2AdtOVlOTIsl0Va0ddx/SMP75HI+lclrpTmBNkrOTnAJsBHaNuSdJmlhL4sihqg4meSdwK7AM2F5V9425LUmaWEsiHACq6hbglnH3MaE8VaelzD+fY5Cq7rqvJGnCLZVrDpKkJcRwkCR1DAdJUsdwkLTkJDk9yWnj7mOSGQ4TLMmZSc5N8rokZ467H022JC9PsjPJHLAHuDPJ4622erzdTR7vVppASX4d+AjwYuDbrbwSeAL4k6r6yrh60+RK8mXgg8BNVfV0qy0DLgHeXVXrxtnfpDEcJlCSe4A/qqo9h9XXAf9UVa8dT2eaZEn2VtWaEx3TT8aSeQhOJ9ULDw8GgKq6PckLx9GQBNyV5FpgB898SvMqYBNw99i6mlAeOUygJNcAvwxcz7P/El4KPFxV7xxXb5pc7XPVNjP/XS4rmP+05n3AvwHXVdUPx9jexDEcJlSSi3j2X8JZYFf7GBNJE85wkLTkJXlLVX163H1MEm9l1bO0L1SSlprfHHcDk8YL0jrcsG/lk06KJK/imdOdxfyXfu2qqq1jbWwCeeSgwz017gY0mZK8B9jJ/H9Q7mD+S8AC3JDk8nH2Nom85qBnSfJIVb183H1o8iT5BvDqqvrRYfVTgPt8zuHk8rTSBEpy75GGAD9GQ+PyY+BlwLcOq5/VxnQSGQ6T6UzgQuDAYfUA/3ny25EAeDdwW5K9PPP8zcuBVwA+e3OSGQ6T6dPAi6rqnsMHknzx5LcjQVV9LsmvAOfx7Odv7jz0WUs6ebzmIEnqeLeSJKljOEiSOoaDJKljOEiSOoaDJKnzf+yqdA9zNK7XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "# it is a real imbalance dataset, and it may influence the result\n",
    "df['output'].value_counts().plot.bar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15259"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "targets=[]\n",
    "for t in df['target']:\n",
    "    targets.append(t)\n",
    "len(targets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique tasks  : 111\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of unique tasks  : {}\".format( df['target'].nunique()))\n",
    "tasks = df['target'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using RDkit to get feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import MACCSkeys, AllChem\n",
    "from rdkit.Avalon import pyAvalonTools as fpAvalon\n",
    "from rdkit.Chem.AtomPairs import Pairs, Torsions\n",
    "from rdkit.Chem.Fingerprints import FingerprintMols\n",
    "from rdkit.Chem.ChemicalFeatures import BuildFeatureFactory\n",
    "from rdkit.Chem import rdMolDescriptors\n",
    "import pandas as pd\n",
    "import argparse\n",
    "import textwrap\n",
    "import getopt\n",
    "from argparse import ArgumentParser\n",
    "import time\n",
    "trains=[]\n",
    "c = 0\n",
    "def CalculateFP_mod(smiles, N_BITS=512):\n",
    "    m = Chem.MolFromSmiles(smiles)\n",
    "    if m is None:\n",
    "        raise ValueError('SMILES cannot be converted to a RDKit molecules:', smiles)\n",
    "\n",
    "    return AllChem.GetMorganFingerprintAsBitVect(m, 2, nBits=N_BITS)\n",
    "for smile in df['smiles']:\n",
    "    temp = CalculateFP_mod(smile)\n",
    "    strFP = temp.ToBitString()\n",
    "    trains.append(strFP)\n",
    "\n",
    "\n",
    "\n",
    "          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "data=[]\n",
    "for train in trains:\n",
    "    data.append(list(train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dataset length: 15259\n",
      "feature length:512\n"
     ]
    }
   ],
   "source": [
    "print('dataset length: {}'.format(len(data)))\n",
    "print('feature length:{}'.format(len(data[0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['smiles'] = data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ids</th>\n",
       "      <th>ic50_val</th>\n",
       "      <th>unit</th>\n",
       "      <th>target</th>\n",
       "      <th>smiles</th>\n",
       "      <th>count</th>\n",
       "      <th>output</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>17840</th>\n",
       "      <td>CHEMBL1479490</td>\n",
       "      <td>50000.0</td>\n",
       "      <td>nM</td>\n",
       "      <td>CHEMBL1614432</td>\n",
       "      <td>[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...</td>\n",
       "      <td>1128</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 ids  ic50_val unit         target  \\\n",
       "17840  CHEMBL1479490   50000.0   nM  CHEMBL1614432   \n",
       "\n",
       "                                                  smiles  count  output  \n",
       "17840  [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...   1128     1.0  "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-Task Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_comp_data = []\n",
    "all_comp_test = []\n",
    "from sklearn.model_selection import train_test_split\n",
    "#build dataset\n",
    "%matplotlib inline\n",
    "from sklearn.datasets.base import Bunch\n",
    "for t in range(len(tasks)):\n",
    "    ID = t\n",
    "    y = np.array(df[df.target.str.contains(tasks[t])]['output'])\n",
    "    x = np.array(df[df.target.str.contains(tasks[t])]['smiles'])\n",
    "    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=12)\n",
    "    \n",
    "    for row, target in zip(X_train, y_train):\n",
    "        instance = {}\n",
    "        instance[\"data\"] =np.array(row).reshape((1,512))\n",
    "        instance[\"subject\"] = ID\n",
    "        instance[\"target\"] = [target]\n",
    "        all_comp_data.append(instance)\n",
    "        \n",
    "    for row, target in zip(X_test, y_test):\n",
    "        instance = {}\n",
    "        instance[\"data\"] = np.array(row).reshape(( 1 ,512))\n",
    "        instance[\"subject\"] = ID\n",
    "        instance[\"target\"] = [target]\n",
    "        all_comp_test.append(instance)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From D:\\tool\\lib\\site-packages\\tensorflow\\python\\framework\\op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Colocations handled automatically by placer.\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_1 (InputLayer)         (None, 512)               0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 64)                32832     \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 64)                4160      \n",
      "_________________________________________________________________\n",
      "dense_223 (Dense)            (None, 10)                650       \n",
      "_________________________________________________________________\n",
      "dense_224 (Dense)            (None, 1)                 11        \n",
      "=================================================================\n",
      "Total params: 37,653\n",
      "Trainable params: 37,653\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "from keras.layers import Input, Dense\n",
    "from keras.models import Model\n",
    "\n",
    "inputs = Input(shape=(512,))\n",
    "shared_mlp = Dense(64, activation='relu')(inputs)\n",
    "shared_mlp = Dense(64, activation='relu')(shared_mlp)\n",
    "\n",
    "models_mtl = []\n",
    "\n",
    "for _ in range(len(tasks)):\n",
    "    personal_layer = Dense(10, activation='relu')(shared_mlp)\n",
    "    output_layer = Dense(1, activation='sigmoid')(personal_layer)\n",
    "    model = Model(inputs=inputs, outputs=output_layer)\n",
    "    model.compile(optimizer='rmsprop', loss='binary_crossentropy')\n",
    "    models_mtl.append(model)\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "random.shuffle(all_comp_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From D:\\tool\\lib\\site-packages\\tensorflow\\python\\ops\\math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use tf.cast instead.\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "from keras.backend.tensorflow_backend import set_session\n",
    "config = tf.ConfigProto()\n",
    "config.gpu_options.allocator_type = 'BFC' #A \"Best-fit with coalescing\" algorithm, simplified from a version of dlmalloc.\n",
    "config.gpu_options.per_process_gpu_memory_fraction = 0.3\n",
    "config.gpu_options.allow_growth = True\n",
    "set_session(tf.Session(config=config))\n",
    "\n",
    "EPOCHS = 10\n",
    "i=0\n",
    "for _underscore in range(10):\n",
    "    for instance in all_comp_data:\n",
    "        models_mtl[instance[\"subject\"]].train_on_batch(x = instance[\"data\"], y = instance[\"target\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i , model in enumerate(models_mtl):\n",
    "    model.save('./model/model{}.h5'.format(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prediction and AUC\n",
    "ROC's full name is 'Receiver Operating Characteristic'. and AUC(Area Under the Curve) is the area under that ROC Curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9914954957598522\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import roc_auc_score\n",
    "computer_ya = [ [] for x in range(190)]\n",
    "computer_yhata = [ [] for x in range(190)]\n",
    "# generate prediction for data\n",
    "preds = []\n",
    "y_tests = []\n",
    "for test in all_comp_test:\n",
    "    X_test = test[\"data\"]\n",
    "    y_test = test[\"target\"]\n",
    "    \n",
    "    mlp = models_mtl[test[\"subject\"]]\n",
    "    \n",
    "    yhat_mlp = mlp.predict(X_test)\n",
    "    preds.append(yhat_mlp[0])\n",
    "    y_tests.append(y_test)\n",
    "auc  = roc_auc_score(y_tests, preds)\n",
    "print(auc)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Enrichment Factor\n",
    "Enrichment Factor is defined as the ratio of correctly classified active ligand within a predefined early recognition threshold compared to the total ligands selected by the model, further normalized by the expected random selection of the ligands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#source code from rdkit.ML.scoring.scoring\n",
    "import math\n",
    "from collections import namedtuple\n",
    "\n",
    "def CalcEnrichment(scores, col, fractions):\n",
    "  \"\"\" Determines the enrichment factor for a set of fractions \"\"\"\n",
    "  numMol = len(scores)\n",
    "  if numMol == 0:\n",
    "    raise ValueError('score list is empty')\n",
    "  if len(fractions) == 0:\n",
    "    raise ValueError('fraction list is empty')\n",
    "  for i in fractions:\n",
    "    if i > 1 or i < 0:\n",
    "      raise ValueError('fractions must be between [0,1]')\n",
    "\n",
    "  numPerFrac = [math.ceil(numMol * f) for f in fractions]\n",
    "  numPerFrac.append(numMol)\n",
    "  numActives = 0\n",
    "  enrich = []\n",
    "\n",
    "  # loop over score list\n",
    "  for i in range(numMol):\n",
    "    if i > (numPerFrac[0] - 1) and i > 0:\n",
    "      enrich.append(1.0 * numActives * numMol / i)\n",
    "      numPerFrac.pop(0)\n",
    "    active = scores[i][col]\n",
    "    if active:\n",
    "      numActives += 1\n",
    "\n",
    "  if numActives > 0:  # check that there are actives\n",
    "    enrich = [e / numActives for e in enrich]\n",
    "  else:\n",
    "    enrich = [0.0] * len(fractions)\n",
    "  return enrich"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction = []\n",
    "for pred in preds:\n",
    "    prediction.append(pred[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "tests=[]\n",
    "for y_test in y_tests:\n",
    "    tests.append(y_test[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "scores = list(zip(prediction,tests))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.2060653188180404]"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CalcEnrichment( scores,1,[0.01,])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BEDROC \n",
    "original defined here: Truchon, J. & Bayly, C.I. Evaluating Virtual Screening Methods: Good and Bad Metric for the “Early Recognition” Problem. J. Chem. Inf. Model. 47, 488-508 (2007)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "#source code from rdkit.ML.scoring.scoring\n",
    "def _RIEHelper(scores, col, alpha):\n",
    "  numMol = len(scores)\n",
    "  alpha = float(alpha)\n",
    "  if numMol == 0:\n",
    "    raise ValueError('score list is empty')\n",
    "  if alpha <= 0.0:\n",
    "    raise ValueError('alpha must be greater than zero')\n",
    "\n",
    "  denom = 1.0 / numMol * ((1 - math.exp(-alpha)) / (math.exp(alpha / numMol) - 1))\n",
    "  numActives = 0\n",
    "  sum_exp = 0\n",
    "\n",
    "  # loop over score list\n",
    "  for i in range(numMol):\n",
    "    active = scores[i][col]\n",
    "    if active:\n",
    "      numActives += 1\n",
    "      sum_exp += math.exp(-(alpha * (i + 1)) / numMol)\n",
    "\n",
    "  if numActives > 0:  # check that there are actives\n",
    "    RIE = sum_exp / (numActives * denom)\n",
    "  else:\n",
    "    RIE = 0.0\n",
    "\n",
    "  return RIE, numActives\n",
    "\n",
    "def CalcBEDROC(scores, col, alpha):\n",
    "  \"\"\" BEDROC original defined here:\n",
    "    Truchon, J. & Bayly, C.I.\n",
    "    Evaluating Virtual Screening Methods: Good and Bad Metric for the \"Early Recognition\"\n",
    "    Problem. J. Chem. Inf. Model. 47, 488-508 (2007).\n",
    "    ** Arguments**\n",
    "      - scores: 2d list or numpy array\n",
    "             0th index representing sample\n",
    "             scores must be in sorted order with low indexes \"better\"\n",
    "             scores[sample_id] = vector of sample data\n",
    "      -  col: int\n",
    "             Index of sample data which reflects true label of a sample\n",
    "             scores[sample_id][col] = True iff that sample is active\n",
    "      -  alpha: float\n",
    "             hyper parameter from the initial paper for how much to enrich the top\n",
    "     **Returns**\n",
    "       float BedROC score\n",
    "    \"\"\"\n",
    "  # calculate RIE\n",
    "  RIE, numActives = _RIEHelper(scores, col, alpha)\n",
    "\n",
    "  if numActives > 0:\n",
    "    numMol = len(scores)\n",
    "    ratio = 1.0 * numActives / numMol\n",
    "    RIEmax = (1 - math.exp(-alpha * ratio)) / (ratio * (1 - math.exp(-alpha)))\n",
    "    RIEmin = (1 - math.exp(alpha * ratio)) / (ratio * (1 - math.exp(alpha)))\n",
    "\n",
    "    if RIEmax != RIEmin:\n",
    "      BEDROC = (RIE - RIEmin) / (RIEmax - RIEmin)\n",
    "    else:  # numActives = numMol\n",
    "      BEDROC = 1.0\n",
    "  else:\n",
    "    BEDROC = 0.0\n",
    "\n",
    "  return BEDROC\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9997158541728963"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CalcBEDROC(scores, 1, 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  another fingerprint for each ligand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem.AtomPairs import Pairs\n",
    "\n",
    "trains=[]\n",
    "c = 0\n",
    "def Calculate_mod(smiles, N_BITS=512):\n",
    "    m = Chem.MolFromSmiles(smiles)\n",
    "    if m is None:\n",
    "        raise ValueError('SMILES cannot be converted to a RDKit molecules:', smiles)\n",
    "\n",
    "    return Pairs.GetAtomPairFingerprint(m)\n",
    "for smile in df['smiles']:\n",
    "    temp = Calculate_mod(smile)\n",
    "#    strFP = temp.ToBitString()\n",
    " #   trains.append(strFP)\n"
   ]
  }
 ],
 "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.8"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}