spot_detection_tutorial_beginner.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spot detection with napari\n",
    "\n",
    "\n",
    "### Overview\n",
    "In this activity, we will perform spot detection on some in situ sequencing data ([Feldman and Singh et al., Cell, 2019](https://www.cell.com/cell/fulltext/S0092-8674(19)31067-0s)). In doing so, we will combine methods from [scipy](https://www.scipy.org/), [scikit-image](https://scikit-image.org/), and [cellpose](https://github.com/MouseLand/cellpose). The goal is to familiarize you with performing analysis that integrates the scientific python ecosystem and napari.\n",
    "\n",
    "### Data source\n",
    "\n",
    "The data were downloaded from the [OpticalPooledScreens github repository](https://github.com/feldman4/OpticalPooledScreens).\n",
    "\n",
    "### Next steps\n",
    "\n",
    "Following this activity, we will use the workflow generated in this activity to create a napari spot detection plugin. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load the data\n",
    "\n",
    "In the cells below load the data using the scikit-image `imread()` function. For more information about the `imread()` function, please see the [scikit-image docs](https://scikit-image.org/docs/dev/api/skimage.io.html#skimage.io.imread). We are loading two images:\n",
    "\n",
    "- `nuclei`: an image of cell nuclei\n",
    "- `spots`: an image of in situ sequencing spots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import io\n",
    "\n",
    "nuclei = io.imread('./data/nuclei_cropped.tif')\n",
    "spots = io.imread('./data/spots_cropped.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# View the data\n",
    "\n",
    "We will use napari to view our data. To do so, we first must create the viewer. Once the Viewer is created, we can add images to the viewer via the Viewer's `add_image()` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import napari\n",
    "\n",
    "# create the napari viewer\n",
    "viewer = napari.Viewer();\n",
    "\n",
    "# add the nuclei image to the viewer\n",
    "viewer.add_image(nuclei);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cell below, add the spots image to the viewer as was done above for the nuclei image. After loading the data, inspect it in the viewer and adjust the layer settings to your liking (e.g., contrast limits, colormap). You can pan/zoom around the image by click/dragging to pan and scrolling with your mousewheel or trackpad to zoom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add the spots image to the viewer\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create an image filter\n",
    "\n",
    "You may have noticed the the spots image contains background and autofluorescence from the cells. To improve spot detection, we will apply a high pass filter to improve the contrast of the spots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import ndimage as ndi\n",
    "\n",
    "def gaussian_high_pass(image: np.ndarray, sigma: float = 2):\n",
    "    \"\"\"Apply a gaussian high pass filter to an image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    image : np.ndarray\n",
    "        The image to be filtered.\n",
    "    sigma : float\n",
    "        The sigma (width) of the gaussian filter to be applied.\n",
    "        The default value is 2.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    high_passed_im : np.ndarray\n",
    "        The image with the high pass filter applied\n",
    "    \"\"\"\n",
    "    low_pass = ndi.gaussian_filter(image, sigma)\n",
    "    high_passed_im = image - low_pass\n",
    "    \n",
    "    return high_passed_im"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cell below, apply the gaussian high pass filter to the `spots` image and add the image to the viewer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the gaussian_high_pass function to filter the spots image\n",
    "\n",
    "\n",
    "# add the filtered image to the viewer\n",
    "# hint: set the opacity < 1 in order to see the layers underneath\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Detect spots\n",
    "\n",
    "Next, we will create a function to detect the spots in the spot image. This function should take the raw image, apply the gaussian high pass filter from above and then use one of the blob detection algorithms from sci-kit image to perform the blob detection. The `detect_spots()` function should return a numpy array containing the coordinates of each spot and a numpy array containing the diameter of each spot.\n",
    "\n",
    "Some hints:\n",
    "- See the [blob detection tutorial from scikit-image](https://scikit-image.org/docs/dev/auto_examples/features_detection/plot_blob.html). - We recommend the [blob_log detector](https://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.blob_log), but feel free to experiment!\n",
    "- See the \"Note\" from the blob_log docs: \"The radius of each blob is approximately $\\sqrt{2}\\sigma$ for a 2-D image\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from skimage.feature import blob_log\n",
    "\n",
    "def detect_spots(\n",
    "    image: np.ndarray,\n",
    "    high_pass_sigma: float = 2,\n",
    "    spot_threshold: float = 0.01,\n",
    "    blob_sigma: float = 2\n",
    "):\n",
    "    \"\"\"Apply a gaussian high pass filter to an image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    image : np.ndarray\n",
    "        The image in which to detect the spots.\n",
    "    high_pass_sigma : float\n",
    "        The sigma (width) of the gaussian filter to be applied.\n",
    "        The default value is 2.\n",
    "    spot_threshold : float\n",
    "        The threshold to be passed to the blob detector.\n",
    "        The default value is 0.01.\n",
    "    blob_sigma: float\n",
    "        The expected sigma (width) of the spots. This parameter\n",
    "        is passed to the \"max_sigma\" parameter of the blob\n",
    "        detector.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    points_coords : np.ndarray\n",
    "        An NxD array with the coordinate for each detected spot.\n",
    "        N is the number of spots and D is the number of dimensions.\n",
    "    sizes : np.ndarray\n",
    "        An array of size N, where N is the number of detected spots\n",
    "        with the diameter of each spot.\n",
    "    \n",
    "    \"\"\"\n",
    "    # filter the image with the gaussian_high_pass filter\n",
    "    filtered_spots = gaussian_high_pass(image, high_pass_sigma)\n",
    "\n",
    "    # detect the spots on the filtered image\n",
    "    blobs_log = blob_log(\n",
    "        filtered_spots,\n",
    "        max_sigma=blob_sigma,\n",
    "        num_sigma=1,\n",
    "        threshold=spot_threshold\n",
    "    )\n",
    "    \n",
    "    # convert the output of the blob detector to the \n",
    "    # desired points_coords and sizes arrays\n",
    "    # (see the docstring for details)\n",
    "    points_coords = blobs_log[:, 0:2]\n",
    "    sizes = 3 * blobs_log[:, 2]\n",
    "\n",
    "    return points_coords, sizes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cell below, apply `detect_spots()` to our `spots` image. To visualize the results, add the spots to the viewer as a [Points layer](https://napari.org/tutorials/fundamentals/points.html). If you would like to see an example of using a points layer, see [this example](https://github.com/napari/napari/blob/master/examples/add_points.py). To test out your function, vary the detection parameters and see how they affect the results. Note that each time you run the cell, the new results are added as an addition Points layer, allowing you to compare results from different parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# detect the spots\n",
    "\n",
    "\n",
    "# add the detected spots to the viewer as a Points layer\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Segment nuclei\n",
    "\n",
    "To segment the nuclei, we will use the [cellpose napari plugin](https://github.com/MouseLand/cellpose-napari). Through this and the following section, we will see how one can integrate plugins developed both others into their analysis workflow. Please perform the segmentation using the instructions below. For \n",
    "more information on cellpose, please see the [paper](https://www.nature.com/articles/s41592-020-01018-x) and [repository](https://github.com/MouseLand/cellpose).\n",
    "\n",
    "1. Start the cellpose plugin. From the menu bar, click Plugins->cellpose-napari: cellpose. You should see the plugin added to the right side of the viewer.\n",
    "<img src=\"./resources/cellpose_plugin.png\">\n",
    "2. Select the \"nuclei\" image layer. \n",
    "<img src=\"./resources/cellpose_screenshots_image_selection.png\">\n",
    "3. Set the model type to \"nuclei\"\n",
    "<img src=\"./resources/cellpose_screenshots_model_selection.png\">\n",
    "4. We need to give cellpose an estimate of the size of the nuclei so it can properly scale the data. We can do so using a napari Shapes layer. With the Shapes layer, we will outline some nuclei and then cellpose will use those annotations to estimate the size of the nuclei.\n",
    "    1. Click the \"add Shapes\" layer button in the viewer. This will create and select a new layer called \"Shapes\".\n",
    "    <img src=\"./resources/cellpose_screenshots_add_shape.png\">\n",
    "    2. Set the mode to \"Ellipse\" by clicking the button in the layer controls.\n",
    "    3. In the canvas, click and drag to add an ellipse that around a \"representative\" nucleus. For the purpose of this demo, this is enough, but for other data you may need to give more examples to make a better estimate of the cell diameter. If you need to pan/zoom while adding an ellipse, holding the spacebar will allow you to pan/zoom using your mouse (pan via click/drag, zoom by scrolling).\n",
    "    4. If you would like to edit or move an ellipse, you can switch to \"Select shapes\" mode in the viewer. Shapes can now be moved by clicking on them and then dragging. They can be resized by selecting them and then dragging the control points.\n",
    "    <img src=\"./resources/cellpose_screenshots_select_shape.png\">\n",
    "    5. Once you are happy with your annotations, you can click the \"compute diameter from shape layer\" button and you will see the \"diameter\" value populated. For this demo, the value is typically around 10 pixels.\n",
    "    <img src=\"./resources/cellpose_screenshots_diameter.png\">\n",
    "5. For this demo, we recommend de-selecting \"average 4 nets\"(potentially less accurate, but faster segmentation) and otherwise using the default settings. If you would like to learn more about the cellpose settings, please see the [cellpose plugin documentation](https://cellpose-napari.readthedocs.io/en/latest/settings.html).\n",
    "<img src=\"./resources/cellpose_screenshots_settings.png\">\n",
    "6. Now you are ready to run the segmentation! Click the \"run segmentation\" button. Segmentation for this demo typically takes ~1 minute.\n",
    "<img src=\"./resources/cellpose_screenshots_run.png\">\n",
    "7. When the segmentation is completed, you will see some new layers added to the layer list. Of particular interest is \"nuclei_cp_masks_000\", which contains our segmentation mask added as a Labels layer.\n",
    "<img src=\"./resources/cellpose_screenshots_results.png\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extension activity: assign detected spots to cells\n",
    "We have provided an optional exercise that demonstrates some more advanced napari features. If you do not have time during the practical session for this exercise, please feel free to perform it at a later time and ask the instructors any questions you may have via the class forum.\n",
    "\n",
    "In this exercise, we will assign the detected spots to the nearest nucleus since as a proxy for assigning detected spots to cells.\n",
    "\n",
    "\n",
    "### Step 1: get the locations of the nuclei\n",
    "As we saw at the end of the cellpose section above, the segmentations are stored in the Labels layer called `nuclei_cp_masks_000`. In the labels image, each each nucleus is given a unique integer label. We can access the label image from viewer as shown below.\n",
    "\n",
    "The Viewer object has a property calls `layers`, which is the `LayersList`. As the name suggests, the `Viewer.layers` contains each layer that is being displayed in the viewer and the layers can be accessed by their name. From the Labels layer object (called `nuclei_labels_layer` below), we can then access label image from the `data` property (stored as `nuclei_labels` below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the nuclei masks from the output of cellpose\n",
    "nuclei_labels_layer = viewer.layers['nuclei_cp_masks_000']\n",
    "nuclei_labels = nuclei_labels_layer.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the label image for the segmented nuclei, we can determine the centroid for each nucleus. First, use regionprops_table to find the centroids of all detected nuclei. Then, create an N x 2 numpy array containing the coordinate of the centroid of each nucleus. You can visualize the coordinates in a Points layer to verify they are correct.\n",
    "\n",
    "For more details on the regionprops_table see the scikit-image [docs](https://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.regionprops_table) and [examples](https://scikit-image.org/docs/stable/auto_examples/segmentation/plot_regionprops.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from skimage.measure import regionprops_table\n",
    "\n",
    "# create a regionprops table from the nuclei_labels\n",
    "# containing the label values and the centroids\n",
    "# of the nuclei\n",
    "rp_table = regionprops_table(\n",
    "    nuclei_labels,\n",
    "    properties=('label', 'centroid')\n",
    ")\n",
    "\n",
    "# create a numpy array containing the coordinates\n",
    "# of the centroids of the nuclei. Each row should\n",
    "# correspond to a nucleus and the columns should be\n",
    "# for the coordinate in the first and second axis\n",
    "nuclei_centroids = np.column_stack(\n",
    "    (rp_table['centroid-0'], rp_table['centroid-1'])\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cell below, add the nuclei centroids to the viewer as a Points layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Find the nearest nucleus for each spot\n",
    "\n",
    "In this step, we will use a KDTree to quickly determine which nucleus (`nuclei_centroids`) each spot (`spot_coords`) is closest to. For more information on the scipy cKDTree, please see [the docs](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html) and [query examples](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.query.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial import cKDTree\n",
    "\n",
    "# create the KDTree containing the nuclei locations\n",
    "kdt = cKDTree(nuclei_centroids)\n",
    "\n",
    "# query kdtree to get the nearest nucleus for each detected spot\n",
    "_, nearest_nucleus_ind = kdt.query(spot_coords, k=1)\n",
    "\n",
    "# convert the index to label value in the segmentation image\n",
    "nucleus_label_values = rp_table['label']\n",
    "nearest_nucleus_label = nucleus_label_values[nearest_nucleus_ind]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: visualize the results\n",
    "\n",
    "In this step, we will visualize how the detected spots were assigned to nuclei. To do so, we will only view the nuclei with spots and color the points to match the color of their respective nuclei.\n",
    "\n",
    "To do so, we will make use of the properties feature of the napari Points layer. `Points.properties` is a table where each row corresponds with a point (index matched to the coordinates) and each column is a feature of that point. The color of points can be mapped to their properties and we will use this to give all points assigned to the same nucleus the same color. To see an implementation of coloring points with properties, please see [this example](https://github.com/napari/napari/blob/master/examples/add_points_with_properties.py).\n",
    "\n",
    "The properties are stored as a dictionary where each key is the name of the column and each value is a numpy array containing the values for that column. For example, if we had 3 points with properties table as shown in the table below,\n",
    "\n",
    "\n",
    "| snr | confidence |\n",
    "| --- | --- |\n",
    "| 1 | 0.1 |\n",
    "| 3 | 0.5 |\n",
    "| 10 | 0.9 |\n",
    "\n",
    "we would create the properties as follows:\n",
    "\n",
    "```python\n",
    "properties = {\n",
    "    'snr': np.array([1, 3, 10]),\n",
    "    'confidence': np.array([0.1, 0.5, 0.9])\n",
    "}\n",
    "```\n",
    "\n",
    "In this case, we will create a properties table with a `nucleus_label` column, which is the label value for the nucleus to which that spot has been assigned. That is the key is `nucleus_label` and the value is `nearest_nucleus_label`, which was defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a properties table that maps each detected spot\n",
    "# to a nucleus label.\n",
    "spot_properties = {'nucleus_label': nearest_nucleus_label}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we define the face color parameters that we will use to create the new points layer. We will color the points using a color cycle. A color cycle is a list of colors that is applied in order and will repeat once all colors have been used.\n",
    "\n",
    "We define the face color options using a dictionary where the key is the setting name and the value is the setting value. See below for an explanation of the settings.\n",
    "\n",
    "- `colors`: these are the colors that will be applied to the faces of the points. Here, we are providing the property name (i.e., key in the properties dictionary) that we will map the colors to. This can also be a single color that is applied to all points or an array of colors (one for each point).\n",
    "- `color_mode`: select mode for how the face color is set. Here, we choose `cycle`, so that the face color is set by color cycle.\n",
    "- `categorical_colormap`: we provide the color cycle that will be used to set the face color."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the color cycle as a list of RGBA colors\n",
    "color_cycle = [\n",
    "    [0.12156863, 0.46666667, 0.70588235, 1.],\n",
    "    [1.        , 0.49803922, 0.05490196, 1.],\n",
    "    [0.17254902, 0.62745098, 0.17254902, 1.],\n",
    "    [0.83921569, 0.15294118, 0.15686275, 1.],\n",
    "    [0.58039216, 0.40392157, 0.74117647, 1.],\n",
    "    [0.54901961, 0.3372549 , 0.29411765, 1.],\n",
    "    [0.89019608, 0.46666667, 0.76078431, 1.],\n",
    "    [0.49803922, 0.49803922, 0.49803922, 1.],\n",
    "    [0.7372549 , 0.74117647, 0.13333333, 1.],\n",
    "    [0.09019608, 0.74509804, 0.81176471, 1.]\n",
    "]\n",
    "\n",
    "# define the face color options\n",
    "face_color = {\n",
    "    'colors': 'nucleus_label',\n",
    "    'color_mode': 'cycle',\n",
    "    'categorical_colormap': color_cycle\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can add the points to the viewer. As before, we use the `add_points()` method. However, this time, we have additional keyword arguments to pass the properties and face color settings we defined in the cells above.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pts = viewer.add_points(\n",
    "    spot_coords,\n",
    "    properties=spot_properties,\n",
    "    face_color=face_color,\n",
    "    symbol='ring',\n",
    "    size=4,\n",
    "    name='assigned_spots'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another thing to note is that `viewer.add_points()` returned the variable we assigned to `pts`. From the line below, we see that `pts` is the Points layer object that was created. This can be useful when further operations need to be performed on a newly created layer. In fact, all of the `viewer.add_*()` methods (e.g., `add_image()`) return the new layer object!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pts is viewer.layers['assigned_spots'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cell below, we inspect the colormap that is used to set the face colors of the points layer. We see that is it is a dictionary where the keys are the nucleus labels and the values are the corresponding colors. In the following cell, we will use this same mapping to make the labels image of the segmented nuclei use the same colormap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the mapping of nucleus label to color\n",
    "# from our points layer\n",
    "spots_cmap = pts._face.categorical_colormap.colormap\n",
    "spots_cmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the final step, we set the colors of the nuclei from our segmentation (`nuclei_cp_masks_000`) to match the spot colors. To do so, we use the colormapping we got from the points layer in the cell above (`spots_cmap`) and apply it to our nuclei labels layer via the `Labels.color` property. Setting the `color` property with a dictionary mapping label values to colors (as in `spots_cmap`) will apply that color map to the displayed labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply the same color mapping to our labels layer\n",
    "nuclei_labels_layer = viewer.layers['nuclei_cp_masks_000']\n",
    "nuclei_labels_layer.color = spots_cmap"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}