notebooks/DARTFISH.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproduce DARTFISH results with Starfish\n",
    "\n",
    "DARTFISH is a multiplexed image based transcriptomics assay from the [Zhang lab](http://genome-tech.ucsd.edu/ZhangLab/). As of this writing, this assay is not published yet. Nevertheless, here we demonstrate that Starfish can be used to process the data from raw images into spatially resolved gene expression profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import os\n",
    "\n",
    "from starfish import display\n",
    "from starfish import data, FieldOfView\n",
    "from starfish.types import Axes, Features, Levels\n",
    "\n",
    "from starfish import IntensityTable\n",
    "\n",
    "from starfish.image import Filter\n",
    "from starfish.spots import DetectPixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data into Starfish from the Cloud\n",
    "\n",
    "The example data here corresopond to DARTFISHv1 2017. The group is actively working on improving the protocol. The data represent human brain tissue from the human occipital cortex from 1 field of view (FOV) of larger experiment. The data from one field of view correspond to 18 images from 6 imaging rounds (r) 3 color channels (c) and 1 z-plane (z). Each image is 988x988 (y,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_test_data = os.getenv(\"USE_TEST_DATA\") is not None\n",
    "experiment = data.DARTFISH(use_test_data=use_test_data)\n",
    "\n",
    "imgs = experiment.fov().get_image(FieldOfView.PRIMARY_IMAGES)\n",
    "\n",
    "print(imgs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize codebook\n",
    "\n",
    "The DARTFISH codebook maps pixel intensities across the rounds and channels to the corresponding barcodes and genes that those pixels code for. For this example dataset, the codebook specifies 96 possible barcodes. The codebook used in this experiment has 3 color channels and one blank channel, each of which contribute to codes. The presence of the blank channel will be important later when the filtering is described."
   ]
  },
  {
   "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": [
    "# %gui qt5\n",
    "# display(stack)\n",
    "\n",
    "single_plane = imgs.sel({Axes.ROUND: 0, Axes.CH: 0, Axes.ZPLANE: 0})\n",
    "single_plane = single_plane.xarray.squeeze()\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(single_plane, cmap='gray', clim=list(np.percentile(single_plane, [1, 99.9])))\n",
    "plt.title('Round: 0, Chanel:0')\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter and scale raw data before decoding into spatially resolved gene expression\n",
    "\n",
    "First, we equalize the intensity of the images by scaling each image by its maximum intensity, which is equivalent to scaling by the 100th percentile value of the pixel values in each image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc_filt = Filter.Clip(p_max=100, level_method=Levels.SCALE_BY_CHUNK)\n",
    "norm_imgs = sc_filt.run(imgs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, for each imaging round, and each pixel location, we zero the intensity values across all three color channels if the magnitude of this 3 vector is below a threshold. As such, the code value associated with these pixels will be the blank. This is necessary to support euclidean decoding for codebooks that include blank values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_filt = Filter.ZeroByChannelMagnitude(thresh=.05, normalize=False)\n",
    "filtered_imgs = z_filt.run(norm_imgs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decode the processed data into spatially resolved gene expression profiles\n",
    "\n",
    "Here, starfish decodes each pixel value, across all rounds and channels, into the corresponding target (gene) it corresponds too. Contiguous pixels that map to the same target gene are called as one RNA molecule. Intuitively, pixel vectors are matched to the codebook by computing the euclidean distance between the pixel vector and all codewords. The minimal distance gene target is then selected, if it is within `distance_threshold` of any code.\n",
    "\n",
    "This decoding operation requires some parameter tuning, which is described below. First, we look at a distribution of pixel vector barcode magnitudes to determine the minimum magnitude threshold at which we will attempt to decode the pixel vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_magnitudes(stack, norm_order=2):\n",
    "\n",
    "    pixel_intensities = IntensityTable.from_image_stack(stack)\n",
    "    feature_traces = pixel_intensities.stack(traces=(Axes.CH.value, Axes.ROUND.value))\n",
    "    norm = np.linalg.norm(feature_traces.values, ord=norm_order, axis=1)\n",
    "\n",
    "    return norm\n",
    "\n",
    "mags = compute_magnitudes(filtered_imgs)\n",
    "\n",
    "plt.hist(mags, bins=20);\n",
    "sns.despine(offset=3)\n",
    "plt.xlabel('Barcode magnitude')\n",
    "plt.ylabel('Number of pixels')\n",
    "plt.yscale('log');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we decode the the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# how much magnitude should a barcode have for it to be considered by decoding? this was set by looking at\n",
    "# the plot above\n",
    "magnitude_threshold = 0.5\n",
    "# how big do we expect our spots to me, min/max size. this was set to be equivalent to the parameters\n",
    "# determined by the Zhang lab.\n",
    "area_threshold = (5, 30)\n",
    "# how close, in euclidean space, should the pixel barcode be to the nearest barcode it was called to?\n",
    "# here, I set this to be a large number, so I can inspect the distribution of decoded distances below\n",
    "distance_threshold = 3\n",
    "\n",
    "psd = DetectPixels.PixelSpotDecoder(\n",
    "    codebook=experiment.codebook,\n",
    "    metric='euclidean',\n",
    "    distance_threshold=distance_threshold,\n",
    "    magnitude_threshold=magnitude_threshold,\n",
    "    min_area=area_threshold[0],\n",
    "    max_area=area_threshold[1]\n",
    ")\n",
    "\n",
    "initial_spot_intensities, results = psd.run(filtered_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spots_df = initial_spot_intensities.to_features_dataframe()\n",
    "spots_df['area'] = np.pi*spots_df['radius']**2\n",
    "spots_df = spots_df.loc[spots_df[Features.PASSES_THRESHOLDS]]\n",
    "spots_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare to benchmark results\n",
    "\n",
    "The below plot aggregates gene copy number across cells in the field of view and compares the results to the same copy numbers from the authors' pipeline. This can likely be improved by tweaking parameters in the above algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load results from authors' pipeline\n",
    "# instead of just calling read_csv with the url, we are using python requests to load it to avoid a\n",
    "# SSL certificate error on some platforms.\n",
    "import io, requests\n",
    "cnts_benchmark = pd.read_csv(io.BytesIO(requests.get('https://d2nhj9g34unfro.cloudfront.net/20181005/DARTFISH/fov_001/counts.csv').content))\n",
    "cnts_benchmark.head()\n",
    "\n",
    "# select spots with distance less than a threshold, and count the number of each target gene\n",
    "min_dist = 0.6\n",
    "cnts_starfish = spots_df[spots_df.distance<=min_dist].groupby('target').count()['area']\n",
    "cnts_starfish = cnts_starfish.reset_index(level=0)\n",
    "cnts_starfish.rename(columns = {'target':'gene', 'area':'cnt_starfish'}, inplace=True)\n",
    "\n",
    "benchmark_comparison = pd.merge(cnts_benchmark, cnts_starfish, on='gene', how='left')\n",
    "benchmark_comparison.head(20)\n",
    "\n",
    "x = benchmark_comparison.dropna().cnt.values\n",
    "y = benchmark_comparison.dropna().cnt_starfish.values\n",
    "r = np.corrcoef(x, y)\n",
    "r = r[0,1]\n",
    "\n",
    "plt.scatter(x, y, 50,zorder=2)\n",
    "\n",
    "plt.xlabel('Gene copy number Benchmark')\n",
    "plt.ylabel('Gene copy number Starfish')\n",
    "plt.title('r = {}'.format(r))\n",
    "\n",
    "sns.despine(offset=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize results\n",
    "\n",
    "This image applies a pseudo-color to each gene channel to visualize the position and size of all called spots in the test image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exclude spots that don't meet our area thresholds\n",
    "area_lookup = lambda x: 0 if x == 0 else results.region_properties[x - 1].area\n",
    "vfunc = np.vectorize(area_lookup)\n",
    "mask = np.squeeze(vfunc(results.label_image))\n",
    "new_image = np.squeeze(results.decoded_image)*(mask > area_threshold[0])*(mask < area_threshold[1])\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(new_image, cmap = 'nipy_spectral');\n",
    "plt.axis('off');\n",
    "plt.title('Coded rolonies');\n",
    "\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib.patches import Rectangle\n",
    "\n",
    "rect = [Rectangle((100, 600), width=200, height=200)]\n",
    "pc = PatchCollection(rect, facecolor='none', alpha=1.0, edgecolor='w', linewidth=1.5)\n",
    "plt.gca().add_collection(pc)\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(new_image[600:800, 100:300], cmap = 'nipy_spectral');\n",
    "plt.axis('off');\n",
    "plt.title('Coded rolonies, zoomed in');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter and QC analysis\n",
    "\n",
    "Here, we further investigate reasonable choices for each of the parameters used by the PixelSpotDecoder. By tuning these parameters, one can achieve different results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,3))\n",
    "\n",
    "plt.subplot(131)\n",
    "plt.hist(mags, bins=100);\n",
    "plt.yscale('log')\n",
    "plt.xlabel('barcode magnitude')\n",
    "plt.ylabel('number of pixels')\n",
    "sns.despine(offset=2)\n",
    "plt.vlines(magnitude_threshold, ymin=plt.gca().get_ylim()[0], ymax=plt.gca().get_ylim()[1])\n",
    "plt.title('Set magnitude threshod')\n",
    "\n",
    "plt.subplot(132)\n",
    "spots_df['area'] = np.pi*spots_df.radius**2\n",
    "spots_df.area.hist(bins=30);\n",
    "plt.xlabel('area')\n",
    "plt.ylabel('number of spots')\n",
    "sns.despine(offset=2)\n",
    "plt.title('Set area threshold')\n",
    "\n",
    "plt.subplot(133)\n",
    "spots_df.distance.hist(bins=30)\n",
    "plt.xlabel('min distance to code');\n",
    "plt.vlines(min_dist, ymin=plt.gca().get_ylim()[0], ymax=plt.gca().get_ylim()[1])\n",
    "sns.despine(offset=2)\n",
    "plt.title('Set minimum distance threshold');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_threshold = min_dist\n",
    "\n",
    "psd = DetectPixels.PixelSpotDecoder(\n",
    "    codebook=experiment.codebook,\n",
    "    metric='euclidean',\n",
    "    distance_threshold=distance_threshold,\n",
    "    magnitude_threshold=magnitude_threshold,\n",
    "    min_area=area_threshold[0],\n",
    "    max_area=area_threshold[1]\n",
    ")\n",
    "\n",
    "spot_intensities, results = psd.run(filtered_imgs)\n",
    "spot_intensities = IntensityTable(spot_intensities.where(spot_intensities[Features.PASSES_THRESHOLDS], drop=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we:\n",
    "\n",
    "1. Pick a rolony that was succesfully decoded to a gene.\n",
    "2. Pull out the average pixel trace for that rolony\n",
    "3. Plot that pixel trace against the barcode of that gene\n",
    "\n",
    "In order to assess, visually, how close decoded barcodes match their targets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reshape the spot intensity table into a RxC barcode vector\n",
    "pixel_traces = spot_intensities.stack(traces=(Axes.ROUND.value, Axes.CH.value))\n",
    "\n",
    "# extract dataframe from spot intensity table for indexing purposes\n",
    "pixel_traces_df = pixel_traces.to_features_dataframe()\n",
    "pixel_traces_df['area'] = np.pi*pixel_traces_df.radius**2\n",
    "\n",
    "# pick index of a barcode that was read and decoded from the ImageStack\n",
    "ind = 4\n",
    "\n",
    "# get the the corresponding gene this barcode was decoded to\n",
    "gene = pixel_traces_df.loc[ind].target\n",
    "\n",
    "# query the codebook for the actual barcode corresponding to this gene\n",
    "real_barcode = experiment.codebook[experiment.codebook.target==gene].stack(traces=(Axes.ROUND.value, Axes.CH.value)).values[0]\n",
    "read_out_barcode = pixel_traces[ind,:]\n",
    "\n",
    "plt.plot(real_barcode, 'ok')\n",
    "plt.stem(read_out_barcode)\n",
    "sns.despine(offset=2)\n",
    "plt.xticks(range(18))\n",
    "plt.title(gene)\n",
    "plt.xlabel('Index into R (0:5) and C(0:2)');"
   ]
  }
 ],
 "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
}