{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n\n# Spot-Based Decoding\n\nSpot-based decoding is the approach of finding spots in images from each round first and then\ndecoding them. The alternative, pixel-based decoding, decodes pixels first and then connects them\ninto spots after.\n\n\n\nStarfish provides multiple options for each component of spot-based decoding:\n\n## Spot Finding\n\n.. list-table:: :py:class:`.FindSpotsAlgorithm`\n :widths: auto\n :header-rows: 1\n\n * - Method\n - Description\n - Works in 3D\n - Finds Threshold\n - Finds Sigma\n - Anisotropic Sigma\n * - :py:class:`.BlobDetector`\n - Wrapper of classic kernel convolution blob detection algorithms in :py:mod:`skimage.feature`\n such as LoG, which uses the Laplacian of Gaussian filter\n - |yes|\n - |no|\n - |yes|\n - |yes|\n * - :py:class:`.LocalMaxPeakFinder`\n - Wrapper of :py:mod:`skimage.feature.peak_local_max`, which finds local maxima pixel\n intensities in an image\n - |yes|\n - |yes|\n - |no|\n - |no|\n * - :py:class:`.TrackpyLocalMaxPeakFinder`\n - Wrapper for :py:mod:`trackpy.locate`, which implements a version of the Crocker-Grier\n algorithm originally developed for particle tracking\n - |yes|\n - |no|\n - |no|\n - |yes|\n\n:py:class:`.BlobDetector` and :py:class:`.LocalMaxPeakFinder` should usually be chosen over\n:py:class:`.TrackpyLocalMaxPeakFinder`, and :py:class:`.BlobDetector` should be favored over\n:py:class:`.LocalMaxPeakFinder` if you are unsure of the size of the spot and the spots are\nuniformly gaussian in shape. :py:class:`.LocalMaxPeakFinder`, by contrast, can help find the correct\nminimum peak intensity threshold. A demonstration of their differences can be found at the end of\nthis tutorial.\n\nDetected spots are returned in :py:class:`.SpotFindingResults`, which can be\n`visually assessed ` before decoding.\n\n* `howto_blobdetector`\n* `howto_localmaxpeakfinder`\n* `howto_trackpylocalmaxpeakfinder`\n\n## Trace Building\n\n.. list-table:: ``TraceBuildingStrategy``\n :widths: auto\n :header-rows: 1\n\n * - Method\n - Description\n - Reference Image\n * - ``SEQUENTIAL``\n - Build traces for every detected spot by setting intensity values to zero for all rounds\n and channels the spot was not found in (i.e. every trace will have only 1 non-zero value)\n - Incompatible\n * - ``EXACT_MATCH``\n - Build traces by combining intensity values of spots from every rounds and channel in the\n exact same location as spots in ``reference_image``\n - Required\n * - ``NEAREST_NEIGHBOR``\n - Build traces by combining intensity values of spots from rounds and channels nearest to the\n spots in the ``anchor_round``\n - Not recommended; will have same result as EXACT_MATCH\n\nThe first step to decoding :py:class:`.SpotFindingResults` is identifying spots from\ndifferent imaging rounds as the same spot. In starfish this is referred to as building traces and it\ntransforms :py:class:`.SpotFindingResults` to an :py:class:`.IntensityTable`. Trace building is\nhidden in the :py:class:`.DecodeSpotsAlgorithm` but it requires the user to select a\n``TraceBuildingStrategy``. `howto_tracebuildingstrategies`. goes further in depth and shows\nhow to build traces independent of the decoding step.\n\n\n## Spot Decoding\n\n.. list-table:: :py:class:`.DecodeSpotsAlgorithm`\n :widths: auto\n :header-rows: 1\n\n * - Method\n - Description\n - Works in 3D\n - Codebook Design\n - TraceBuildingStrategy\n - Returns Quality Score\n * - :py:class:`.SimpleLookupDecoder`\n - Description\n - |yes|\n - Linearly multiplexed\n - Sequential\n - |no|\n * - :py:class:`.PerRoundMaxChannel`\n - Description\n - |yes|\n - One hot exponentially multiplexed\n - Sequential, Exact_Match or Nearest_Neighbor\n - |yes|\n * - :py:class:`.MetricDistance`\n - Description\n - |yes|\n - Exponentially multiplexed\n - Exact_Match or Nearest_Neighbor\n - |yes|\n\n.. |yes| unicode:: U+2705 .. White Heavy Check Mark\n.. |no| unicode:: U+274C .. Cross Mark\n\nStarfish decoding is done by running a :py:class:`.DecodeSpotsAlgorithm` on\n:py:class:`.SpotFindingResults` to return a :py:class:`.DecodedIntensityTable`.\n:py:class:`.PerRoundMaxChannel` should generally be used rather than the\nthe other two decoding algorithms if possible. :py:class:`.MetricDistance` is necessary for\n:term:`codebooks` that contain :term:`codewords` without exactly one hot\nchannel in every round. This is used for error-robustness (e.g. MERFISH) and/or reducing optical\ncrowding in each round.\n\n* `howto_simplelookupdecoder`\n* `howto_perroundmaxchannel`\n* `howto_metricdistance`\n\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparison of :py:class:`.FindSpotsAlgorithm`\\s\nThis tutorial demonstrates usage and characteristics of the three available\n:py:class:`.FindSpotsAlgorithm`\\s on 3-Dimensional STARmap images. Parameters were roughly\ntuned, but these results are not reflective of the best possible performance of each\n:py:class:`.FindSpotsAlgorithm`.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# load STARmap data\nfrom starfish import data, display\nfrom starfish.image import ApplyTransform, LearnTransform\nfrom starfish.spots import FindSpots\nfrom starfish.types import Axes\nfrom starfish.util.plot import imshow_plane\nfrom starfish.core.spots.DecodeSpots.trace_builders import build_spot_traces_exact_match\n\nexperiment = data.STARmap(use_test_data=True)\nimgs = experiment['fov_000'].get_image('primary')\n\n# register rounds\nprojection = imgs.reduce({Axes.CH, Axes.ZPLANE}, func=\"max\")\nreference_image = projection.sel({Axes.ROUND: 0})\nltt = LearnTransform.Translation(reference_stack=reference_image, axes=Axes.ROUND, upsampling=1000)\ntransforms = ltt.run(projection)\nwarp = ApplyTransform.Warp()\nimgs = warp.run(stack=imgs, transforms_list=transforms)\n\n# make reference round for decoding\ndots = imgs.reduce({Axes.CH, Axes.ROUND}, func=\"max\")\n\n# view a cropped region of image for spot finding\nimshow_plane(dots, sel={Axes.ZPLANE: 15, Axes.X: (400, 600), Axes.Y: (400, 600)})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The ``dots`` reference image shows spots are approximately 10 pixels in diameter and can be\ntightly packed together. This helped inform the parameter settings of the spot finders below.\nHowever, the accuracy can always be improved by further tuning parameters. For example,\nif low intensity background noise is being detected as spots, increasing values for\n``threshold``, ``stringency``, ``min_mass``, and ``percentile`` can remove them. If large spots\nare not missed, increasing values such as ``max_sigma``, ``max_obj_area``, ``spot_diameter``,\nand ``max_size`` could include them. Moreover, signal enhancement and background reduction\nprior to this step can also improve accuracy of spot finding.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"bd = FindSpots.BlobDetector(\n min_sigma=2,\n max_sigma=6,\n num_sigma=20,\n threshold=0.1,\n is_volume=True,\n measurement_type='mean',\n)\n\nlmp = FindSpots.LocalMaxPeakFinder(\n min_distance=2,\n stringency=8,\n min_obj_area=6,\n max_obj_area=600,\n is_volume=True\n)\n\ntlmpf = FindSpots.TrackpyLocalMaxPeakFinder(\n spot_diameter=11,\n min_mass=0.2,\n max_size=8,\n separation=3,\n preprocess=False,\n percentile=80,\n verbose=True,\n)\n\n# crop imagestacks\ncrop_selection = {Axes.X: (300, 700), Axes.Y: (300, 700)}\ncropped_imgs= imgs.sel(crop_selection)\ncropped_dots = dots.sel(crop_selection)\n\n# find spots on cropped images\nbd_spots = bd.run(image_stack=cropped_imgs, reference_image=cropped_dots)\nlmp_spots = lmp.run(image_stack=cropped_imgs, reference_image=cropped_dots)\ntlmpf_spots = tlmpf.run(image_stack=cropped_imgs, reference_image=cropped_dots)\n\n# build spot traces into intensity table\nbd_table = build_spot_traces_exact_match(bd_spots)\nlmp_table = build_spot_traces_exact_match(lmp_spots)\ntlmpf_table = build_spot_traces_exact_match(tlmpf_spots)\n\n# plot spots found\nimport matplotlib\nimport matplotlib.pyplot as plt\n\n# get x, y coords and spot size from intensity table\ndef get_cropped_coords(table, x_min, x_max, y_min, y_max):\n df = table.to_features_dataframe()\n df = df.loc[df['x'].between(x_min, x_max) & df['y'].between(y_min, y_max)]\n return df['x'].values-x_min, df['y'].values-y_min, df['radius'].values.astype(int)\nbd_x, bd_y, bd_s = get_cropped_coords(bd_table, 200, 250, 200, 250)\nlmp_x, lmp_y, lmp_s = get_cropped_coords(lmp_table, 200, 250, 200, 250)\ntlmpf_x, tlmpf_y, tlmpf_s = get_cropped_coords(tlmpf_table, 200, 250, 200, 250)\n\nmatplotlib.rcParams[\"figure.dpi\"] = 150\nf, (ax1, ax2, ax3) = plt.subplots(ncols=3)\n\n# Plot cropped region of max z-projected dots image\nimshow_plane(cropped_dots.reduce({Axes.ZPLANE}, func=\"max\"), sel={Axes.X: (200, 250), Axes.Y: (200, 250)}, ax=ax1, title='BlobDetector')\nimshow_plane(cropped_dots.reduce({Axes.ZPLANE}, func=\"max\"), sel={Axes.X: (200, 250), Axes.Y: (200, 250)}, ax=ax2, title='LocalMaxPeak')\nimshow_plane(cropped_dots.reduce({Axes.ZPLANE}, func=\"max\"), sel={Axes.X: (200, 250), Axes.Y: (200, 250)}, ax=ax3, title='Trackpy')\n# Overlay spots found by each FindSpotsAlgorithm\nax1.scatter(bd_x, bd_y, marker='o', facecolors='none', edgecolors='r', s=bd_s*10)\nax2.scatter(lmp_x, lmp_y, marker='o', facecolors='none', edgecolors='r', s=lmp_s*10)\nax3.scatter(tlmpf_x, tlmpf_y, marker='o', facecolors='none', edgecolors='r', s=tlmpf_s*10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These images show :py:class:`.BlobDetector` is good at finding round Gaussian spots with\nseparation from one another. However, multiple spots packed closely together may be detected as\none large spot. :py:class:`.LocalMaxPeakFinder` and :py:class:`.TrackpyLocalMaxPeakFinder` are\nmore sensitive to peaks of intensity, which allows them to detect multiple partially\noverlapping spots but also makes them more vulnerable to noise. Smoothing the image with a low\nband pass filter can mitigate that issue. The last thing to note is spots found by\n:py:class:`.LocalMaxPeakFinder` all have a radius of one, because unlike the other two\nalgorithms it isn't technically finding spots but just peaks of intensity.\n\nThe multi-dimensional spot finding results can also be viewed interactively with napari using\nthe code below:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# uncomment to view with napari\n# %gui qt\n# BlobDetector\n# viewer = display(stack=cropped_dots, spots=bd_table)\n\n# LocalMaxPeakFinder\n# display(stack=cropped_dots, spots=lmp_table)\n\n# TrackpyLocalMaxPeakFinder\n# display(stack=cropped_dots, spots=tlmpf_table)"
]
}
],
"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.6.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}