03_MAP4SVM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Pandarallel will run on 8 workers.\n",
      "INFO: Pandarallel will use Memory file system to transfer data between the main process and workers.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from rdkit import Chem \n",
    "from rdkit.Chem import Lipinski, Descriptors, rdMolDescriptors, AllChem, PandasTools\n",
    "from pandarallel import pandarallel\n",
    "import numpy as np\n",
    "pandarallel.initialize(progress_bar=False)\n",
    "import tmap as tm\n",
    "from map4 import MAP4Calculator\n",
    "import os\n",
    "import joblib\n",
    "from sklearn import metrics\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn import svm\n",
    "from sklearn.metrics import roc_curve\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "folder = \"/data/coconut/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "coconut = pd.read_pickle(folder + \"coconut_prop.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "origins = [\"plants\", \"fungi\", \"bacteria\", \"animal\", \"Homo_sapiens\", \"marine\", \"other\"]\n",
    "def origin(simple_tax):\n",
    "    if simple_tax in origins:\n",
    "        return origins.index(simple_tax)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "coconut[\"origin\"] = coconut.simple_tax.map(origin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm_smiles(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    smiles_n = Chem.MolToSmiles(mol, isomericSmiles=False)\n",
    "    return smiles_n\n",
    "\n",
    "coconut[\"norm_smiles\"] = coconut[\"SMILES\"].map(norm_smiles) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "coconut = coconut.query(\"origin != 6 and origin != 5 and origin != 4 and origin != 3\").sort_values(\"origin\", ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "coconut = coconut.drop_duplicates(subset=\"norm_smiles\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "33821"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(coconut.query(\"origin==0\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60333"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(coconut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15693"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(coconut.query(\"origin==1\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10819"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(coconut.query(\"origin==2\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(coconut.query(\"origin==3\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(coconut.query(\"origin==4\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(coconut.query(\"origin==5\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n"
     ]
    }
   ],
   "source": [
    "coconut[\"Set\"] = \"test\"\n",
    "for i in range(6):\n",
    "    print(i)\n",
    "    coconut_training = coconut.query(\"origin == @i\").sample(frac=0.50, random_state=12)\n",
    "    coconut.loc[coconut_training.index, \"Set\"] = \"training\"\n",
    "\n",
    "\n",
    "pd.to_pickle(coconut, folder + \"coconut_4classifier_less_classes.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def map4_kernel_SVM(a, b=None):\n",
    "    if b is None:\n",
    "        b=a\n",
    "    JS_all_pairs = np.zeros((len(a),len(b)))\n",
    "    for i,fp1 in enumerate(a):\n",
    "        for j,fp2 in enumerate(b):\n",
    "            JS_all_pairs[i,j] = np.float(np.count_nonzero(fp1 == fp2)) / np.float(len(fp1))\n",
    "    return JS_all_pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alice/miniconda3/envs/aipep/lib/python3.6/site-packages/pandas/core/missing.py:47: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  mask = arr == x\n",
      "/home/alice/miniconda3/envs/aipep/lib/python3.6/site-packages/pandas/core/missing.py:58: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  mask |= arr == x\n"
     ]
    }
   ],
   "source": [
    "coconut.replace([\"NaN\", 'NaT'], np.nan, inplace = True)\n",
    "coconut = coconut.dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preparing features for training (X) and test (x)\n",
    "X = np.array(coconut[coconut[\"Set\"] == \"training\"].MAP4.to_list())\n",
    "x = np.array(coconut[coconut[\"Set\"] == \"test\"].MAP4.to_list())\n",
    "# generate prediction for training (Y) and test(y) sets\n",
    "Y = coconut[coconut[\"Set\"] == \"training\"].origin.values\n",
    "y = coconut[coconut[\"Set\"] == \"test\"].origin.values\n",
    "Y = np.array(Y, dtype=np.float)\n",
    "y = np.array(y, dtype=np.float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = ['molecular_weight','number_of_carbons',\\\n",
    "       'number_of_nitrogens', 'number_of_oxygens',\\\n",
    "       'total_atom_number', 'bond_count', 'topoPSA', 'fcsp3', 'HBA',\\\n",
    "       'HBD', 'aLogP']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_clf = make_pipeline(svm.SVC(probability=False, class_weight=\"balanced\", cache_size=1900, random_state=12))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(folder+\"MAP4_classifier_eval.pkl\"):\n",
    "    # SVM optimization for ROC auc\n",
    "    param_grid = {'svc__C': [0.1,1, 10, 100, 1000], \\\n",
    "                  'svc__kernel': [map4_kernel_SVM]}\n",
    "    grid = GridSearchCV(base_clf, param_grid, scoring='balanced_accuracy',\\\n",
    "                        refit=True,verbose=2, n_jobs=7)\n",
    "    grid.fit(X,Y)\n",
    "    print(grid.best_estimator_)\n",
    "    \n",
    "    joblib.dump(grid, folder + \"SVM-MAP4_classifier_eval.pkl\")\n",
    "else:\n",
    "    grid = joblib.load(folder + \"SVM-MAP4_classifier_eval.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_predicted = grid.predict(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.utils import check_matplotlib_support\n",
    "from sklearn.utils.multiclass import unique_labels\n",
    "from sklearn.utils.validation import _deprecate_positional_args\n",
    "from sklearn.base import is_classifier\n",
    "import matplotlib\n",
    "\n",
    "class ConfusionMatrixDisplay:\n",
    "    \"\"\"Confusion Matrix visualization.\n",
    "    It is recommend to use :func:`~sklearn.metrics.plot_confusion_matrix` to\n",
    "    create a :class:`ConfusionMatrixDisplay`. All parameters are stored as\n",
    "    attributes.\n",
    "    Read more in the :ref:`User Guide <visualizations>`.\n",
    "    Parameters\n",
    "    ----------\n",
    "    confusion_matrix : ndarray of shape (n_classes, n_classes)\n",
    "        Confusion matrix.\n",
    "    display_labels : ndarray of shape (n_classes,)\n",
    "        Display labels for plot.\n",
    "    Attributes\n",
    "    ----------\n",
    "    im_ : matplotlib AxesImage\n",
    "        Image representing the confusion matrix.\n",
    "    text_ : ndarray of shape (n_classes, n_classes), dtype=matplotlib Text, \\\n",
    "            or None\n",
    "        Array of matplotlib axes. `None` if `include_values` is false.\n",
    "    ax_ : matplotlib Axes\n",
    "        Axes with confusion matrix.\n",
    "    figure_ : matplotlib Figure\n",
    "        Figure containing the confusion matrix.\n",
    "    \"\"\"\n",
    "    def __init__(self, confusion_matrix, display_labels):\n",
    "        self.confusion_matrix = confusion_matrix\n",
    "        self.display_labels = display_labels\n",
    "\n",
    "    def plot(self, include_values=True, cmap='viridis',\n",
    "             xticks_rotation='horizontal', values_format=None, ax=None):\n",
    "        \"\"\"Plot visualization.\n",
    "        Parameters\n",
    "        ----------\n",
    "        include_values : bool, default=True\n",
    "            Includes values in confusion matrix.\n",
    "        cmap : str or matplotlib Colormap, default='viridis'\n",
    "            Colormap recognized by matplotlib.\n",
    "        xticks_rotation : {'vertical', 'horizontal'} or float, \\\n",
    "                         default='horizontal'\n",
    "            Rotation of xtick labels.\n",
    "        values_format : str, default=None\n",
    "            Format specification for values in confusion matrix. If `None`,\n",
    "            the format specification is '.2g'.\n",
    "        ax : matplotlib axes, default=None\n",
    "            Axes object to plot on. If `None`, a new figure and axes is\n",
    "            created.\n",
    "        Returns\n",
    "        -------\n",
    "        display : :class:`~sklearn.metrics.ConfusionMatrixDisplay`\n",
    "        \"\"\"\n",
    "        check_matplotlib_support(\"ConfusionMatrixDisplay.plot\")\n",
    "        import matplotlib.pyplot as plt\n",
    "\n",
    "        if ax is None:\n",
    "            fig, ax = plt.subplots()\n",
    "        else:\n",
    "            fig = ax.figure\n",
    "\n",
    "        cm = self.confusion_matrix\n",
    "        n_classes = cm.shape[0]\n",
    "        self.im_ = ax.imshow(cm, interpolation='nearest', cmap=cmap, vmin=0, vmax=1,)\n",
    "        self.text_ = None\n",
    "\n",
    "        cmap_min, cmap_max = self.im_.cmap(0), self.im_.cmap(256)\n",
    "\n",
    "        if include_values:\n",
    "            self.text_ = np.empty_like(cm, dtype=object)\n",
    "            if values_format is None:\n",
    "                values_format = '.0%'\n",
    "\n",
    "            # print text with appropriate color depending on background\n",
    "            thresh = (cm.max() + cm.min()) / 2.0\n",
    "            for i, j in product(range(n_classes), range(n_classes)):\n",
    "                color = cmap_max if cm[i, j] < thresh else cmap_min\n",
    "                self.text_[i, j] = ax.text(j, i,\n",
    "                                           format(cm[i, j], values_format),\n",
    "                                           ha=\"center\", va=\"center\",\n",
    "                                           color=color, size= 6 )\n",
    "\n",
    "        \n",
    "        colorbar = fig.colorbar(self.im_, ax=ax)\n",
    "        colorbar.ax.tick_params(labelsize=6) \n",
    "        colorbar.set_ticks([0, .25, 0.5, .75, 1])\n",
    "        colorbar.set_ticklabels(['0%', '25%', '50%','75%', '100%'])\n",
    "        #colorbar.set_label(\"Predicted value/True Value\")\n",
    "        ax.set(xticks=np.arange(n_classes),\n",
    "               yticks=np.arange(n_classes),\n",
    "               xticklabels=self.display_labels,\n",
    "               yticklabels=self.display_labels,\n",
    "               ylabel=\"True label\",\n",
    "               xlabel=\"Predicted label\")\n",
    "\n",
    "        \n",
    "        ax.set_ylim((n_classes - 0.5, -0.5))\n",
    "        plt.setp(ax.get_xticklabels(), rotation=xticks_rotation, fontsize=6)\n",
    "        plt.setp(ax.get_yticklabels(), fontsize=6)\n",
    "        ax.xaxis.get_label().set_fontsize(6)\n",
    "        ax.yaxis.get_label().set_fontsize(6)\n",
    "        \n",
    "        self.figure_ = fig\n",
    "        self.ax_ = ax\n",
    "\n",
    "        return self\n",
    "\n",
    "\n",
    "\n",
    "def plot_confusion_matrix(y_pred, y_true, *, labels=None,\n",
    "                          sample_weight=None, normalize=None,\n",
    "                          display_labels=None, include_values=True,\n",
    "                          xticks_rotation='horizontal',\n",
    "                          values_format=None,\n",
    "                          cmap='viridis', ax=None, colorbar=False):\n",
    "\n",
    "    check_matplotlib_support(\"plot_confusion_matrix\")\n",
    "\n",
    "\n",
    "    cm = confusion_matrix(y_true, y_pred, sample_weight=sample_weight,\n",
    "                          labels=labels, normalize=normalize)\n",
    "\n",
    "    if display_labels is None:\n",
    "        display_labels = labels\n",
    "\n",
    "    disp = ConfusionMatrixDisplay(confusion_matrix=cm,\n",
    "                                  display_labels=display_labels)\n",
    "    return disp.plot(include_values=include_values,\n",
    "                     cmap=cmap, ax=ax, xticks_rotation=xticks_rotation,\n",
    "                     values_format=values_format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "#sns.set_context(\"paper\")\n",
    "cm = 1/2.54  # centimeters in inches\n",
    "fig, ax = plt.subplots(figsize=(6*cm, 4*cm))\n",
    "m = plot_confusion_matrix(np.array(y_predicted, dtype=int), np.array(y, dtype=int), display_labels=[\"plants\", \"fungi\", \"bacteria\"], normalize = \"true\", ax=ax)\n",
    "plt.title(\"MAP4 SVM Classifier\\nNormalized Confusion Matrix\", fontsize=6)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"plots/confusion_matrix_normtrue_right_size.png\", dpi=500)\n",
    "plt.savefig(\"plots/confusion_matrix_normtrue_right_size.svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics.balanced_accuracy_score(y_predicted, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics.matthews_corrcoef(y_predicted, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics.f1_score(y_predicted, y, average=\"micro\")"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}