assay_prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Random Forest with Fingerprint\n",
    "\n",
    "In this baseline model, we want to use compound's fingerprint to predict the outcome of assay. This is a single task, and the prediction performance is measured by the test set with cross-validation.\n",
    "\n",
    "The fingerprint is extracted from all compounds with length `1024`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle\n",
    "from json import load, dump\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from collections import Counter\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from sklearn.model_selection import cross_validate, StratifiedKFold, GridSearchCV\n",
    "from sklearn.linear_model import LogisticRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9404, 141)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the output matrix\n",
    "output_matrix_npz = np.load(\"./resource/output_matrix.npz\")\n",
    "compound_inchi = output_matrix_npz[\"compound_inchi\"]\n",
    "compound_broad_id = output_matrix_npz[\"compound_broad_id\"]\n",
    "output_matrix = output_matrix_npz[\"output_matrix\"]\n",
    "output_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the fingerprint\n",
    "fingerprint = np.load(\"./resource/fp_features.npz\")\n",
    "\n",
    "# Build the feature array\n",
    "fp_feature = np.zeros((len(compound_broad_id), fingerprint[\"features\"].shape[1]))\n",
    "\n",
    "# Choose selected compounds and arrange them by the output matrix order\n",
    "fp_index = dict(zip(fingerprint[\"names\"].astype(str), range(len(fingerprint[\"names\"]))))\n",
    "\n",
    "for b in range(len(compound_broad_id)):\n",
    "    bid = compound_broad_id[b]\n",
    "    fp_feature[b, :] = fingerprint[\"features\"][fp_index[bid], :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It took one hour to create this $9404 \\times 1024$ feature matrix with finger print on condor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9404, 1024)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fp_feature = np.load(\"./resource/extracted_fp_feature.npz\")[\"feature\"]\n",
    "fp_feature.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the PriA-SSB study, some good parameters for a random forests are selected with a different assay dataset. We expect the best parameters (`RF_h`) to perform well in this dataset.\n",
    "\n",
    "For each assay, we only train and predict on compounds that give non-NA results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](https://user-images.githubusercontent.com/3532898/46812713-e8a93280-cd3a-11e8-949a-3469236a3943.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_rf_on_assay(assay_array, rf, kfold=5, n_jobs=4):\n",
    "    \"\"\"\n",
    "    Train and measure a random forest model in a cross validation fasshion.\n",
    "    The model is only trained on compounds which give non-NA results.\n",
    "    \"\"\"\n",
    "    # Filter out NA from the assay_array\n",
    "    y_index = [False if a == -1 else True for a in assay_array ]\n",
    "    y = assay_array[y_index]\n",
    "    x = fp_feature[y_index, :]\n",
    "    \n",
    "    if Counter(y)[1] < kfold:\n",
    "        print(\"Warning: the number of total postives is less than kfold\")\n",
    "        \n",
    "    # One-hot encode y array\n",
    "    #y = np.vstack([[1, 0] if i == 1 else [0, 1] for i in y])\n",
    "    \n",
    "    # Build the cross validation scheme\n",
    "    # Since some assays have extremely small number of positives,\n",
    "    # I will use StratifiedKFold to preserve the proportion of positive in test set\n",
    "    \n",
    "    my_kfold_gen = StratifiedKFold(n_splits=kfold, shuffle=True)\n",
    "    scoring = [\"f1\", \"accuracy\", \"precision\", \"recall\", \"average_precision\", \"roc_auc\"]\n",
    "    cv = cross_validate(rf, x, y, scoring=scoring, cv=my_kfold_gen, n_jobs=n_jobs,\n",
    "                        return_train_score=False)\n",
    "    \n",
    "    cv[\"total_count\"] = Counter(y)\n",
    "    \n",
    "    return cv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {}\n",
    "\n",
    "rf_classifier = RandomForestClassifier(n_estimators=8000, max_features=\"log2\",\n",
    "                                       min_samples_leaf=1, class_weight=\"balanced\")\n",
    "\n",
    "for i in range(output_matrix.shape[1]):\n",
    "    print(i)\n",
    "    results[i] = train_rf_on_assay(output_matrix[:,i], rf_classifier, kfold=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do a stratified 5-fold cross validation for each assay and record metrics on 5 test sets.\n",
    "\n",
    "It took 2 hour to train these 141 single tasks. We can visualize the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = np.load(\"./resource/results.npz\")['results'].item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mean_score(results):\n",
    "    \"\"\"\n",
    "    Aggregate each score over 5 test sets.\n",
    "    \"\"\"\n",
    "    mean_df = {\n",
    "        \"f1\": [],\n",
    "        \"accuracy\": [],\n",
    "        \"average_precision\": [],\n",
    "        \"roc_auc\": [],\n",
    "        \"precision\": [],\n",
    "        \"recall\": [],\n",
    "        \"pos_num\": [],\n",
    "        \"neg_num\": []\n",
    "    }\n",
    "    \n",
    "    for k, r in results.items():\n",
    "        mean_df[\"f1\"].append(np.mean(r[\"test_f1\"]))\n",
    "        mean_df[\"accuracy\"].append(np.mean(r[\"test_accuracy\"]))\n",
    "        mean_df[\"average_precision\"].append(np.mean(r[\"test_average_precision\"]))\n",
    "        mean_df[\"roc_auc\"].append(np.mean(r[\"test_roc_auc\"]))\n",
    "        mean_df[\"precision\"].append(np.mean(r[\"test_precision\"]))\n",
    "        mean_df[\"recall\"].append(np.mean(r[\"test_recall\"]))\n",
    "        mean_df[\"pos_num\"].append(r[\"total_count\"][1])\n",
    "        mean_df[\"neg_num\"].append(r[\"total_count\"][0])\n",
    "    \n",
    "    return pd.DataFrame(mean_df)\n"
   ]
  },
  {
   "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>f1</th>\n",
       "      <th>accuracy</th>\n",
       "      <th>average_precision</th>\n",
       "      <th>roc_auc</th>\n",
       "      <th>precision</th>\n",
       "      <th>recall</th>\n",
       "      <th>pos_num</th>\n",
       "      <th>neg_num</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.080000</td>\n",
       "      <td>0.989659</td>\n",
       "      <td>0.249031</td>\n",
       "      <td>0.753840</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>0.050000</td>\n",
       "      <td>17</td>\n",
       "      <td>1530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.759369</td>\n",
       "      <td>0.682051</td>\n",
       "      <td>0.781926</td>\n",
       "      <td>0.746656</td>\n",
       "      <td>0.676644</td>\n",
       "      <td>0.867821</td>\n",
       "      <td>196</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.183729</td>\n",
       "      <td>0.781891</td>\n",
       "      <td>0.520559</td>\n",
       "      <td>0.776232</td>\n",
       "      <td>0.533333</td>\n",
       "      <td>0.121905</td>\n",
       "      <td>74</td>\n",
       "      <td>256</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.158615</td>\n",
       "      <td>0.828373</td>\n",
       "      <td>0.400695</td>\n",
       "      <td>0.698030</td>\n",
       "      <td>0.805556</td>\n",
       "      <td>0.090476</td>\n",
       "      <td>210</td>\n",
       "      <td>932</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.030199</td>\n",
       "      <td>0.930673</td>\n",
       "      <td>0.168273</td>\n",
       "      <td>0.654066</td>\n",
       "      <td>0.400000</td>\n",
       "      <td>0.015692</td>\n",
       "      <td>128</td>\n",
       "      <td>1747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.660623</td>\n",
       "      <td>0.701898</td>\n",
       "      <td>0.792594</td>\n",
       "      <td>0.774347</td>\n",
       "      <td>0.704681</td>\n",
       "      <td>0.622222</td>\n",
       "      <td>180</td>\n",
       "      <td>206</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.992296</td>\n",
       "      <td>0.026350</td>\n",
       "      <td>0.657967</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>24</td>\n",
       "      <td>3221</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.580000</td>\n",
       "      <td>0.470000</td>\n",
       "      <td>0.777778</td>\n",
       "      <td>0.633333</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.733333</td>\n",
       "      <td>11</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.250102</td>\n",
       "      <td>0.678020</td>\n",
       "      <td>0.509820</td>\n",
       "      <td>0.675496</td>\n",
       "      <td>0.553435</td>\n",
       "      <td>0.162791</td>\n",
       "      <td>215</td>\n",
       "      <td>431</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.019755</td>\n",
       "      <td>0.886010</td>\n",
       "      <td>0.264240</td>\n",
       "      <td>0.685354</td>\n",
       "      <td>0.566667</td>\n",
       "      <td>0.010095</td>\n",
       "      <td>396</td>\n",
       "      <td>3078</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         f1  accuracy  average_precision   roc_auc  precision    recall  \\\n",
       "0  0.080000  0.989659           0.249031  0.753840   0.200000  0.050000   \n",
       "1  0.759369  0.682051           0.781926  0.746656   0.676644  0.867821   \n",
       "2  0.183729  0.781891           0.520559  0.776232   0.533333  0.121905   \n",
       "3  0.158615  0.828373           0.400695  0.698030   0.805556  0.090476   \n",
       "4  0.030199  0.930673           0.168273  0.654066   0.400000  0.015692   \n",
       "5  0.660623  0.701898           0.792594  0.774347   0.704681  0.622222   \n",
       "6  0.000000  0.992296           0.026350  0.657967   0.000000  0.000000   \n",
       "7  0.580000  0.470000           0.777778  0.633333   0.500000  0.733333   \n",
       "8  0.250102  0.678020           0.509820  0.675496   0.553435  0.162791   \n",
       "9  0.019755  0.886010           0.264240  0.685354   0.566667  0.010095   \n",
       "\n",
       "   pos_num  neg_num  \n",
       "0       17     1530  \n",
       "1      196      143  \n",
       "2       74      256  \n",
       "3      210      932  \n",
       "4      128     1747  \n",
       "5      180      206  \n",
       "6       24     3221  \n",
       "7       11       10  \n",
       "8      215      431  \n",
       "9      396     3078  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_df = get_mean_score(results)\n",
    "mean_df.to_csv(\"./mean_df.csv\", index=False)\n",
    "mean_df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Across 141 assays, the average f1=10.03%, accuracy=91.96%, ap=30.35%, auc=72.19%, precision=27.76%, recall=7.83%.\n"
     ]
    }
   ],
   "source": [
    "print((\"Across 141 assays, the average f1={:.2f}%, accuracy={:.2f}%, ap={:.2f}%, auc={:.2f}%, \" + \\\n",
    "       \"precision={:.2f}%, recall={:.2f}%.\").format(\n",
    "    np.mean(mean_df[\"f1\"]) * 100,\n",
    "    np.mean(mean_df[\"accuracy\"]) * 100,\n",
    "    np.mean(mean_df[\"average_precision\"]) * 100,\n",
    "    np.mean(mean_df[\"roc_auc\"]) * 100,\n",
    "    np.mean(mean_df[\"precision\"]) * 100,\n",
    "    np.mean(mean_df[\"recall\"]) * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![random_forest_plot_1.png](./plots/random_forest_plot_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. Comments\n",
    "\n",
    "- Accuracies are pretty high except some outliers with few assays. It is due to the skewness of positive samples.\n",
    "- This baseline model performs poorly based on metrics `F1`, `AP`, `Precision`, `Recall`.\n",
    "- However, this model has high `AUC`, which is the main metric used in the ICCR paper.\n",
    "- Using `AUC`, finterprint baseline gives a similar result as ICCR paper's advanced model (attached below).\n",
    "- `AP` has a very large variance. Low sample size (both pos and neg) tends to give high `AP` values.\n",
    "\n",
    "![](https://i.imgur.com/Z790iMB.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. Random Forest + Finger Print using Converted InChI\n",
    "\n",
    "I have built an output matrix for the intersected compounds using converted InChI strings. I have also excluded some problematic entries (conversion collision) from the intersection.\n",
    "\n",
    "This output matrix uses 209 assays (increasing from 141), and 26638 compounds (increasing from 9404). It should give us a more fair comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26638, 209)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the output matrix\n",
    "output_matrix_npz = np.load(\"./resource/output_matrix_convert.npz\")\n",
    "compound_inchi = output_matrix_npz[\"compound_inchi\"]\n",
    "compound_broad_id = output_matrix_npz[\"compound_broad_id\"]\n",
    "output_matrix = output_matrix_npz[\"output_matrix\"]\n",
    "output_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the fingerprint\n",
    "fingerprint = np.load(\"./resource/fp_features.npz\")\n",
    "\n",
    "# Build the feature array\n",
    "fp_feature = np.zeros((len(compound_broad_id), fingerprint[\"features\"].shape[1]))\n",
    "\n",
    "# Choose selected compounds and arrange them by the output matrix order\n",
    "fp_index = dict(zip(fingerprint[\"names\"].astype(str), range(len(fingerprint[\"names\"]))))\n",
    "\n",
    "for b in range(len(compound_broad_id)):\n",
    "    bid = compound_broad_id[b]\n",
    "    fp_feature[b, :] = fingerprint[\"features\"][fp_index[bid], :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(209, 8)\n"
     ]
    },
    {
     "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>f1</th>\n",
       "      <th>accuracy</th>\n",
       "      <th>average_precision</th>\n",
       "      <th>roc_auc</th>\n",
       "      <th>precision</th>\n",
       "      <th>recall</th>\n",
       "      <th>pos_num</th>\n",
       "      <th>neg_num</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.080000</td>\n",
       "      <td>0.988705</td>\n",
       "      <td>0.408457</td>\n",
       "      <td>0.873482</td>\n",
       "      <td>0.400000</td>\n",
       "      <td>0.044444</td>\n",
       "      <td>45</td>\n",
       "      <td>3762</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.762071</td>\n",
       "      <td>0.689145</td>\n",
       "      <td>0.785920</td>\n",
       "      <td>0.742719</td>\n",
       "      <td>0.682989</td>\n",
       "      <td>0.864217</td>\n",
       "      <td>383</td>\n",
       "      <td>283</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.998618</td>\n",
       "      <td>0.168920</td>\n",
       "      <td>0.628251</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>13</td>\n",
       "      <td>9395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.712121</td>\n",
       "      <td>0.591429</td>\n",
       "      <td>0.643214</td>\n",
       "      <td>0.411111</td>\n",
       "      <td>0.604286</td>\n",
       "      <td>0.900000</td>\n",
       "      <td>18</td>\n",
       "      <td>14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.818182</td>\n",
       "      <td>0.222150</td>\n",
       "      <td>0.388889</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>10</td>\n",
       "      <td>45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.273907</td>\n",
       "      <td>0.784866</td>\n",
       "      <td>0.512232</td>\n",
       "      <td>0.736592</td>\n",
       "      <td>0.660476</td>\n",
       "      <td>0.175484</td>\n",
       "      <td>154</td>\n",
       "      <td>506</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.064461</td>\n",
       "      <td>0.865974</td>\n",
       "      <td>0.391856</td>\n",
       "      <td>0.733942</td>\n",
       "      <td>0.866667</td>\n",
       "      <td>0.034040</td>\n",
       "      <td>383</td>\n",
       "      <td>2400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.012698</td>\n",
       "      <td>0.930264</td>\n",
       "      <td>0.195145</td>\n",
       "      <td>0.671045</td>\n",
       "      <td>0.300000</td>\n",
       "      <td>0.006504</td>\n",
       "      <td>309</td>\n",
       "      <td>4122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.502233</td>\n",
       "      <td>0.685106</td>\n",
       "      <td>0.675068</td>\n",
       "      <td>0.728110</td>\n",
       "      <td>0.657874</td>\n",
       "      <td>0.407571</td>\n",
       "      <td>297</td>\n",
       "      <td>462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.992969</td>\n",
       "      <td>0.056756</td>\n",
       "      <td>0.633252</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>58</td>\n",
       "      <td>8333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         f1  accuracy  average_precision   roc_auc  precision    recall  \\\n",
       "0  0.080000  0.988705           0.408457  0.873482   0.400000  0.044444   \n",
       "1  0.762071  0.689145           0.785920  0.742719   0.682989  0.864217   \n",
       "2  0.000000  0.998618           0.168920  0.628251   0.000000  0.000000   \n",
       "3  0.712121  0.591429           0.643214  0.411111   0.604286  0.900000   \n",
       "4  0.000000  0.818182           0.222150  0.388889   0.000000  0.000000   \n",
       "5  0.273907  0.784866           0.512232  0.736592   0.660476  0.175484   \n",
       "6  0.064461  0.865974           0.391856  0.733942   0.866667  0.034040   \n",
       "7  0.012698  0.930264           0.195145  0.671045   0.300000  0.006504   \n",
       "8  0.502233  0.685106           0.675068  0.728110   0.657874  0.407571   \n",
       "9  0.000000  0.992969           0.056756  0.633252   0.000000  0.000000   \n",
       "\n",
       "   pos_num  neg_num  \n",
       "0       45     3762  \n",
       "1      383      283  \n",
       "2       13     9395  \n",
       "3       18       14  \n",
       "4       10       45  \n",
       "5      154      506  \n",
       "6      383     2400  \n",
       "7      309     4122  \n",
       "8      297      462  \n",
       "9       58     8333  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = np.load(\"./resource/results_converted.npz\")['results'].item()\n",
    "mean_df = get_mean_score(results)\n",
    "mean_df.to_csv(\"./mean_df_converted.csv\", index=False)\n",
    "print(mean_df.shape)\n",
    "mean_df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Across 209 assays, the average f1=13.18%, accuracy=90.72%, ap=31.74%, auc=71.22%, precision=30.44%, recall=11.43%.\n"
     ]
    }
   ],
   "source": [
    "print((\"Across 209 assays, the average f1={:.2f}%, accuracy={:.2f}%, ap={:.2f}%, auc={:.2f}%, \" + \\\n",
    "       \"precision={:.2f}%, recall={:.2f}%.\").format(\n",
    "    np.mean(mean_df[\"f1\"]) * 100,\n",
    "    np.mean(mean_df[\"accuracy\"]) * 100,\n",
    "    np.mean(mean_df[\"average_precision\"]) * 100,\n",
    "    np.mean(mean_df[\"roc_auc\"]) * 100,\n",
    "    np.mean(mean_df[\"precision\"]) * 100,\n",
    "    np.mean(mean_df[\"recall\"]) * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./plots/random_forest_plot_2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.1. Comment\n",
    "\n",
    "We can see the same pattern as the random forest results from 141 assays. In conclusion, `AUC` is not a good metrics, and it is good that we have this fingerprint baseline."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Logistic Regression with Inception Feature\n",
    "\n",
    "We can use the Inception V3 extracted features from U20S images to build end-to-end single task models. It can tell us more about our future image-to-assay model.\n",
    "\n",
    "In this section, we will use the data with ~~209~~ 212 assays. Due to the size limit of Gluster, we only have 244 plates. The paper reports that they have 375 plates, and their shared download script implies there are 406 plates. Therefore, we will not have features for all used compounds in the 212 assays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(27241, 212)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the output matrix\n",
    "output_matrix_npz = np.load(\"./resource/output_matrix_convert_collision.npz\")\n",
    "compound_inchi = output_matrix_npz[\"compound_inchi\"]\n",
    "compound_broad_id = output_matrix_npz[\"compound_broad_id\"]\n",
    "output_matrix = output_matrix_npz[\"output_matrix\"]\n",
    "output_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is slow to work with DataFrame, we can make a dictionary with BID -> PID+WID\n",
    "mean_well_df = pd.read_csv(\"./resource/merged_mean_table.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From our 244 plates, we found 17915 compounds over 27241 total intersected compounds.\n"
     ]
    }
   ],
   "source": [
    "all_pids = mean_well_df[\"Metadata_Plate\"].tolist()\n",
    "all_wids = mean_well_df[\"Metadata_Well\"].tolist()\n",
    "all_bids = mean_well_df[\"Metadata_broad_sample\"].tolist()\n",
    "\n",
    "mean_well_dict = {}\n",
    "intersect_bids = set(compound_broad_id)\n",
    "for i in range(len(all_bids)):\n",
    "    cur_bid = all_bids[i]\n",
    "    if cur_bid in intersect_bids:\n",
    "        if cur_bid not in mean_well_dict:\n",
    "            mean_well_dict[cur_bid] = {\"pid\": [all_pids[i]],\n",
    "                                       \"wid\": [all_wids[i]]}\n",
    "        else:\n",
    "            mean_well_dict[cur_bid][\"pid\"].append(all_pids[i])\n",
    "            mean_well_dict[cur_bid][\"wid\"].append(all_wids[i])\n",
    "     \n",
    "print(\"From our 244 plates, we found {} compounds over {} total intersected compounds.\".format(\n",
    "    len(mean_well_dict), len(compound_broad_id)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reorder pids and wids to match output matrix row order\n",
    "ordered_pids, ordered_wids = [], []\n",
    "compound_replicate = {}\n",
    "\n",
    "for b in compound_broad_id:\n",
    "    if b in mean_well_dict:\n",
    "        ordered_pids.extend(mean_well_dict[b][\"pid\"])\n",
    "        ordered_wids.extend(mean_well_dict[b][\"wid\"])\n",
    "        compound_replicate[b] = len(mean_well_dict[b][\"pid\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({3: 2107, 4: 14806, 2: 544, 1: 457, 8: 1})"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Counter(compound_replicate.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our dataset, most compounds have 4 replicates (tested on 4 wells). For each well, we have 6-9 field of view shots, so we have about 24 samples for one compound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract features for those plates & wells\n",
    "pwid = [\"{}_{}\".format(ordered_pids[i], ordered_wids[i]) for i in range(len(ordered_wids))]\n",
    "pwid_set = set(pwid)\n",
    "\n",
    "combined_feature = np.load(\"/Users/JayWong/Downloads/combined_feature.npz\")\n",
    "combined_feature_pwids = combined_feature[\"names\"]\n",
    "\n",
    "# Get matching indices\n",
    "matched_indices = []\n",
    "matched_pwids = []\n",
    "for i in range(len(combined_feature_pwids)):\n",
    "    if combined_feature_pwids[i] in pwid_set:\n",
    "        matched_indices.append(i)\n",
    "        matched_pwids.append(combined_feature_pwids[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_feature = combined_feature[\"features\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "matched_feature = all_feature[matched_indices,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "matched_names = combined_feature_pwids[matched_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savez(\"./resource/matched_collision_raw_features.npz\", features=matched_feature, names=matched_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "matched_features = np.load(\"./resource/matched_184_raw_features.npz\")[\"features\"]\n",
    "matched_names = np.load(\"./resource/matched_184_raw_features.npz\")[\"names\"]\n",
    "output_matrix = np.load(\"./resource/output_matrix_inception.npz\")[\"output_matrix\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(553044, 4096)\n",
      "15192 compounds found in our extracted features (we only extracted 184 plates).\n"
     ]
    }
   ],
   "source": [
    "# Rearrange output matrix to match our feature matrix\n",
    "print(matched_feature.shape)\n",
    "\n",
    "pwd_to_b_dict = {}\n",
    "for k, v in mean_well_dict.items():\n",
    "    for i in range(len(v[\"pid\"])):\n",
    "        pwd = \"{}_{}\".format(v[\"pid\"][i], v[\"wid\"][i])\n",
    "        if pwd in pwd_to_b_dict:\n",
    "            print(pwd)\n",
    "        pwd_to_b_dict[pwd] = k\n",
    "        \n",
    "output_compound_b_to_index = dict(zip(compound_broad_id, range(len(compound_broad_id))))\n",
    "\n",
    "# Create the indices to select entries from output matrix\n",
    "select_bids = []\n",
    "select_indices = []\n",
    "for pw in matched_names:\n",
    "    bid = pwd_to_b_dict[pw]\n",
    "    select_bids.append(bid)\n",
    "    select_indices.append(output_compound_b_to_index[bid])\n",
    "\n",
    "print(\"{} compounds found in our extracted features (we only extracted 184 plates).\".format(len(set(select_bids))))\n",
    "\n",
    "# Rearrange output matrix\n",
    "rearranged_output_matrix = output_matrix[select_indices,:]\n",
    "rearranged_output_matrix.shape\n",
    "np.savez(\"./resource/output_matrix_collision_inception.npz\", output_matrix=rearranged_output_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_lr_on_assay(aid, kfold=5, n_jobs=4):\n",
    "    \"\"\"\n",
    "    (1). Use cross validation to find the best regularization parameter\n",
    "    (2). Report the best parameter, and scores on the test set\n",
    "\n",
    "    The model is only trained on compounds which give non-NA results.\n",
    "    \"\"\"\n",
    "    assay_array = rearranged_output_matrix[:, aid]\n",
    "\n",
    "    # Filter out NA from the assay_array\n",
    "    y_index = [False if a == -1 else True for a in assay_array ]\n",
    "    y = assay_array[y_index]\n",
    "    x = matched_features[y_index, :]\n",
    "    \n",
    "    print(x.shape, y.shape)\n",
    "    return\n",
    "    \n",
    "    if Counter(y)[1] < kfold:\n",
    "        print(\"Warning: the number of total postives is less than kfold\")\n",
    "        \n",
    "    \n",
    "    # Check if the training samples are too small\n",
    "    if Counter(y)[1] < 10 or Counter(y)[0] < 10:\n",
    "        print(\"Not enough training samples for assay {}\".formt(aid))\n",
    "        return\n",
    "    \n",
    "    # Build the cross validation scheme\n",
    "    # Since some assays have extremely small number of positives,\n",
    "    # I will use StratifiedKFold to preserve the proportion of positive in test set\n",
    "    \n",
    "    my_kfold_gen = StratifiedKFold(n_splits=kfold, shuffle=True)\n",
    "    scoring = [\"f1\", \"accuracy\", \"precision\", \"recall\", \"average_precision\", \"roc_auc\"]\n",
    "    gs = GridSearchCV(LogisticRegression(penalty='l1', verbose=0,\n",
    "                                         solver='liblinear'),\n",
    "                      {'C': [0.01, 0.1, 1, 10]},\n",
    "                      cv=my_kfold_gen,\n",
    "                      scoring=scoring,\n",
    "                      refit=False,\n",
    "                      n_jobs=n_jobs, verbose=2)\n",
    "    \n",
    "    gs.fit(x, y)\n",
    "    result = gs.cv_results_\n",
    "    result[\"count\"] = Counter(y)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It takes 12 hours to train 3 assays due to the size of training data ($546690 \\times 4096$).\n",
    "\n",
    "- We can use L1 logistic regression to select features first (reduce to 512 features from 4096)\n",
    "- We also can use UMAP to do the feature selection\n",
    "- Set `max_iter` to a smaller number (100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "assay_array = output_matrix[:, 5].copy()\n",
    "\n",
    "# Filter out NA from the assay_array\n",
    "y_index = [False if a == -1 else True for a in assay_array ]\n",
    "y = assay_array[y_index].copy()\n",
    "x = matched_features[y_index, :].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = LogisticRegression(penalty='l1', verbose=1, solver='saga', max_iter=1000)\n",
    "lr.fit(x,y)\n",
    "\n",
    "abs_coef = np.abs(lr.coef_[0])\n",
    "largest_coef_index = np.argpartition(abs_coef, -512)[-512:]\n",
    "\n",
    "np.savez(\"feature_index.npz\", index=largest_coef_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We finished training with 212 assays (ignoring collision) with full 4096 features using Condor batch submission. There are 4 assays jobs requiring 200+GB memory and more time to train.\n",
    "\n",
    "- `SAGA` solver works better than `Liblinear`\n",
    "- Across 212-4=208 assays, the average f1=34.87%, accuracy=84.48%, ap=33.12%, auc=85.22%, precision=29.77%, recall=75.17%.\n",
    "\n",
    "![](./plots/lr_plot_1.png)\n",
    "\n",
    "- The average performance is better than fingerprint random forest.\n",
    "- Accuracy has a larger variance (matching AUC).\n",
    "- Similar to the fingerprint random forest, this model struggles with high positive size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Logistic Regression with Inception Features for 349 Plates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(27241, 212)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the output matrix\n",
    "output_matrix_npz = np.load(\"./resource/output_matrix_convert_collision.npz\")\n",
    "compound_inchi = output_matrix_npz[\"compound_inchi\"]\n",
    "compound_broad_id = output_matrix_npz[\"compound_broad_id\"]\n",
    "output_matrix = output_matrix_npz[\"output_matrix\"]\n",
    "output_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is slow to work with DataFrame, we can make a dictionary with BID -> PID+WID\n",
    "mean_well_df = pd.read_csv(\"./resource/merged_mean_table_349.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From our 349 plates, we found 24857 compounds over 27241 total intersected compounds.\n"
     ]
    }
   ],
   "source": [
    "all_pids = mean_well_df[\"Metadata_Plate\"].tolist()\n",
    "all_wids = mean_well_df[\"Metadata_Well\"].tolist()\n",
    "all_bids = mean_well_df[\"Metadata_broad_sample\"].tolist()\n",
    "\n",
    "mean_well_dict = {}\n",
    "intersect_bids = set(compound_broad_id)\n",
    "for i in range(len(all_bids)):\n",
    "    cur_bid = all_bids[i]\n",
    "    if cur_bid in intersect_bids:\n",
    "        if cur_bid not in mean_well_dict:\n",
    "            mean_well_dict[cur_bid] = {\"pid\": [all_pids[i]],\n",
    "                                       \"wid\": [all_wids[i]]}\n",
    "        else:\n",
    "            mean_well_dict[cur_bid][\"pid\"].append(all_pids[i])\n",
    "            mean_well_dict[cur_bid][\"wid\"].append(all_wids[i])\n",
    "     \n",
    "print(\"From our 349 plates, we found {} compounds over {} total intersected compounds.\".format(\n",
    "    len(mean_well_dict), len(compound_broad_id)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reorder pids and wids to match output matrix row order\n",
    "ordered_pids, ordered_wids = [], []\n",
    "compound_replicate = {}\n",
    "\n",
    "for b in compound_broad_id:\n",
    "    if b in mean_well_dict:\n",
    "        ordered_pids.extend(mean_well_dict[b][\"pid\"])\n",
    "        ordered_wids.extend(mean_well_dict[b][\"wid\"])\n",
    "        compound_replicate[b] = len(mean_well_dict[b][\"pid\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({3: 3041, 4: 19680, 8: 768, 1: 972, 2: 392, 7: 1, 6: 3})"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Counter(compound_replicate.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Logistic Regression with Normalized Inception Features for 406 Plates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(27241, 212)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the output matrix\n",
    "output_matrix_npz = np.load(\"./resource/output_matrix_convert_collision.npz\")\n",
    "compound_inchi = output_matrix_npz[\"compound_inchi\"]\n",
    "compound_broad_id = output_matrix_npz[\"compound_broad_id\"]\n",
    "output_matrix = output_matrix_npz[\"output_matrix\"]\n",
    "output_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is slow to work with DataFrame, we can make a dictionary with BID -> PID+WID\n",
    "mean_well_df = pd.read_csv(\"./resource/merged_table_406.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From all 406 plates, we found 26939 compounds over 27241 total intersected compounds.\n"
     ]
    }
   ],
   "source": [
    "all_pids = mean_well_df[\"Metadata_Plate\"].tolist()\n",
    "all_wids = mean_well_df[\"Metadata_Well\"].tolist()\n",
    "all_bids = mean_well_df[\"Metadata_broad_sample\"].tolist()\n",
    "\n",
    "mean_well_dict = {}\n",
    "intersect_bids = set(compound_broad_id)\n",
    "for i in range(len(all_bids)):\n",
    "    cur_bid = all_bids[i]\n",
    "    if cur_bid in intersect_bids:\n",
    "        if cur_bid not in mean_well_dict:\n",
    "            mean_well_dict[cur_bid] = {\"pid\": [all_pids[i]],\n",
    "                                       \"wid\": [all_wids[i]]}\n",
    "        else:\n",
    "            mean_well_dict[cur_bid][\"pid\"].append(all_pids[i])\n",
    "            mean_well_dict[cur_bid][\"wid\"].append(all_wids[i])\n",
    "     \n",
    "print(\"From all 406 plates, we found {} compounds over {} total intersected compounds.\".format(\n",
    "    len(mean_well_dict), len(compound_broad_id)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_well_df = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reorder pids and wids to match output matrix row order\n",
    "ordered_pids, ordered_wids = [], []\n",
    "compound_replicate = {}\n",
    "\n",
    "for b in compound_broad_id:\n",
    "    if b in mean_well_dict:\n",
    "        ordered_pids.extend(mean_well_dict[b][\"pid\"])\n",
    "        ordered_wids.extend(mean_well_dict[b][\"wid\"])\n",
    "        compound_replicate[b] = len(mean_well_dict[b][\"pid\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({3: 2558, 4: 22026, 8: 1522, 7: 238, 2: 392, 1: 200, 6: 3})"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Counter(compound_replicate.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "659267 images out of 913930 images are used\n"
     ]
    }
   ],
   "source": [
    "# Extract features for those plates & wells\n",
    "pwid = [\"{}_{}\".format(ordered_pids[i], ordered_wids[i]) for i in range(len(ordered_wids))]\n",
    "pwid_set = set(pwid)\n",
    "\n",
    "combined_feature = np.load(\"./resource/normed_features.npz\")\n",
    "combined_feature_pwids = combined_feature[\"names\"]\n",
    "\n",
    "# Get matching indices\n",
    "matched_indices = []\n",
    "matched_pwids = []\n",
    "for i in range(len(combined_feature_pwids)):\n",
    "    if combined_feature_pwids[i] in pwid_set:\n",
    "        matched_indices.append(i)\n",
    "        matched_pwids.append(combined_feature_pwids[i])\n",
    "        \n",
    "print(\"{} images out of {} images are used\".format(len(matched_indices), len(combined_feature_pwids)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_feature = combined_feature[\"features\"]\n",
    "matched_feature = all_feature[matched_indices,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matched_names = combined_feature_pwids[matched_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savez(\"./resource/matched_collision_raw_features.npz\", features=matched_feature, names=matched_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Use Mean Well Features to Predict Activity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(27241, 212)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the output matrix\n",
    "output_matrix_npz = np.load(\"./resource/output_matrix_convert_collision.npz\")\n",
    "compound_inchi = output_matrix_npz[\"compound_inchi\"]\n",
    "compound_broad_id = output_matrix_npz[\"compound_broad_id\"]\n",
    "output_matrix = output_matrix_npz[\"output_matrix\"]\n",
    "output_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is slow to work with DataFrame, we can make a dictionary with BID -> PID+WID\n",
    "mean_well_df = pd.read_csv(\"./resource/merged_table_406.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is too slow to operate on the dataframe, we read line by line\n",
    "with open(\"./resource/merged_table_406.csv\", 'r') as fp:\n",
    "    for line in fp:\n",
    "        attrs = fp.split(',')\n",
    "        pid = int(attrs[0])\n",
    "        wid = attrs[1]\n",
    "        bid = attrs[6]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From all 406 plates, we found 26939 compounds over 27241 total intersected compounds.\n"
     ]
    }
   ],
   "source": [
    "all_pids = mean_well_df[\"Metadata_Plate\"].tolist()\n",
    "all_wids = mean_well_df[\"Metadata_Well\"].tolist()\n",
    "all_bids = mean_well_df[\"Metadata_broad_sample\"].tolist()\n",
    "\n",
    "mean_well_dict = {}\n",
    "intersect_bids = set(compound_broad_id)\n",
    "for i in range(len(all_bids)):\n",
    "    cur_bid = all_bids[i]\n",
    "    if cur_bid in intersect_bids:\n",
    "        if cur_bid not in mean_well_dict:\n",
    "            mean_well_dict[cur_bid] = {\"pid\": [all_pids[i]],\n",
    "                                       \"wid\": [all_wids[i]]}\n",
    "        else:\n",
    "            mean_well_dict[cur_bid][\"pid\"].append(all_pids[i])\n",
    "            mean_well_dict[cur_bid][\"wid\"].append(all_wids[i])\n",
    "     \n",
    "print(\"From all 406 plates, we found {} compounds over {} total intersected compounds.\".format(\n",
    "    len(mean_well_dict), len(compound_broad_id)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reorder pids and wids to match output matrix row order\n",
    "ordered_pids, ordered_wids = [], []\n",
    "compound_replicate = {}\n",
    "\n",
    "for b in compound_broad_id:\n",
    "    if b in mean_well_dict:\n",
    "        ordered_pids.extend(mean_well_dict[b][\"pid\"])\n",
    "        ordered_wids.extend(mean_well_dict[b][\"wid\"])\n",
    "        compound_replicate[b] = len(mean_well_dict[b][\"pid\"])\n",
    "\n",
    "ordered_pid_wid_pairs = set(zip(ordered_pids, ordered_wids))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract overlapping rows from the well df\n",
    "matched_indices = []\n",
    "\n",
    "for i, r in mean_well_df.iterrows():\n",
    "    if i == 20:\n",
    "        print(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is too slow to operate on the dataframe, we read line by line\n",
    "with open(\"./resource/merged_table_406.csv\", 'r') as fp:\n",
    "    with open(\"./resource/merged_table_406_intersect.csv\", 'w') as fo:\n",
    "        fo.write(fp.readline())\n",
    "        for line in fp:\n",
    "            attrs = line.split(',')\n",
    "            pid = int(attrs[0])\n",
    "            wid = attrs[1]\n",
    "            pid_wid_pair = (pid, wid)\n",
    "            if pid_wid_pair in ordered_pid_wid_pairs:\n",
    "                fo.write(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_well_df = pd.read_csv(\"./resource/merged_table_406_intersect.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>Metadata_Plate</th>\n",
       "      <th>Metadata_Well</th>\n",
       "      <th>Metadata_Assay_Plate_Barcode</th>\n",
       "      <th>Metadata_Plate_Map_Name</th>\n",
       "      <th>Metadata_well_position</th>\n",
       "      <th>Metadata_ASSAY_WELL_ROLE</th>\n",
       "      <th>Metadata_broad_sample</th>\n",
       "      <th>Metadata_mmoles_per_liter</th>\n",
       "      <th>Metadata_solvent</th>\n",
       "      <th>Metadata_pert_id</th>\n",
       "      <th>...</th>\n",
       "      <th>Nuclei_Texture_Variance_DNA_5_0</th>\n",
       "      <th>Nuclei_Texture_Variance_ER_10_0</th>\n",
       "      <th>Nuclei_Texture_Variance_ER_3_0</th>\n",
       "      <th>Nuclei_Texture_Variance_ER_5_0</th>\n",
       "      <th>Nuclei_Texture_Variance_Mito_10_0</th>\n",
       "      <th>Nuclei_Texture_Variance_Mito_3_0</th>\n",
       "      <th>Nuclei_Texture_Variance_Mito_5_0</th>\n",
       "      <th>Nuclei_Texture_Variance_RNA_10_0</th>\n",
       "      <th>Nuclei_Texture_Variance_RNA_3_0</th>\n",
       "      <th>Nuclei_Texture_Variance_RNA_5_0</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>25855</td>\n",
       "      <td>a01</td>\n",
       "      <td>25855</td>\n",
       "      <td>H-CBLF-002-4</td>\n",
       "      <td>a01</td>\n",
       "      <td>treated</td>\n",
       "      <td>BRD-K14087339-001-01-6</td>\n",
       "      <td>4.759276</td>\n",
       "      <td>DMSO</td>\n",
       "      <td>BRD-K14087339</td>\n",
       "      <td>...</td>\n",
       "      <td>3.377425</td>\n",
       "      <td>1.258327</td>\n",
       "      <td>1.328166</td>\n",
       "      <td>1.288077</td>\n",
       "      <td>1.436918</td>\n",
       "      <td>1.645024</td>\n",
       "      <td>1.571905</td>\n",
       "      <td>2.358518</td>\n",
       "      <td>2.383887</td>\n",
       "      <td>2.380844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>25855</td>\n",
       "      <td>a02</td>\n",
       "      <td>25855</td>\n",
       "      <td>H-CBLF-002-4</td>\n",
       "      <td>a02</td>\n",
       "      <td>treated</td>\n",
       "      <td>BRD-K53903148-001-01-7</td>\n",
       "      <td>4.924897</td>\n",
       "      <td>DMSO</td>\n",
       "      <td>BRD-K53903148</td>\n",
       "      <td>...</td>\n",
       "      <td>3.349516</td>\n",
       "      <td>1.308108</td>\n",
       "      <td>1.367728</td>\n",
       "      <td>1.343064</td>\n",
       "      <td>1.488288</td>\n",
       "      <td>1.612786</td>\n",
       "      <td>1.543314</td>\n",
       "      <td>3.132551</td>\n",
       "      <td>2.955101</td>\n",
       "      <td>2.992260</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>25855</td>\n",
       "      <td>a03</td>\n",
       "      <td>25855</td>\n",
       "      <td>H-CBLF-002-4</td>\n",
       "      <td>a03</td>\n",
       "      <td>treated</td>\n",
       "      <td>BRD-K37357048-001-01-8</td>\n",
       "      <td>5.138267</td>\n",
       "      <td>DMSO</td>\n",
       "      <td>BRD-K37357048</td>\n",
       "      <td>...</td>\n",
       "      <td>3.815022</td>\n",
       "      <td>1.584991</td>\n",
       "      <td>1.748497</td>\n",
       "      <td>1.705850</td>\n",
       "      <td>1.439318</td>\n",
       "      <td>1.577664</td>\n",
       "      <td>1.529118</td>\n",
       "      <td>3.023558</td>\n",
       "      <td>2.832276</td>\n",
       "      <td>2.911847</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>25855</td>\n",
       "      <td>a04</td>\n",
       "      <td>25855</td>\n",
       "      <td>H-CBLF-002-4</td>\n",
       "      <td>a04</td>\n",
       "      <td>treated</td>\n",
       "      <td>BRD-K25385069-001-01-7</td>\n",
       "      <td>4.916720</td>\n",
       "      <td>DMSO</td>\n",
       "      <td>BRD-K25385069</td>\n",
       "      <td>...</td>\n",
       "      <td>3.486581</td>\n",
       "      <td>1.385492</td>\n",
       "      <td>1.472731</td>\n",
       "      <td>1.435720</td>\n",
       "      <td>1.588143</td>\n",
       "      <td>1.813048</td>\n",
       "      <td>1.776891</td>\n",
       "      <td>3.168433</td>\n",
       "      <td>2.967504</td>\n",
       "      <td>3.046684</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>25855</td>\n",
       "      <td>a05</td>\n",
       "      <td>25855</td>\n",
       "      <td>H-CBLF-002-4</td>\n",
       "      <td>a05</td>\n",
       "      <td>treated</td>\n",
       "      <td>BRD-K63140065-001-01-3</td>\n",
       "      <td>5.226743</td>\n",
       "      <td>DMSO</td>\n",
       "      <td>BRD-K63140065</td>\n",
       "      <td>...</td>\n",
       "      <td>4.170932</td>\n",
       "      <td>1.832773</td>\n",
       "      <td>1.935458</td>\n",
       "      <td>1.886178</td>\n",
       "      <td>1.626151</td>\n",
       "      <td>1.706171</td>\n",
       "      <td>1.688343</td>\n",
       "      <td>2.773610</td>\n",
       "      <td>2.698169</td>\n",
       "      <td>2.705382</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 1800 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   Metadata_Plate Metadata_Well  Metadata_Assay_Plate_Barcode  \\\n",
       "0           25855           a01                         25855   \n",
       "1           25855           a02                         25855   \n",
       "2           25855           a03                         25855   \n",
       "3           25855           a04                         25855   \n",
       "4           25855           a05                         25855   \n",
       "\n",
       "  Metadata_Plate_Map_Name Metadata_well_position Metadata_ASSAY_WELL_ROLE  \\\n",
       "0            H-CBLF-002-4                    a01                  treated   \n",
       "1            H-CBLF-002-4                    a02                  treated   \n",
       "2            H-CBLF-002-4                    a03                  treated   \n",
       "3            H-CBLF-002-4                    a04                  treated   \n",
       "4            H-CBLF-002-4                    a05                  treated   \n",
       "\n",
       "    Metadata_broad_sample  Metadata_mmoles_per_liter Metadata_solvent  \\\n",
       "0  BRD-K14087339-001-01-6                   4.759276             DMSO   \n",
       "1  BRD-K53903148-001-01-7                   4.924897             DMSO   \n",
       "2  BRD-K37357048-001-01-8                   5.138267             DMSO   \n",
       "3  BRD-K25385069-001-01-7                   4.916720             DMSO   \n",
       "4  BRD-K63140065-001-01-3                   5.226743             DMSO   \n",
       "\n",
       "  Metadata_pert_id  ... Nuclei_Texture_Variance_DNA_5_0  \\\n",
       "0    BRD-K14087339  ...                        3.377425   \n",
       "1    BRD-K53903148  ...                        3.349516   \n",
       "2    BRD-K37357048  ...                        3.815022   \n",
       "3    BRD-K25385069  ...                        3.486581   \n",
       "4    BRD-K63140065  ...                        4.170932   \n",
       "\n",
       "  Nuclei_Texture_Variance_ER_10_0  Nuclei_Texture_Variance_ER_3_0  \\\n",
       "0                        1.258327                        1.328166   \n",
       "1                        1.308108                        1.367728   \n",
       "2                        1.584991                        1.748497   \n",
       "3                        1.385492                        1.472731   \n",
       "4                        1.832773                        1.935458   \n",
       "\n",
       "  Nuclei_Texture_Variance_ER_5_0 Nuclei_Texture_Variance_Mito_10_0  \\\n",
       "0                       1.288077                          1.436918   \n",
       "1                       1.343064                          1.488288   \n",
       "2                       1.705850                          1.439318   \n",
       "3                       1.435720                          1.588143   \n",
       "4                       1.886178                          1.626151   \n",
       "\n",
       "  Nuclei_Texture_Variance_Mito_3_0 Nuclei_Texture_Variance_Mito_5_0  \\\n",
       "0                         1.645024                         1.571905   \n",
       "1                         1.612786                         1.543314   \n",
       "2                         1.577664                         1.529118   \n",
       "3                         1.813048                         1.776891   \n",
       "4                         1.706171                         1.688343   \n",
       "\n",
       "   Nuclei_Texture_Variance_RNA_10_0  Nuclei_Texture_Variance_RNA_3_0  \\\n",
       "0                          2.358518                         2.383887   \n",
       "1                          3.132551                         2.955101   \n",
       "2                          3.023558                         2.832276   \n",
       "3                          3.168433                         2.967504   \n",
       "4                          2.773610                         2.698169   \n",
       "\n",
       "   Nuclei_Texture_Variance_RNA_5_0  \n",
       "0                         2.380844  \n",
       "1                         2.992260  \n",
       "2                         2.911847  \n",
       "3                         3.046684  \n",
       "4                         2.705382  \n",
       "\n",
       "[5 rows x 1800 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_well_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "matched_feature = mean_well_df.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_bids = mean_well_df[\"Metadata_broad_sample\"].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(110622, 1800)\n"
     ]
    }
   ],
   "source": [
    "# Rearrange output matrix to match our feature matrix\n",
    "print(matched_feature.shape)\n",
    "output_compound_b_to_index = dict(zip(compound_broad_id, range(len(compound_broad_id))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26939 compounds found in our extracted features.\n"
     ]
    }
   ],
   "source": [
    "# Create the indices to select entries from output matrix\n",
    "select_bids = []\n",
    "select_indices = []\n",
    "for bid in all_bids:\n",
    "    select_bids.append(bid)\n",
    "    select_indices.append(output_compound_b_to_index[bid])\n",
    "\n",
    "print(\"{} compounds found in our extracted features.\".format(len(set(select_bids))))\n",
    "\n",
    "# Rearrange output matrix\n",
    "rearranged_output_matrix = output_matrix[select_indices,:]\n",
    "rearranged_output_matrix.shape\n",
    "np.savez(\"./resource/output_matrix_collision_meanwell.npz\", output_matrix=rearranged_output_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(110622, 212)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rearranged_output_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1787\n"
     ]
    }
   ],
   "source": [
    "# Ditch string (categorical) features\n",
    "float_cols = []\n",
    "for c in range(matched_feature.shape[1]):\n",
    "    col = matched_feature[:, c]\n",
    "    try:\n",
    "        temp = col.astype(np.float64)\n",
    "        float_cols.append(c)\n",
    "    except ValueError:\n",
    "        pass\n",
    "    \n",
    "print(len(float_cols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "float_matrix = matched_feature[:, float_cols]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "float_matrix = float_matrix.astype(np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_cols = [i for i in range(1787) if i != 3]\n",
    "float_matrix = float_matrix[:, selected_cols]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savez(\"meanwell_train_features.npz\", features=float_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Metadata_Well', 'Metadata_Plate_Map_Name',\n",
       "       'Metadata_well_position', 'Metadata_ASSAY_WELL_ROLE',\n",
       "       'Metadata_broad_sample', 'Metadata_solvent', 'Metadata_pert_id',\n",
       "       'Metadata_pert_mfc_id', 'Metadata_pert_well', 'Metadata_cell_id',\n",
       "       'Metadata_broad_sample_type', 'Metadata_pert_vehicle',\n",
       "       'Metadata_pert_type'], dtype='<U51')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# and Metadata_pert_id_vendor\n",
    "np.array(list(mean_well_df))[list(set(range(1800)) - set(float_cols))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(110622, 1786)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "float_matrix.shape"
   ]
  },
  {
   "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.7.2"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}