notebooks/ISS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproduce In-situ Sequencing results with Starfish\n",
    "\n",
    "In Situ Sequencing (ISS) is an image based transcriptomics technique that can spatially resolve hundreds RNA species and their expression levels in-situ. The protocol and data analysis are described in this [publication](https://www.ncbi.nlm.nih.gov/pubmed/23852452). This notebook walks through how to use Starfish to process the raw images from an ISS experiment into a spatially resolved cell by gene expression matrix. We verify that Starfish can accurately reproduce the results from the authors' original [pipeline](https://cellprofiler.org/previous_examples/#sequencing-rna-molecules-in-situ-combining-cellprofiler-with-imagej-plugins)\n",
    "\n",
    "Please see [documentation](https://spacetx-starfish.readthedocs.io/en/stable/) for detailed descriptions of all the data structures and methods used here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pprint\n",
    "\n",
    "from starfish import data, FieldOfView\n",
    "from starfish.types import Axes, Features, FunctionSource\n",
    "from starfish.util.plot import imshow_plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams[\"figure.dpi\"] = 150"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Data into Starfish from the Cloud\n",
    "\n",
    "The primary data from one field of view correspond to 16 images from 4 hybridzation rounds (r) 4 color channels (c) one z plane (z). Each image is 1044 x 1390 (y, x). These data arise from human breast tissue. O(10) transcripts are barcoded for subsequent spatial resolution. Average pixel intensity values for one 'spot' in the image, across all rounds and channels, can be decoded into the nearest barcode, thus resolving each pixel into a particular gene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_test_data = os.getenv(\"USE_TEST_DATA\") is not None\n",
    "\n",
    "# An experiment contains a codebook, primary images, and auxiliary images\n",
    "experiment = data.ISS(use_test_data=use_test_data)\n",
    "pp = pprint.PrettyPrinter(indent=2)\n",
    "pp.pprint(experiment._src_doc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fov = experiment.fov()\n",
    "\n",
    "# note the structure of the 5D tensor containing the raw imaging data\n",
    "imgs = fov.get_image(FieldOfView.PRIMARY_IMAGES)\n",
    "print(imgs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize Codebook\n",
    "\n",
    " The ISS codebook maps each barcode to a gene. This protocol asserts that genes are encoded with\n",
    " a length 4 quatenary barcode that can be read out from the images. Each round encodes a position in the codeword.\n",
    " The maximum signal in each color channel (columns in the above image) corresponds to a letter in the codeword.\n",
    " The channels, in order, correspond to the letters: 'T', 'G', 'C', 'A'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment.codebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize raw data\n",
    "\n",
    "A nice way to page through all this data is to use the display command. We have commented this out for now, because it will not render in Github. Instead, we simply show an image from the first round and color channel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Display all the data in an interactive pop-up window. Uncomment to have this version work.\n",
    "# %gui qt5\n",
    "# display(imgs)\n",
    "\n",
    "# Display a single plane of data\n",
    "sel={Axes.ROUND: 0, Axes.CH: 0, Axes.ZPLANE: 0}\n",
    "single_plane = imgs.sel(sel)\n",
    "imshow_plane(single_plane, title=\"Round: 0, Channel: 0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "'dots' is a general stain for all possible transcripts. This image should correspond to the maximum projcection of all color channels within a single imaging round. This auxiliary image is useful for registering images from multiple imaging rounds to this reference image. We'll see an example of this further on in the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from starfish.image import Filter\n",
    "\n",
    "dots = fov.get_image(\"dots\")\n",
    "dots_single_plane = dots.reduce({Axes.ROUND, Axes.CH, Axes.ZPLANE}, func=\"max\")\n",
    "imshow_plane(dots_single_plane, title=\"Anchor channel, all RNA molecules\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is a DAPI image, which specifically marks nuclei. This is useful cell segmentation later on in the processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nuclei = fov.get_image(\"nuclei\")\n",
    "nuclei_single_plane = nuclei.reduce({Axes.ROUND, Axes.CH, Axes.ZPLANE}, func=\"max\")\n",
    "imshow_plane(nuclei_single_plane, title=\"Nuclei (DAPI) channel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter raw data before decoding into spatially resolved gene expression\n",
    "\n",
    "A White-Tophat filter can be used to enhance spots while minimizing background autoflourescence. The ```masking_radius``` parameter specifies the expected radius, in pixels, of each spot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter raw data\n",
    "masking_radius = 15\n",
    "filt = Filter.WhiteTophat(masking_radius, is_volume=False)\n",
    "\n",
    "filtered_imgs = filt.run(imgs, verbose=True, in_place=False)\n",
    "filt.run(dots, verbose=True, in_place=True)\n",
    "filt.run(nuclei, verbose=True, in_place=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_plane_filtered = filtered_imgs.sel(sel)\n",
    "\n",
    "f, (ax1, ax2) = plt.subplots(ncols=2)\n",
    "vmin, vmax = np.percentile(single_plane.xarray.values.data, [5, 99])\n",
    "imshow_plane(\n",
    "    single_plane, ax=ax1, vmin=vmin, vmax=vmax,\n",
    "    title=\"Original data\\nRound: 0, Channel: 0\"\n",
    ")\n",
    "vmin, vmax = np.percentile(single_plane_filtered.xarray.values.data, [5, 99])\n",
    "imshow_plane(\n",
    "    single_plane_filtered, ax=ax2, vmin=vmin, vmax=vmax,\n",
    "    title=\"Filtered data\\nRound: 0, Channel: 0\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Register data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Images may have shifted between imaging rounds. This needs to be corrected for before decoding, since this shift in the images will corrupt the barcodes, thus hindering decoding accuracy. A simple procedure can correct for this shift. For each imaging round, the max projection across color channels should look like the dots stain. Below, we simply shift all images in each round to match the dots stain by learning the shift that maximizes the cross-correlation between the images and the dots stain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from starfish.image import ApplyTransform, LearnTransform\n",
    "\n",
    "learn_translation = LearnTransform.Translation(reference_stack=dots, axes=Axes.ROUND, upsampling=1000)\n",
    "transforms_list = learn_translation.run(imgs.reduce({Axes.CH, Axes.ZPLANE}, func=\"max\"))\n",
    "warp = ApplyTransform.Warp()\n",
    "registered_imgs = warp.run(filtered_imgs, transforms_list=transforms_list, in_place=False, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decode the processed data into spatially resolved gene expression profiles\n",
    "\n",
    "To decode, first we find spots, and record, for reach spot, the average pixel intensities across rounds and channels. This spot detection can be achieved by the ```BlobDetector``` algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "from starfish.spots import FindSpots, DecodeSpots\n",
    "\n",
    "bd = FindSpots.BlobDetector(\n",
    "    min_sigma=1,\n",
    "    max_sigma=10,\n",
    "    num_sigma=30,\n",
    "    threshold=0.01,\n",
    "    measurement_type='mean',\n",
    ")\n",
    "\n",
    "dots_max = dots.reduce((Axes.ROUND, Axes.ZPLANE), func=FunctionSource.np(\"max\"))\n",
    "spots = bd.run(image_stack=registered_imgs, reference_image=dots_max)\n",
    "\n",
    "decoder = DecodeSpots.PerRoundMaxChannel(codebook=experiment.codebook)\n",
    "decoded = decoder.run(spots=spots)\n",
    "\n",
    "# Besides house keeping genes, VIM and HER2 should be most highly expessed, which is consistent here.\n",
    "genes, counts = np.unique(decoded.loc[decoded[Features.PASSES_THRESHOLDS]][Features.TARGET], return_counts=True)\n",
    "table = pd.Series(counts, index=genes).sort_values(ascending=False)\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Segment Cells and create Cell by Gene Expression Matrix\n",
    "\n",
    "After calling spots and decoding their gene information, cells must be segmented to assign genes to cells. This paper used a seeded watershed approach to segment the cells, which we also use here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from starfish.morphology import Binarize, Filter, Merge, Segment\n",
    "from starfish.types import Levels\n",
    "\n",
    "dapi_thresh = .18  # binary mask for cell (nuclear) locations\n",
    "stain_thresh = .22  # binary mask for overall cells // binarization of stain\n",
    "min_dist = 57\n",
    "min_allowed_size = 10\n",
    "max_allowed_size = 10000\n",
    "\n",
    "mp = registered_imgs.reduce({Axes.CH, Axes.ZPLANE}, func=\"max\")\n",
    "stain = mp.reduce(\n",
    "    {Axes.ROUND},\n",
    "    func=\"mean\",\n",
    "    level_method=Levels.SCALE_BY_IMAGE)\n",
    "\n",
    "nuclei_mp_scaled = nuclei.reduce(\n",
    "    {Axes.ROUND, Axes.CH, Axes.ZPLANE},\n",
    "    func=\"max\",\n",
    "    level_method=Levels.SCALE_BY_IMAGE)\n",
    "\n",
    "binarized_nuclei = Binarize.ThresholdBinarize(dapi_thresh).run(nuclei_mp_scaled)\n",
    "labeled_masks = Filter.MinDistanceLabel(min_dist, 1).run(binarized_nuclei)\n",
    "watershed_markers = Filter.AreaFilter(min_area=min_allowed_size, max_area=max_allowed_size).run(labeled_masks)\n",
    "thresholded_stain = Binarize.ThresholdBinarize(stain_thresh).run(stain)\n",
    "markers_and_stain = Merge.SimpleMerge().run([thresholded_stain, watershed_markers])\n",
    "watershed_mask = Filter.Reduce(\n",
    "    \"logical_or\",\n",
    "    lambda shape: np.zeros(shape=shape, dtype=bool)\n",
    ").run(markers_and_stain)\n",
    "\n",
    "segmenter = Segment.WatershedSegment(connectivity=np.ones((1, 3, 3), dtype=bool))\n",
    "masks = segmenter.run(\n",
    "    stain,\n",
    "    watershed_markers,\n",
    "    watershed_mask,\n",
    ")\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from showit import image\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "\n",
    "plt.subplot(321)\n",
    "nuclei_numpy = nuclei_mp_scaled._squeezed_numpy(Axes.ROUND, Axes.CH, Axes.ZPLANE)\n",
    "image(nuclei_numpy, ax=plt.gca(), size=20, bar=True)\n",
    "plt.title('Nuclei')\n",
    "\n",
    "plt.subplot(322)\n",
    "image(\n",
    "    stain._squeezed_numpy(Axes.ROUND, Axes.CH, Axes.ZPLANE),\n",
    "    ax=plt.gca(), size=20, bar=True)\n",
    "plt.title('Stain')\n",
    "\n",
    "plt.subplot(323)\n",
    "image(\n",
    "    binarized_nuclei.uncropped_mask(0).squeeze(Axes.ZPLANE.value).values,\n",
    "    bar=False,\n",
    "    ax=plt.gca(),\n",
    ")\n",
    "plt.title('Nuclei Thresholded')\n",
    "\n",
    "plt.subplot(324)\n",
    "image(\n",
    "    watershed_mask.to_label_image().xarray.squeeze(Axes.ZPLANE.value).values,\n",
    "    bar=False,\n",
    "    ax=plt.gca(),\n",
    ")\n",
    "plt.title('Watershed Mask')\n",
    "\n",
    "plt.subplot(325)\n",
    "image(\n",
    "    watershed_markers.to_label_image().xarray.squeeze(Axes.ZPLANE.value).values,\n",
    "    size=20,\n",
    "    cmap=plt.cm.nipy_spectral,\n",
    "    ax=plt.gca(),\n",
    ")\n",
    "plt.title('Found: {} cells'.format(len(watershed_markers)))\n",
    "\n",
    "plt.subplot(326)\n",
    "image(\n",
    "    masks.to_label_image().xarray.squeeze(Axes.ZPLANE.value).values,\n",
    "    size=20,\n",
    "    cmap=plt.cm.nipy_spectral,\n",
    "    ax=plt.gca(),\n",
    ")\n",
    "plt.title('Segmented Cells')\n",
    "plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that cells have been segmented, we can assign spots to cells in order to create a cell x gene count matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from starfish.spots import AssignTargets\n",
    "from starfish import ExpressionMatrix\n",
    "\n",
    "al = AssignTargets.Label()\n",
    "labeled = al.run(masks, decoded)\n",
    "cg = labeled.to_expression_matrix()\n",
    "cg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Compare to results from paper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This FOV was selected to make sure that we can visualize the tumor/stroma boundary, below this is described by pseudo-coloring HER2 (tumor) and vimentin (VIM, stroma). This distribution matches the one described in the original paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.color import rgb2gray\n",
    "\n",
    "GENE1 = 'HER2'\n",
    "GENE2 = 'VIM'\n",
    "\n",
    "rgb = np.zeros(registered_imgs.tile_shape + (3,))\n",
    "nuclei_numpy = nuclei.reduce({Axes.ROUND, Axes.CH, Axes.ZPLANE}, func=\"max\")._squeezed_numpy(Axes.ROUND, Axes.CH, Axes.ZPLANE)\n",
    "rgb[:,:,0] = nuclei_numpy\n",
    "dots_numpy = dots.reduce({Axes.ROUND, Axes.CH, Axes.ZPLANE}, func=\"max\")._squeezed_numpy(Axes.ROUND, Axes.CH, Axes.ZPLANE)\n",
    "rgb[:,:,1] = dots_numpy\n",
    "do = rgb2gray(rgb)\n",
    "do = do/(do.max())\n",
    "\n",
    "plt.imshow(do,cmap='gray')\n",
    "plt.axis('off');\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter('ignore', FutureWarning)\n",
    "    is_gene1 = decoded.where(decoded[Features.AXIS][Features.TARGET] == GENE1, drop=True)\n",
    "    is_gene2 = decoded.where(decoded[Features.AXIS][Features.TARGET] == GENE2, drop=True)\n",
    "\n",
    "plt.plot(is_gene1.x, is_gene1.y, 'or', markersize=3)\n",
    "plt.plot(is_gene2.x, is_gene2.y, 'ob', markersize=3)\n",
    "plt.title(f'Red: {GENE1}, Blue: {GENE2}');"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}