RUN_SpotDetection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RUN: smFISH Spot Detection\n",
    "\n",
    "This notebook implements the spot detection for the pea3 smFISH data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes\n",
    "\n",
    "Spots are detected and counted using `skimage.blob.blob_log` from scikit-image. However, since my production version of scikit-image (`0.13.0`) does not support 3D blob detection, I copied and adapted the newer version (`0.15.dev0`) from github (see `katachi\\utilities\\skimage_blob.py`).\n",
    "\n",
    "Note that this approach is quite sensitive to image quality and parameters, so it will need to be adjusted for different smFISH experiments!"
   ]
  },
  {
   "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",
    "# smFISH\n",
    "import katachi.utilities.loading as ld\n",
    "import katachi.utilities.skimage_blob as blob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load data\n",
    "\n",
    "# Prep loader\n",
    "dirpath = r'data\\experimentB\\image_data'\n",
    "loader = ld.DataLoader(dirpath, recurse=True)\n",
    "\n",
    "# Load smFISH image data\n",
    "smf_stacks, prim_IDs, _ = loader.load_dataset(r\"pea3smFISH.tif\")\n",
    "\n",
    "# Load segmentations\n",
    "seg_stacks, _, _ = loader.load_dataset(r\"lynEGFP_seg.tif\", IDs=prim_IDs)\n",
    "\n",
    "# Get number of cells\n",
    "n_cells = [len(np.unique(seg_stacks[prim_ID]))-1 for prim_ID in prim_IDs]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some Testing"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Parameter screening (parallel)\n",
    "\n",
    "# WARNING: For some reason (leak, improper memory management on my part, or just inherently), this consumes way\n",
    "#          too much memory and pretty swiftly crashes the machine. Do not use as is!\n",
    "\n",
    "# Screen params\n",
    "n_test_imgs = 5\n",
    "n_processes = 5\n",
    "test_spec   = (0.15, 0.25, 20)\n",
    "\n",
    "# Generate required lists\n",
    "images_list = [smf_stacks[prim_ID] for prim_ID in prim_IDs[:n_test_imgs] * test_spec[2]]\n",
    "test_list   = np.tile(np.linspace(*test_spec), 8)\n",
    "\n",
    "# Assemble parameter dict\n",
    "param_dict = {'min_sigma' : 1,\n",
    "              'max_sigma' : 4,\n",
    "              'num_sigma' : 10,\n",
    "              'threshold' : test_list,\n",
    "              'overlap'   : 0.5,\n",
    "              'log_scale' : False}\n",
    "\n",
    "# Run the screen\n",
    "print \"Running...\"\n",
    "screen_results = blob.blob_log_multi(images_list, n_processes, param_dict)\n",
    "print \"Done!\""
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Parameter screening (serial)\n",
    "\n",
    "# WARNING: This takes quite a while...\n",
    "\n",
    "# Screen params\n",
    "n_test_imgs = 10\n",
    "#test_spec   = (0.15, 0.25, 20)  # threshold\n",
    "test_spec   = (0.50, 1.50, 20)  # min_sigma\n",
    "\n",
    "# Assemble parameter dict\n",
    "#param_dict = {'min_sigma' : 1,\n",
    "#              'max_sigma' : 4,\n",
    "#              'num_sigma' : 10,\n",
    "#              'overlap'   : 0.5,\n",
    "#              'log_scale' : False}\n",
    "param_dict = {'threshold' : 0.1,\n",
    "              'max_sigma' : 4,\n",
    "              'num_sigma' : 10,\n",
    "              'overlap'   : 0.5,\n",
    "              'log_scale' : False}\n",
    "\n",
    "# Run the screen\n",
    "print \"Running...\"\n",
    "screen_results = np.empty((test_spec[2], n_test_imgs))\n",
    "for t, test in enumerate(np.linspace(*test_spec)):\n",
    "    print \"  Started test value\", t+1, \"of\", test_spec[2]\n",
    "    for n, prim_ID in enumerate(prim_IDs[:n_test_imgs]):\n",
    "        #result = blob.blob_log(smf_stacks[prim_ID][:,-250:,:350], threshold=test, **param_dict)  # threshold\n",
    "        result = blob.blob_log(smf_stacks[prim_ID][:,-250:,:350], min_sigma=test, **param_dict)  # min_sigma\n",
    "        screen_results[t,n] = result.shape[0]\n",
    "        print \"    Completed prim\", n+1, \"of\", n_test_imgs\n",
    "    print \"  Finished test value\", t+1, \"of\", test_spec[2]\n",
    "print \"Done!\""
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Save screening results\n",
    "with open('other/smFISH_screen_threshold.pkl','wb') as outfile:\n",
    "    pickle.dump(screen_results, outfile, pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Save screening results\n",
    "with open('other/smFISH_screen_minsigma.pkl.pkl','wb') as outfile:\n",
    "    pickle.dump(screen_results, outfile, pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Load screening results\n",
    "with open('other/smFISH_screen_threshold.pkl','rb') as infile:\n",
    "    screen_results = pickle.load(infile)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Show screening results\n",
    "\n",
    "for i in range(test_spec[2]):\n",
    "    plt.scatter(n_cells[:n_test_imgs], screen_results[i,:], \n",
    "                c=[i for j in range(n_test_imgs)], vmin=0, vmax=test_spec[2])\n",
    "plt.show()\n",
    "\n",
    "for i in range(n_test_imgs):\n",
    "    plt.plot(np.linspace(*test_spec), screen_results[:,i])\n",
    "plt.show()\n",
    "\n",
    "plt.plot(np.linspace(*test_spec), screen_results.mean(axis=1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running Spot Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "### Run spot detection\n",
    "\n",
    "# Relevant parameters [see help(blob_log) for details]\n",
    "min_sigma      = 1     # Minimum standard deviation for Gaussian kernel\n",
    "max_sigma      = 4     # Maximum standard deviation for Gaussian kernel\n",
    "num_sigma      = 10    # Number of standard deviation steps to run\n",
    "spot_threshold = 0.22  # Determines which spots to keep (lower is more inclusive, higher is more stringent)\n",
    "overlap        = 0.50  # Smaller blob eliminated when area of two blobs overlaps by this much (separation resolution)\n",
    "\n",
    "# For each prim...\n",
    "print \"Running spot detection...\"\n",
    "for prim_ID in prim_IDs:\n",
    "    \n",
    "    # Report\n",
    "    print \"  Working on prim\", prim_ID\n",
    "    \n",
    "    # Run spot detection\n",
    "    blobs = blob.blob_log(smf_stacks[prim_ID], \n",
    "                          min_sigma=min_sigma, max_sigma=max_sigma, num_sigma=num_sigma, \n",
    "                          threshold=spot_threshold, overlap=overlap, log_scale=False)\n",
    "    \n",
    "    # Get corresponding cell label for each spot\n",
    "    blob_coords = blobs[:, :-1].astype(np.int)\n",
    "    hit_labels  = seg_stacks[prim_ID][blob_coords[:,0], \n",
    "                                      blob_coords[:,1], \n",
    "                                      blob_coords[:,2]]\n",
    "    \n",
    "    # Get counts\n",
    "    cell_counts = np.zeros(np.unique(seg_stacks[prim_ID]).size, dtype=np.int)\n",
    "    hit_labels_id, hit_counts  = np.unique(hit_labels, return_counts=True)\n",
    "    cell_counts[hit_labels_id] = hit_counts\n",
    "    cell_counts = cell_counts[1:]\n",
    "    \n",
    "    # Save spots file\n",
    "    base_path  = [p for p in loader.data[prim_ID] \n",
    "                  if p.endswith(r\"pea3smFISH.tif\")][0]\n",
    "    spots_path = base_path[:-4] + '_RNAspots.npy'\n",
    "    np.save(spots_path, blobs)\n",
    "    \n",
    "    # Save counts file\n",
    "    counts_path = base_path[:-4] + '_RNAcounts.npy'\n",
    "    np.save(counts_path, cell_counts)\n",
    "    \n",
    "    # Report\n",
    "    print \"    Detected\", blobs.shape[0], \"spots!\"\n",
    "    \n",
    "# Report\n",
    "print \"Processing complete!\""
   ]
  }
 ],
 "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
}