examples/Basic usage.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Object-oriented FISSA interface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains a step-by-step example of how to use the object-oriented (class-based) interface to the [FISSA](https://github.com/rochefort-lab/fissa) toolbox.\n",
    "\n",
    "For more details about the methodology behind FISSA, please see our paper:\n",
    "> S. W. Keemink, S. C. Lowe, J. M. P. Pakan, E. Dylda, M. C. W. van Rossum, and N. L. Rochefort. FISSA: A neuropil decontamination toolbox for calcium imaging signals, *Scientific Reports*, **8**(1):3493, 2018. doi: [10.1038/s41598-018-21640-2](https://www.doi.org/10.1038/s41598-018-21640-2).\n",
    "\n",
    "See [basic_usage.py](https://github.com/rochefort-lab/fissa/blob/master/examples/basic_usage.py) (or [basic_usage_windows.py](https://github.com/rochefort-lab/fissa/blob/master/examples/basic_usage_windows.py) for Windows users) for a short example script outside of a notebook interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the FISSA toolbox\n",
    "import fissa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will be plotting our results using [HoloViews](http://holoviews.org/), but you can use any plotting library you prefer ([matplotlib](https://matplotlib.org/), [bokeh](https://docs.bokeh.org/), etc)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import our plotting toolbox, and enable options for embedded notebook figures\n",
    "import holoviews as hv\n",
    "\n",
    "%load_ext holoviews.ipython\n",
    "%output widgets=\"embed\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining an experiment\n",
    "\n",
    "To run a separation step with fissa, you need create a `fissa.Experiment` object, which will hold your extraction parameters and results.\n",
    "\n",
    "The inputs to the `fissa.Experiment` instance are the:\n",
    "- experiment images\n",
    "- the regions of interest (ROIs) to extract\n",
    "Define your inputs. All that's necessary to define are the image data and ROIs. \n",
    "\n",
    "Images can be defined as a folder with tiff stacks:\n",
    "```python\n",
    "images = \"folder\"\n",
    "```\n",
    "Where each tiff stack in the folder is a trial with several frames.\n",
    "\n",
    "Or the data can also be given as a list of arrays if not stored as tiffs:\n",
    "```python\n",
    "images = [array1, array2, array3, ...]\n",
    "```\n",
    "\n",
    "For ROIs either a set of ROIs across all images should be defined, or a set of ROIs per image. \n",
    "\n",
    "If the ROIs were defined using ImageJ use ImageJ's export function to save them in a zip. Then, indicate the locations as a list:\n",
    "```python\n",
    "rois = \"rois.zip\"  # for a single set of rois across images\n",
    "rois = [\"rois1.zip\", \"rois2.zip\", ...]  # for a roiset for each image\n",
    "```\n",
    "Defining a different roiset per image can be useful if you need to adjust for motion drift for example.\n",
    "\n",
    "Then, we can define our experiment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the folder where FISSA's outputs will be cached, so they can be\n",
    "# quickly reloaded in the future without having to recompute them.\n",
    "#\n",
    "# This argument is optional; if it is not provided, FISSA will not save its\n",
    "# results for later use.\n",
    "#\n",
    "# Note: you *must* use a different folder for each experiment.\n",
    "#\n",
    "# In this example, we will use the current datetime as the\n",
    "# name of the experiment, but you can name your experiments\n",
    "# however you want to.\n",
    "\n",
    "import datetime\n",
    "\n",
    "output_folder = \"fissa-example_{}\".format(\n",
    "    datetime.datetime.now().strftime(\"%Y-%m-%d_%H-%M-%S\")\n",
    ")\n",
    "\n",
    "print(output_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define image and ROI locations\n",
    "images_location = \"exampleData/20150529\"\n",
    "rois_location = \"exampleData/20150429.zip\"\n",
    "\n",
    "# Define the folder where FISSA's outputs will be stored, so they can be\n",
    "# quickly reloaded in the future without having to recompute them.\n",
    "# Make sure to use a different folder for each experiment.\n",
    "# This argument is optional; if it is not provided, FISSA will not save its\n",
    "# results for later use.\n",
    "output_folder = \"fissa_example\"\n",
    "\n",
    "experiment = fissa.Experiment(images_location, rois_location, output_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previously analyzed experiments in output_folder will be loaded, if they exist, and the next step could be skipped."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting traces and separating them\n",
    "Next, we need to extract the traces and separate them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment.separate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to redo preparation and/or separation you can set:\n",
    "```\n",
    "experiment.separate(redo_prep=True, redo_sep=True)\n",
    "```\n",
    "(If you redo prepartion this will also redo the separation, to make sure these always match up)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing results\n",
    "\n",
    "After running ```experiment.separate()``` the results are stored as follows.\n",
    "\n",
    "### ROI outlines\n",
    "\n",
    "The ROI outlines, as well as the extra neuropil regions, can be found as in ```experiment.roi_polys``` as follows. For cell number ```c``` and tiff number `t`, the set of ROIs for that cell and tiff is at\n",
    "```python\n",
    "experiment.roi_polys[c][t][0][0]  # basic ROI\n",
    "experiment.roi_polys[c][t][n][0]  # n = 1, 2, 3, .... the neuropil regions\n",
    "```\n",
    "Sometimes ROIs cannot be expressed as a single polygon (e.g. a ring-ROI), in those cases several polygons are used to describe it as:\n",
    "```python\n",
    "experiment.roi_polys[c][t][n][i]  # i iterates over the different polygons\n",
    "```\n",
    "\n",
    "As an example, plotting the first region of interest plus its surrounding neuropil subregions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualise the ROI polygon for a sample cell\n",
    "c = 0\n",
    "t = 0\n",
    "\n",
    "neuropil1 = hv.Curve(experiment.roi_polys[c][t][1][0])\n",
    "neuropil2 = hv.Curve(experiment.roi_polys[c][t][2][0])\n",
    "neuropil3 = hv.Curve(experiment.roi_polys[c][t][3][0])\n",
    "neuropil4 = hv.Curve(experiment.roi_polys[c][t][4][0])\n",
    "cell = hv.Curve(experiment.roi_polys[c][t][0][0])\n",
    "\n",
    "neuropil1 * neuropil2 * neuropil3 * neuropil4 * cell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, the ROI (grey) is on the boundary of the image, so the surrounding neuropil subregions are arranged to share the available space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FISSA extracted traces\n",
    "\n",
    "The final extracted traces can be found in ```experiment.result``` as follows. For cell number ```c``` and tiff number `t`, the final extracted trace is given by:\n",
    "```python\n",
    "experiment.result[c][t][0, :]\n",
    "```\n",
    "\n",
    "In `experiment.result` one can find the signals present in the cell ROI, ordered by how strongly they are present (relative to the surrounding regions). `experiment.result[c][t][0, :]` gives the most strongly present signal, and is considered the cell's \"true\" signal. `[i, :]` for `i=1,2,3,...` gives the other signals which are present in the ROI, but driven by other cells or neuropil.\n",
    "\n",
    "### Before decontamination\n",
    "\n",
    "The raw extracted signals can be found in `experiment.raw` in the same way. Now in  `experiment.raw[c][t][i,:]`, `i` indicates the region number, with `i=0` being the cell, and `i=1,2,3,...` indicating the surrounding regions.\n",
    "\n",
    "As an example, plotting the raw and extracted signals for the second trial for the third cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 2\n",
    "t = 1\n",
    "\n",
    "(\n",
    "    hv.Curve(experiment.raw[c][t][0, :], label=\"Raw\")\n",
    "    * hv.Curve(experiment.result[c, t][0, :], label=\"Decontaminated\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### df/f<sub>0</sub>\n",
    "\n",
    "It is often useful to calculate the intensity of a signal relative to the baseline value, df/f<sub>0</sub>, for the traces.\n",
    "This can be done as follows.\n",
    "\n",
    "Note that by default, f<sub>0</sub> is determined as the minimum across all trials (all tiffs) to ensure that results are directly comparable between trials, but you can normalise each trial individually instead if you prefer by setting `across_trials=False`.\n",
    "\n",
    "Since FISSA is very good at removing contamination from the ROI signals, the minimum value on the decontaminated trace will typically be `0.`. Consequently, we use the minimum value of the (smoothed) raw signal to provide the f<sub>0</sub> from the raw trace for both the raw and decontaminated df/f<sub>0</sub>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment.calc_deltaf(freq=10, across_trials=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 2\n",
    "t = 1\n",
    "\n",
    "(\n",
    "    hv.Curve(experiment.deltaf_raw[c][t], label=\"Raw\", vdims=[\"df/f0\"])\n",
    "    * hv.Curve(experiment.deltaf_result[c, t][0, :], label=\"Decontaminated\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exporting to MATLAB\n",
    "\n",
    "The results can easily be exported to a MATLAB-compatible matfile as follows.\n",
    "\n",
    "The output will appear in the `output_folder` we supplied to `experiment` when we created it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment.save_to_matlab()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loading `output_folder/matlab.mat` in MATLAB will give you three structs, `ROIs`, `raw`, and `result`.\n",
    "\n",
    "These interface similarly as `experiment.ROIs`, `experiment.raw`, and `experiment.result` described above. However, Matlab counts from 1 (as opposed to Python counting from 0), such that the ROI, raw trace, and decontaminated trace are all found for cell 0 trial 0 as:\n",
    "\n",
    "```octave\n",
    "ROIs.cell0.trial0{1}  % polygon for the ROI\n",
    "ROIs.cell0.trial0{2}  % polygon for first neuropil subregion\n",
    "result.cell0.trial0(1, :)  % final extracted cell signal\n",
    "result.cell0.trial0(2, :)  % contaminating signal\n",
    "raw.cell0.trial0(1, :)  % raw measured celll signal\n",
    "raw.cell0.trial0(2, :)  % raw signal from first neuropil subregion\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Addendum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the tiff files\n",
    "\n",
    "If you find something noteworthy in one of the traces and need to backreference to the corresponding tiff file, you can look up the path to the tiff file with `experiment.images`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trial_of_interest = 1\n",
    "\n",
    "print(experiment.images[trial_of_interest])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean image data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get the temporal-mean image for a trial with `experiment.means`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Image(experiment.means[trial_of_interest])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also superimpose the cell ROIs on the temporal-mean image as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using holoviews\n",
    "\n",
    "t = trial_of_interest\n",
    "\n",
    "fig = hv.Raster(experiment.means[t])\n",
    "\n",
    "for c in range(experiment.nCell):\n",
    "    roi_poly = experiment.roi_polys[c][t][0][0]\n",
    "    x = roi_poly[:, 1]\n",
    "    y = roi_poly[:, 0]\n",
    "    fig *= hv.Curve(zip(x, y))\n",
    "\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using matplotlib\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "t = trial_of_interest\n",
    "\n",
    "plt.imshow(experiment.means[t], cmap=\"gray\")\n",
    "for roi_poly in experiment.roi_polys:\n",
    "    # Plot border around cells\n",
    "    # plt.plot(roi_poly[t][0][0][:, 1], roi_poly[t][0][0][:, 0], \":c\")\n",
    "    # Fill cells with partially-transparent shaded area\n",
    "    plt.fill(roi_poly[t][0][0][:, 1], roi_poly[t][0][0][:, 0], \":r\", alpha=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FISSA customisation settings\n",
    "\n",
    "FISSA has several user-definable settings, which can be set when defining the `experiment` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FISSA uses multiprocessing to speed up its processing.\n",
    "# By default, it will spawn one worker per CPU core on your machine.\n",
    "# However, if you have a lot of cores and not much memory, you many not\n",
    "# be able to suport so many workers simultaneously.\n",
    "# In particular, this can be problematic during the data preparation step\n",
    "# in which tiffs are loaded into memory.\n",
    "# The default number of cores for the data preparation and separation steps\n",
    "# can be changed as follows.\n",
    "ncores_preparation = 4  # If None, uses all available cores\n",
    "ncores_separation = None  # if None, uses all available cores\n",
    "\n",
    "# By default, FISSA uses 4 subregions for the neuropil region.\n",
    "# If you have very dense data with a lot of different signals per unit area,\n",
    "# you may wish to increase the number of regions.\n",
    "nRegions = 8\n",
    "\n",
    "# By default, each surrounding region has the same area as the central ROI.\n",
    "# i.e. expansion = 1\n",
    "# However, you may wish to increase or decrease this value.\n",
    "expansion = 0.75\n",
    "\n",
    "# The degree of signal sparsity can be controlled with the alpha parameter.\n",
    "alpha = 0.1\n",
    "\n",
    "# Set up a FISSA experiment with these parameters\n",
    "experiment = fissa.Experiment(\n",
    "    images_location,\n",
    "    rois_location,\n",
    "    output_folder,\n",
    "    nRegions=nRegions,\n",
    "    expansion=expansion,\n",
    "    alpha=alpha,\n",
    "    ncores_preparation=ncores_preparation,\n",
    "    ncores_separation=ncores_separation,\n",
    ")\n",
    "# Extract the data with these new parameters.\n",
    "# Note that we are using the same output folder as before. Since FISSA has cached\n",
    "# a result alrady to this directory, its default behaviour is to restore the\n",
    "# previously generated results.\n",
    "# To make sure FISSA runs a fresh decontamination process with the new parameters,\n",
    "# we need to make sure to specify to redo the preparation and separation.\n",
    "# FISSA will then ignore the cached output and overwrite it with new results.\n",
    "experiment.separate(redo_prep=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the new results for our example trace from before. Although we doubled the number of neuropil regions around the cell, very little has changed for this example because there were not many sources of contamination.\n",
    "\n",
    "However, there will be more of a difference if your data has more neuropil sources per unit area within the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the new results\n",
    "c = 2\n",
    "t = 1\n",
    "\n",
    "(\n",
    "    hv.Curve(experiment.raw[c][t][0, :], label=\"Raw\")\n",
    "    * hv.Curve(experiment.result[c, t][0, :], label=\"Decontaminated\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, these settings can be refined after creating the `experiment` object, as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment.ncores_preparation = 8\n",
    "experiment.alpha = 0.1\n",
    "experiment.expansion = 0.75"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading data from large tiff files\n",
    "\n",
    "By default, FISSA loads entire tiff files into memory at once and then manipulates all ROIs within the tiff.\n",
    "This can sometimes be problematic when working with very large tiff files which can not be loaded into memory all at once.\n",
    "If you have out-of-memory problems, you can activate FISSA's low memory mode, which will cause it to manipulate each tiff file frame-by-frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment = fissa.Experiment(\n",
    "    images_location, rois_location, output_folder, lowmemory_mode=True\n",
    ")\n",
    "experiment.separate(redo_prep=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Handling custom formats\n",
    "\n",
    "By default, FISSA can use tiff files or numpy arrays as its input image data, and numpy arrays or ImageJ zip files for the ROI definitions.\n",
    "However, it is also possible to extend this functionality and integrate other data formats into FISSA in order to work with other custom and/or proprietary formats that might be used in your lab.\n",
    "\n",
    "This is done by defining your own DataHandler class.\n",
    "Your custom data handler should be a subclass of [fissa.extraction.DataHandlerAbstract](https://fissa.readthedocs.io/en/stable/source/packages/fissa.extraction.html#fissa.extraction.DataHandlerAbstract), and implement the following methods:\n",
    "\n",
    "- `image2array(image)` takes an image of whatever format and turns it into *data* (typically a numpy.ndarray).\n",
    "- `getmean(data)` calculates the 2D mean for a video defined by *data*.\n",
    "- `rois2masks(rois, data)` creates masks from the rois inputs, of appropriate size *data*.\n",
    "- `extracttraces(data, masks)` applies the *masks* to *data* in order to extract traces.\n",
    "\n",
    "See [fissa.extraction.DataHandlerAbstract](https://fissa.readthedocs.io/en/stable/source/packages/fissa.extraction.html#fissa.extraction.DataHandlerAbstract) for further description for each of the methods.\n",
    "\n",
    "If you only need to handle a new *image* input format, which is converted to a numpy.ndarray, you may find it is easier to create a subclass of the default datahandler, [fissa.extraction.DataHandlerTifffile](https://fissa.readthedocs.io/en/stable/source/packages/fissa.extraction.html#fissa.extraction.DataHandlerTifffile).\n",
    "In this case, only the `image2array` method needs to be overwritten and the other methods can be left as they are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fissa.extraction import DataHandlerTifffile\n",
    "\n",
    "# Define a custom datahandler class.\n",
    "#\n",
    "# By inheriting from DataHandlerTifffile, most methods are defined\n",
    "# appropriately. In this case, we only need to overwrite the\n",
    "# `image2array` method to work with our custom data format.\n",
    "\n",
    "\n",
    "class DataHandlerCustom(DataHandlerTifffile):\n",
    "    @staticmethod\n",
    "    def image2array(image):\n",
    "        \"\"\"Open a given image file as a custom instance.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        image : custom\n",
    "            Your image format (avi, hdf5, etc.)\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        numpy.ndarray\n",
    "            A 3D array containing the data, shaped\n",
    "            ``(frames, y_coordinate, x_coordinate)``.\n",
    "        \"\"\"\n",
    "        # Some custom code\n",
    "        pass\n",
    "\n",
    "\n",
    "# Then pass an instance of this class to fissa.Experiment when creating\n",
    "# a new experiment.\n",
    "datahandler = DataHandlerCustom()\n",
    "\n",
    "experiment = fissa.Experiment(\n",
    "    images_location,\n",
    "    rois_location,\n",
    "    datahandler=datahandler,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For advanced users that want to entirely replace the DataHandler with their own methods, you can also inherit a class from the abstract class, [fissa.extraction.DataHandlerAbstract](https://fissa.readthedocs.io/en/stable/source/packages/fissa.extraction.html#fissa.extraction.DataHandlerAbstract).\n",
    "This can be useful if you want to integrate FISSA into your workflow without changing everything into the numpy array formats that FISSA usually uses internally."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}