notebooks/smFISH.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Single Field of View for sequential single-molecule FISH processed in 3d\n",
    "========================================================================\n",
    "\n",
    "This notebook walks through a work flow that analyzes one field of view of a mouse gene panel from\n",
    "the Allen Institute for Brain Science, using the starfish package.\n",
    "\n",
    "This example processes an experiment with a single round from a single field of view of sequential\n",
    "smFISH data taken from mouse primary visual cortex. These data are unpublished, and were kindly\n",
    "contributed by the Allen Institute for Brain Science as a part of the SpaceTx consortium\n",
    "project.\n",
    "\n",
    "The data consist of 45 images from 1 round, 1 channels, and 33 z-planes. Each image is\n",
    "(2048, 2048) (y, x). There are no test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Optional, Tuple\n",
    "from IPython import get_ipython\n",
    "\n",
    "import starfish\n",
    "import starfish.data\n",
    "from starfish import FieldOfView, DecodedIntensityTable\n",
    "from starfish.types import TraceBuildingStrategies\n",
    "\n",
    "# equivalent to %gui qt\n",
    "ipython = get_ipython()\n",
    "ipython.magic(\"gui qt5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define image filters\n",
    "--------------------\n",
    "The 3d smFISH workflow run by the Allen runs a bandpass filter to remove high and low frequency\n",
    "signal and blurs over z with a 1-pixel gaussian to smooth the signal over the z-axis.\n",
    "\n",
    "low-intensity signal is (stringently) clipped from the images before and after these filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bandpass filter to remove cellular background and camera noise\n",
    "bandpass = starfish.image.Filter.Bandpass(lshort=.5, llong=7, threshold=0.0)\n",
    "\n",
    "# gaussian blur to smooth z-axis\n",
    "glp = starfish.image.Filter.GaussianLowPass(\n",
    "    sigma=(1, 0, 0),\n",
    "    is_volume=True\n",
    ")\n",
    "\n",
    "# pre-filter clip to remove low-intensity background signal\n",
    "clip1 = starfish.image.Filter.Clip(p_min=50, p_max=100)\n",
    "\n",
    "# post-filter clip to eliminate all but the highest-intensity peaks\n",
    "clip2 = starfish.image.Filter.Clip(p_min=99, p_max=100, is_volume=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a spot detection method\n",
    "------------------------------\n",
    "Spots are detected using a spot finder based on trackpy's locate method, which identifies\n",
    "local intensity maxima, and spots are matched to the gene they represent by looking them up in a\n",
    "codebook that records which (round, channel) matches which gene target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tlmpf = starfish.spots.FindSpots.TrackpyLocalMaxPeakFinder(\n",
    "    spot_diameter=5,  # must be odd integer\n",
    "    min_mass=0.02,\n",
    "    max_size=2,  # this is max radius\n",
    "    separation=7,\n",
    "    noise_size=0.65,  # this is not used because preprocess is False\n",
    "    preprocess=False,\n",
    "    percentile=10,  # this is irrelevant when min_mass, spot_diameter, and max_size are set properly\n",
    "    verbose=True,\n",
    "    is_volume=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the pipeline\n",
    "----------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# override print to print to stderr for cromwell\n",
    "from functools import partial\n",
    "import sys\n",
    "print = partial(print, file=sys.stderr)\n",
    "\n",
    "\n",
    "def processing_pipeline(\n",
    "    experiment: starfish.Experiment,\n",
    "    fov_name: str,\n",
    "    n_processes: Optional[int]=None\n",
    ") -> Tuple[starfish.ImageStack, starfish.IntensityTable]:\n",
    "    \"\"\"Process a single field of view of an experiment\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    experiment : starfish.Experiment\n",
    "        starfish experiment containing fields of view to analyze\n",
    "    fov_name : str\n",
    "        name of the field of view to process\n",
    "    n_processes : int\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    starfish.IntensityTable :\n",
    "        decoded IntensityTable containing spots matched to the genes they are hybridized against\n",
    "    \"\"\"\n",
    "\n",
    "    all_intensities = list()\n",
    "    codebook = experiment.codebook\n",
    "\n",
    "    print(\"Loading images...\")\n",
    "    images = enumerate(experiment[fov_name].get_images(FieldOfView.PRIMARY_IMAGES))\n",
    "\n",
    "    decoder = starfish.spots.DecodeSpots.PerRoundMaxChannel(\n",
    "        codebook=codebook,\n",
    "        trace_building_strategy=TraceBuildingStrategies.SEQUENTIAL\n",
    "    )\n",
    "\n",
    "    for image_number, primary_image in images:\n",
    "        print(f\"Filtering image {image_number}...\")\n",
    "        filter_kwargs = dict(\n",
    "            in_place=True,\n",
    "            verbose=True,\n",
    "            n_processes=n_processes\n",
    "        )\n",
    "        print(\"Applying Clip...\")\n",
    "        clip1.run(primary_image, **filter_kwargs)\n",
    "        print(\"Applying Bandpass...\")\n",
    "        bandpass.run(primary_image, **filter_kwargs)\n",
    "        print(\"Applying Gaussian Low Pass...\")\n",
    "        glp.run(primary_image, **filter_kwargs)\n",
    "        print(\"Applying Clip...\")\n",
    "        clip2.run(primary_image, **filter_kwargs)\n",
    "\n",
    "        print(\"Calling spots...\")\n",
    "        spots = tlmpf.run(primary_image)\n",
    "        print(\"Decoding spots...\")\n",
    "        decoded_intensities = decoder.run(spots=spots)\n",
    "        all_intensities.append(decoded_intensities)\n",
    "\n",
    "    decoded = DecodedIntensityTable.concatenate_intensity_tables(all_intensities)\n",
    "    decoded = decoded[decoded[\"total_intensity\"] > .025]\n",
    "\n",
    "    print(\"Processing complete.\")\n",
    "\n",
    "    return primary_image, decoded"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data, run pipeline, display results\n",
    "----------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment = starfish.data.allen_smFISH(use_test_data=True)\n",
    "\n",
    "image, intensities = processing_pipeline(experiment, fov_name='fov_001')\n",
    "\n",
    "# uncomment the below line to visualize the output with the spot calls.\n",
    "# viewer = starfish.display(image, intensities)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "starfish",
   "language": "python",
   "name": "starfish"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}