RUN_Archetypes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RUN: Archetype Prediction\n",
    "\n",
    "This notebook runs the archetype prediction with the shape space as features and manual annotations as target.\n",
    "\n",
    "Hyperparameter optimization was done in `DEV_Archetypes.ipynb`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes\n",
    "----\n",
    "\n",
    "- The paper uses slightly different nomenclature than what is used here:\n",
    "\n",
    "| Paper | Code |\n",
    "| :---: | :--: |\n",
    "| `central rosette cells` | `innerRosetteCells` |\n",
    "| `peripheral rosette cells` | `outerRosetteCells` |\n",
    "| `inter-organ cells` | `betweenRosetteCells` |\n",
    "| `leader cells` | `leaderCells` |\n",
    "\n",
    "\n",
    "- Selection criteria for manual cell annotation\n",
    "    - `innerRosetteCells:` cells directly adjacent to a lumen without contact to the outside of the tissue\n",
    "    - `outerRosetteCells:` cells to the left or right of a lumen with a large area of contact to the outside of the tissue\n",
    "    - `betweenRosetteCells:` cells between two rosettes, both within and on the outside of the tissue\n",
    "    - `leaderCells:` the first few cells from the front, in particular those without a clear apical backwards-polarity\n",
    "    \n",
    " \n",
    "- Manual annotation was originally done with the Fiji `Point` tool and `Ctrl+M` to get measurements\n",
    "    - Results stored in a separate `.csv` file for each class\n",
    "    - This was converted to a single numpy vector (see code below)\n",
    "    - The numeric encoding is as follows:\n",
    "        - `0 : 'unclassified',`\n",
    "        - `1 : 'innerRosetteCells',`\n",
    "        - `2 : 'outerRosetteCells',`\n",
    "        - `3 : 'betweenRosetteCells',`\n",
    "        - `4 : 'leaderCells'`\n",
    "    - This is also stored in the metadata as `archetype_decodedict` (and the inverse as `archetype_encodedict`)\n",
    "\n",
    "\n",
    "- On the IDR, the manual annotations are stored as `.tsv` files like everything else\n",
    "    - The code below allows this data to be loaded and used here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Imports\n",
    "\n",
    "# Generic\n",
    "from __future__ import division\n",
    "import os, sys, pickle\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Specific\n",
    "from sklearn.decomposition import PCA\n",
    "import sklearn.svm as svm\n",
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "# Modules\n",
    "import katachi.utilities.loading as ld"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Prep encoding and decoding dict\n",
    "\n",
    "archetype_decodedict = {0 : 'unclassified',\n",
    "                       1 : 'innerRosetteCells',\n",
    "                       2 : 'outerRosetteCells',\n",
    "                       3 : 'betweenRosetteCells',\n",
    "                       4 : 'leaderCells'}\n",
    "\n",
    "archetype_encodedict = {name:key for key,name \n",
    "                       in archetype_decodedict.iteritems()}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conversion from Fiji Output\n",
    "\n",
    "Conversion of manual annotation files from Fiji into a clean numpy vector.\n",
    "\n",
    "When using data from the IDR, this is not necessary.\n",
    "\n",
    "The filename expected is `<sample_ID>_<archetype>_manual.csv`, where `<archetype>` is one of `innerRosetteCells`, `outerRosetteCells`, `betweenRosetteCells`, `leaderCell`. One file must be present for each archetype.\n",
    "\n",
    "The resulting numpy array is saved as `<sample_ID>_archetype_manual.npy`. It contains all four archetypes numerically encoded now, so there won't be one file for each archetype."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Load metadata\n",
    "\n",
    "rawpath   = r'data\\experimentA\\image_data'\n",
    "rawloader = ld.DataLoader(rawpath, recurse=True)\n",
    "meta, prim_IDs, _ = rawloader.load_dataset(\"stack_metadata.pkl\")\n",
    "print \"Imported metadata for\", len(meta.keys()), \"primordia\""
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Double-check that all manual annotation sets are complete\n",
    "\n",
    "# Report\n",
    "print \"Checking files...\"\n",
    "\n",
    "# For each prim...\n",
    "annotations_found = 0\n",
    "for prim_idx, prim_ID in enumerate(prim_IDs):\n",
    "    \n",
    "    # Find all manual annotation files (if any)\n",
    "    found_manual = False\n",
    "    labels_found = []\n",
    "    for fpath in rawloader.data[prim_ID]:\n",
    "        if '_manual.csv' in fpath:\n",
    "            found_manual = True\n",
    "            annotations_found += 1\n",
    "            labels_found.append( [key for key in archetype_encodedict.keys() if key in fpath][0] )\n",
    "    \n",
    "    # Check if all expected files are present\n",
    "    if found_manual:\n",
    "        if not sorted(labels_found) == sorted(['innerRosetteCells', 'outerRosetteCells',\n",
    "                                               'betweenRosetteCells', 'leaderCells']):\n",
    "            print \"  WARNING -- Annotation set incomplete or incorrect for prim:\", prim_ID\n",
    "            \n",
    "# Report\n",
    "print \"Check completed, found\", annotations_found, \"annotated samples.\""
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "### Perform the conversion\n",
    "\n",
    "# For each prim...\n",
    "print \"\\nConverting annotations...\"\n",
    "for prim_idx, prim_ID in enumerate(prim_IDs):\n",
    "    \n",
    "    # Create multiclass vector\n",
    "    # Note: It's a bit weird not to have the number of cells as a \n",
    "    #       more explicit value in the metadata...\n",
    "    archetypes = np.zeros(meta[prim_ID]['centroids'].shape[0], dtype=np.int8)\n",
    "    \n",
    "    # Find any manual annotation files...\n",
    "    found_manual=False\n",
    "    for fpath in rawloader.data[prim_ID]:\n",
    "        if '_manual.csv' in fpath:\n",
    "            \n",
    "            # Remember\n",
    "            found_manual=True\n",
    "            \n",
    "            # Read the file\n",
    "            with open(fpath, 'r') as infile:\n",
    "                \n",
    "                # Skip header\n",
    "                infile.readline()\n",
    "                \n",
    "                # Get data labels\n",
    "                labels = []\n",
    "                for line in infile.readlines():\n",
    "                    labels.append(int(line.strip().split(',')[1]))\n",
    "                labels = np.array(labels) - 1\n",
    "            \n",
    "            # Encode\n",
    "            archetypes[labels] = archetype_encodedict[ [key for key in archetype_encodedict.keys()\n",
    "                                                      if key in fpath][0] ]\n",
    "    \n",
    "    # If a manual annotation was found...\n",
    "    if found_manual:\n",
    "    \n",
    "        # Save encoded array as separate file\n",
    "        for fpath in rawloader.data[prim_ID]:\n",
    "            if fpath.endswith('_seg.tif'):\n",
    "                outpath = fpath[:-4] + '_archetypes_manual.npy'\n",
    "                np.save(outpath, archetypes)\n",
    "    \n",
    "        # Add the dicts to the metadata\n",
    "        meta[prim_ID]['archetype_decodedict'] = archetype_decodedict\n",
    "        meta[prim_ID]['archetype_encodedict'] = archetype_encodedict\n",
    "\n",
    "        # Save metadata\n",
    "        for fpath in rawloader.data[prim_ID]:\n",
    "            if '_stack_metadata.pkl' in fpath:\n",
    "                with open(fpath, 'wb') as metafile:\n",
    "                    pickle.dump(meta[prim_ID], metafile, pickle.HIGHEST_PROTOCOL)\n",
    "                \n",
    "# Report\n",
    "print \"Conversion complete.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Archetype Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load shape space data from IDR\n",
    "\n",
    "# Path to the data\n",
    "dirpath = r'data/experimentA/extracted_measurements/'\n",
    "\n",
    "# Set whether to use TFOR or CFOR\n",
    "fspace_type = 'TFOR'\n",
    "#fspace_type = 'CFOR'\n",
    "\n",
    "# Load corresponding shape space\n",
    "loader = ld.DataLoaderIDR(dirpath, recurse=True)\n",
    "suffix = 'shape_'+fspace_type+'_raw_measured.tsv'\n",
    "fspace, prim_IDs, fspace_idx = loader.load_dataset(suffix)\n",
    "print \"Imported shape space of shape:\", fspace.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load annotations from IDR\n",
    "\n",
    "archetypes, _, _ = loader.load_dataset(\"archetype_manual_annotations.tsv\")\n",
    "archetypes[np.isnan(archetypes)] = 0\n",
    "print \"Imported manual archetype annotations of shape:\", archetypes.shape"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Load annotations from converted Fiji data (untested!)\n",
    "\n",
    "rawpath = r'data\\experimentA\\image_data'\n",
    "rawloader = ld.DataLoader(rawpath, recurse=True)\n",
    "raw_archetypes, raw_prim_IDs, raw_fspace_idx = rawloader.load_dataset(\"_archetypes_manual.npy\")\n",
    "archetypes = np.zeros(fspace.shape[0])\n",
    "archetypes[np.in1d(np.array(prim_IDs)[fspace_idx], raw_prim_IDs)] = sub_archetypes\n",
    "print \"Imported manual archetype annotations of shape:\", archetypes.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Preprocess data\n",
    "\n",
    "# Reduce the training data to the actually annotated samples\n",
    "fspace_train   = fspace[archetypes!=0]\n",
    "fspace_predict = fspace\n",
    "archetypes_rdy = archetypes[archetypes!=0]\n",
    "\n",
    "# Standardization (beneficial for TFOR according to tests)\n",
    "if fspace_type == 'TFOR':\n",
    "    train_means, train_stds = (np.mean(fspace_train, axis=0), np.std(fspace_train, axis=0))\n",
    "    fspace_train   = (fspace_train   - train_means) / train_stds\n",
    "    fspace_predict = (fspace_predict - train_means) / train_stds\n",
    "\n",
    "# PCA (beneficial for CFOR according to tests)\n",
    "if fspace_type == 'CFOR':\n",
    "    pca = PCA().fit(fspace_train)\n",
    "    fspace_train   = pca.transform(fspace_train)\n",
    "    fspace_predict = pca.transform(fspace_predict)\n",
    "\n",
    "# Report\n",
    "print \"Final shapes:\"\n",
    "print \"  archetypes_rdy: \", archetypes_rdy.shape\n",
    "print \"  fspace_train:   \", fspace_train.shape\n",
    "print \"  fspace_predict: \", fspace_predict.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Run prediction\n",
    "\n",
    "# Prep regressor (Hyperparams optimized in `TEST_archetypes`)\n",
    "if fspace_type == 'TFOR':               \n",
    "    svc = svm.SVC(probability=True, kernel='rbf',\n",
    "                  C=1.0, gamma=1.0 / fspace_train.shape[1]) \n",
    "if fspace_type == 'CFOR':\n",
    "    svc = svm.SVC(probability=True, kernel='rbf',\n",
    "                  C=10.0, gamma=1.0 / fspace_train.shape[1]) \n",
    "\n",
    "# Fit\n",
    "print \"Fitting...\"\n",
    "svc.fit(fspace_train, archetypes_rdy)\n",
    "\n",
    "# Predict\n",
    "print \"Predicting...\"\n",
    "archetypes_pred = svc.predict(fspace_predict)\n",
    "archetypes_prob = svc.predict_proba(fspace_predict)\n",
    "\n",
    "# Report\n",
    "print \"\\nPredicted archetypes_pred of shape:\", archetypes_pred.shape\n",
    "print   \"Predicted archetypes_prob of shape:\", archetypes_prob.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Quick & Simple confusion matrix evaluation on training data\n",
    "\n",
    "# Repredict just the training samples\n",
    "archetypes_train_pred = svc.predict(fspace_train)\n",
    "\n",
    "# Compute confusion matrix\n",
    "cm = confusion_matrix(archetypes_rdy, archetypes_train_pred)\n",
    "\n",
    "# Prep plot\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3,3))\n",
    "\n",
    "# Function for creating and styling a confusion matrix\n",
    "def confmat(ax, cm):\n",
    "    \n",
    "    # Handle axis boundaries...\n",
    "    ax.set_adjustable('box-forced')\n",
    "    \n",
    "    # Show image and text\n",
    "    ax.imshow(cm, interpolation='none', cmap='Blues')\n",
    "    for (i, j), z in np.ndenumerate(cm):\n",
    "        ax.text(j, i, z, ha='center', va='center')\n",
    "    \n",
    "    # Adjust ticks\n",
    "    ax.set_xticks(range(cm.shape[0]))\n",
    "    ax.set_yticks(range(cm.shape[0]))\n",
    "\n",
    "# Plot cm\n",
    "confmat(ax, cm)\n",
    "ax.set_title(\"Quick Check on Training Data\", fontsize=10)\n",
    "ax.set_yticklabels([archetype_decodedict[i] for i in range(1, cm.shape[0]+1)], \n",
    "                    rotation=45, va='top', fontsize=8)\n",
    "ax.set_xticklabels([archetype_decodedict[i] for i in range(1, cm.shape[0]+1)], \n",
    "                    rotation=45, ha='right', fontsize=8)\n",
    "ax.set_ylabel(\"Ground Truth\", labelpad=10)\n",
    "ax.set_xlabel(\"Prediction\", labelpad=10)\n",
    "\n",
    "# Done\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Writing the prediction results\n",
    "#   Note: This writes the results as numpy arrays associated with individual\n",
    "#         samples. The original results provided as a tsv in the IDR are not \n",
    "#         overwritten by this.\n",
    "\n",
    "# Report\n",
    "print \"Saving output...\"\n",
    "\n",
    "# Find samples in the image_data folder\n",
    "rawpath = r'data\\experimentA\\image_data'\n",
    "rawloader = ld.DataLoader(rawpath, recurse=True)\n",
    "\n",
    "# Load metadata\n",
    "meta, prim_IDs, _ = rawloader.load_dataset(\"stack_metadata.pkl\")\n",
    "\n",
    "# For each prim...\n",
    "for prim_ID in prim_IDs:\n",
    "    \n",
    "    # Find the segmentation path...\n",
    "    for fpath in rawloader.data[prim_ID]:\n",
    "        if fpath.endswith('_seg.tif'):\n",
    "    \n",
    "            # Construct prediction output path\n",
    "            outpath = fpath[:-4] + '_archetypes_' + fspace_type + '-PREDICTED.npy'\n",
    "            \n",
    "            # Save prediction data\n",
    "            np.save(outpath, archetypes_pred[fspace_idx==prim_IDs.index(prim_ID)])\n",
    "            \n",
    "            # Construct probabilities output path\n",
    "            outpath = fpath[:-4] + '_archetypes_' + fspace_type + '-PROBA.npy'\n",
    "            \n",
    "            # Save probability data\n",
    "            np.save(outpath, archetypes_prob[fspace_idx==prim_IDs.index(prim_ID)])\n",
    "            \n",
    "    # Add proba classes to metadata (just in case)\n",
    "    meta[prim_ID]['archetype_probaclasses'] = svc.classes_\n",
    "\n",
    "    # Save metadata\n",
    "    for fpath in rawloader.data[prim_ID]:\n",
    "        if '_stack_metadata.pkl' in fpath:\n",
    "            with open(fpath, 'wb') as metafile:\n",
    "                pickle.dump(meta[prim_ID], metafile, pickle.HIGHEST_PROTOCOL)\n",
    "\n",
    "# Save the fitted SVC separately for predictions on other datasets\n",
    "with open('other/archetype_svc_'+fspace_type+'.pkl', 'wb') as outfile:\n",
    "    pickle.dump(svc, outfile, pickle.HIGHEST_PROTOCOL)\n",
    "\n",
    "# Report\n",
    "print \"Done!\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}