examples/Suite2p example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using FISSA with Suite2p\n",
    "\n",
    "[suite2p](https://github.com/MouseLand/suite2p) is a blind source separation toolbox for cell detection and signal extraction. \n",
    "\n",
    "Here we illustrate how to use suite2p to detect cell locations, and then use FISSA to remove neuropil signals from the ROI signals.\n",
    "\n",
    "The suite2p parts of this tutorial are based on their [Jupyter notebook example](https://github.com/MouseLand/suite2p/blob/master/jupyter/run_pipeline_tiffs_or_batch.ipynb).\n",
    "\n",
    "Note that the below results are not representative of either suite2p or FISSA performance, as we are using a very small example dataset.\n",
    "\n",
    "**Reference:**\n",
    "> Pachitariu, M., Stringer, C., Dipoppa, M., Schröder, S., Rossi, L. F., Dalgleish, H., Carandini, M. & Harris, K. D. (2017). Suite2p: beyond 10,000 neurons with standard two-photon microscopy. bioRxiv:&nbsp;[061507](https://www.biorxiv.org/content/10.1101/061507v2); doi:&nbsp;[10.1101/061507](https://doi.org/10.1101/061507)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FISSA toolbox\n",
    "import fissa\n",
    "\n",
    "# suite2p toolbox\n",
    "import suite2p.run_s2p\n",
    "\n",
    "# For plotting our results, use numpy and matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run suite2p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set your options for running\n",
    "ops = suite2p.run_s2p.default_ops()  # populates ops with the default options\n",
    "\n",
    "# Provide an h5 path in 'h5py' or a tiff path in 'data_path'\n",
    "# db overwrites any ops (allows for experiment specific settings)\n",
    "db = {\n",
    "    \"h5py\": [],  # a single h5 file path\n",
    "    \"h5py_key\": \"data\",\n",
    "    \"look_one_level_down\": False,  # whether to look in ALL subfolders when searching for tiffs\n",
    "    \"data_path\": [\"exampleData/20150529\"],  # a list of folders with tiffs\n",
    "    # (or folder of folders with tiffs if look_one_level_down is True,\n",
    "    # or subfolders is not empty)\n",
    "    \"save_path0\": \"./\",  # save path\n",
    "    \"subfolders\": [],  # choose subfolders of 'data_path' to look in (optional)\n",
    "    \"fast_disk\": \"./\",  # string which specifies where the binary file will be stored (should be an SSD)\n",
    "    \"reg_tif\": True,  # save the motion corrected tiffs\n",
    "    \"tau\": 0.7,  # timescale of gcamp6f\n",
    "    \"fs\": 1,  # sampling rate\n",
    "    \"spatial_scale\": 4,\n",
    "}\n",
    "\n",
    "# Run one experiment\n",
    "opsEnd = suite2p.run_s2p.run_s2p(ops=ops, db=db)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the relevant data from the analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the motion corrected tiffs (make sure that the reg_tif option is set to true, see above)\n",
    "images = \"./suite2p/plane0/reg_tif\"\n",
    "\n",
    "# Load the detected regions of interest\n",
    "stat = np.load(\"./suite2p/plane0/stat.npy\", allow_pickle=True)  # cell stats\n",
    "ops = np.load(\"./suite2p/plane0/ops.npy\", allow_pickle=True).item()\n",
    "iscell = np.load(\"./suite2p/plane0/iscell.npy\", allow_pickle=True)[:, 0]\n",
    "\n",
    "# Get image size\n",
    "Lx = ops[\"Lx\"]\n",
    "Ly = ops[\"Ly\"]\n",
    "\n",
    "# Get the cell ids\n",
    "ncells = len(stat)\n",
    "cell_ids = np.arange(ncells)  # assign each cell an ID, starting from 0.\n",
    "cell_ids = cell_ids[iscell == 1]  # only take the ROIs that are actually cells.\n",
    "num_rois = len(cell_ids)\n",
    "\n",
    "# Generate ROI masks in a format usable by FISSA (in this case, a list of masks)\n",
    "rois = [np.zeros((Ly, Lx), dtype=bool) for n in range(num_rois)]\n",
    "\n",
    "for i, n in enumerate(cell_ids):\n",
    "    # i is the position in cell_ids, and n is the actual cell number\n",
    "    ypix = stat[n][\"ypix\"][~stat[n][\"overlap\"]]\n",
    "    xpix = stat[n][\"xpix\"][~stat[n][\"overlap\"]]\n",
    "    rois[i][ypix, xpix] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run FISSA with the defined ROIs and data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_folder = \"fissa_suite2p_example\"\n",
    "exp = fissa.Experiment(images, [rois[:ncells]], output_folder)\n",
    "exp.separate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the resulting ROI signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fetch the colormap object for Cynthia Brewer's Paired color scheme\n",
    "cmap = plt.get_cmap(\"Paired\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select which trial (TIFF index) to plot\n",
    "trial = 0\n",
    "\n",
    "# Plot the mean image and ROIs from the FISSA experiment\n",
    "plt.imshow(experiment.means[trial], cmap=\"gray\")\n",
    "\n",
    "XLIM = plt.xlim()\n",
    "YLIM = plt.ylim()\n",
    "\n",
    "for i_roi, roi_poly in enumerate(experiment.roi_polys):\n",
    "    # Plot border around ROI\n",
    "    plt.plot(\n",
    "        roi_poly[trial][0][0][:, 1],\n",
    "        roi_poly[trial][0][0][:, 0],\n",
    "        color=cmap((i_roi * 2 + 1) % cmap.N),\n",
    "    )\n",
    "\n",
    "# ROI co-ordinates are half a pixel outside the image,\n",
    "# so we reset the axis limits\n",
    "plt.xlim(XLIM)\n",
    "plt.ylim(YLIM)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all ROIs and trials\n",
    "\n",
    "# Get the number of ROIs and trials\n",
    "n_roi = experiment.result.shape[0]\n",
    "n_trial = experiment.result.shape[1]\n",
    "\n",
    "# Find the maximum signal intensities for each ROI\n",
    "roi_max_raw = [\n",
    "    np.max([np.max(experiment.raw[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max_result = [\n",
    "    np.max([np.max(experiment.result[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max = np.maximum(roi_max_raw, roi_max_result)\n",
    "\n",
    "# Plot our figure using subplot panels\n",
    "plt.figure(figsize=(16, 2.5 * n_roi))\n",
    "for i_roi in range(n_roi):\n",
    "    for i_trial in range(n_trial):\n",
    "        # Make subplot axes\n",
    "        i_subplot = 1 + i_roi * n_trial + i_trial\n",
    "        plt.subplot(n_roi, n_trial, i_subplot)\n",
    "        # Plot the data\n",
    "        plt.plot(\n",
    "            experiment.raw[i_roi][i_trial][0, :],\n",
    "            label=\"Raw (suite2p)\",\n",
    "            color=cmap(i_roi * 2 % cmap.N),\n",
    "        )\n",
    "        plt.plot(\n",
    "            experiment.result[i_roi][i_trial][0, :],\n",
    "            label=\"FISSA\",\n",
    "            color=cmap((i_roi * 2 + 1) % cmap.N),\n",
    "        )\n",
    "        # Labels and boiler plate\n",
    "        plt.ylim([-0.05 * roi_max[i_roi], roi_max[i_roi] * 1.05])\n",
    "        if i_trial == 0:\n",
    "            plt.ylabel(\"Signal intensity\\n(candela per unit area)\")\n",
    "        if i_roi == 0:\n",
    "            plt.title(\"Trial {}\".format(i_trial))\n",
    "        if i_roi == n_roi - 1:\n",
    "            plt.xlabel(\"Time (frame number)\")\n",
    "            plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figure shows the raw signal from the ROI identified by suite2p (pale), and after decontaminating with FISSA (dark).\n",
    "The hues match the ROI locations drawn above.\n",
    "Each row shows the results from one of the ROIs detected by suite2p.\n",
    "Each column shows the results from one of the three trials.\n",
    "\n",
    "Note that with the above settings for suite2p it seems to have detected more small local axon signals, instead of cells. This can possibly be improved with manual curation and suite2p setting changes, but as noted above these results should not be seen as indicative for either suite2p or FISSA due to the small dataset size.\n",
    "\n",
    "Also note that the above Suite2P traces are done without suite2p's own neuropil removal algorithm."
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}