3 - Quantification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: Quantification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will look at ways to quantitatively characterize objects as they appear in microscopy images. Many different types of measurements can be extracted from images. The ones that are most relevant is ultimately dependent on the underlying biological question. Here, we will look at commonly-used (non-machine-learning) strategies to quantify information that pertains to image intensity and object shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import imageio.v2 as imageio\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams['figure.dpi'] = 200\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skimage as sk\n",
    "import skimage.measure as skm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this part, we will use pandas to handle numerical data (https://pandas.pydata.org/). Feel free to consult the extensive documentation available on their website if you want to know more about this library!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Data loading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will again work with images from the BBBC datasets introduced in notebook 1 - Handling Image Data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Load data\n",
    "\n",
    "Using what you learned in notebook 1 - Handling Image Data, load the following image files:\n",
    "- data/Part_3/BBBC010/raw.tif\n",
    "- data/Part_3/BBBC020/membranes.tif\n",
    "\n",
    "*Do not forget to convert the images into floating-point arrays with img_as_float!*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BBBC010\n",
    "##### Add you code here ##### \n",
    "\n",
    "# BBBC020\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Load segmentations\n",
    "\n",
    "To help us with the quantification, we also have instance segmentation masks provided for each image in a folder called \"instances\". Run the lines below to load these masks into an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BBBC010\n",
    "bbbc010_instances = {}\n",
    "for file in os.listdir('data/Part_3/BBBC010/instances'):\n",
    "    if file.endswith(\".png\"):\n",
    "        im = imageio.imread(os.path.join('data/Part_3/BBBC010/instances',file))\n",
    "        idx = file[:6]\n",
    "        bbbc010_instances[idx] = im\n",
    "                            \n",
    "# BBBC020\n",
    "bbbc020_instances = {}\n",
    "for file in os.listdir('data/Part_3/BBBC020/instances'):\n",
    "    if file.endswith(\".TIF\"):\n",
    "        im = imageio.imread(os.path.join('data/Part_3/BBBC020/instances',file))\n",
    "        idx = file[-9:-4]\n",
    "        bbbc020_instances[idx] = im"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Intensity-based measurements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most direct measurements one can make are those based on pixel values. When dealing with raw (unprocessed) images, pixel values are related to the amount of light (photons) that reached the camera during the acquisition process - although the exact nature of this relationship may be quite complicated. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Intensity stats\n",
    "\n",
    "The BBBC020 image corresponds to the signal of a fluorescently-tagged cell-surface protein (CD11b/APC). Run the code below to extract some statistics about the intensity distributions for each objects in the image using the provided instance segmentation masks. Are the results surprising? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract intensity features\n",
    "bbbc020_intensity_feats = { 'instance_id': [],\n",
    "                            'median': [],\n",
    "                            'standard_deviation': [] }\n",
    "\n",
    "for idx in bbbc020_instances.keys():\n",
    "    mask = bbbc020_instances[idx]\n",
    "    \n",
    "    intensities = bbbc020[mask>0] \n",
    "    bbbc020_intensity_feats['instance_id'].append(idx)\n",
    "    bbbc020_intensity_feats['median'].append(np.median(intensities))\n",
    "    bbbc020_intensity_feats['standard_deviation'].append(np.std(intensities))\n",
    "    \n",
    "bbbc020_intensity_feats = pd.DataFrame(bbbc020_intensity_feats)\n",
    "bbbc020_intensity_feats.set_index('instance_id', inplace = True)\n",
    "\n",
    "display(bbbc020_intensity_feats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the distribution of each feature as an histogram\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].hist(bbbc020_intensity_feats['median'].values) \n",
    "axes[0].set_title('Median')\n",
    "\n",
    "axes[1].hist(bbbc020_intensity_feats['standard_deviation'].values)  \n",
    "axes[1].set_title('Standard Deviation')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Analysis\n",
    "\n",
    "The BBBC010 image is a brightfield view of *C. elegans* worms treated with an antibiotic. The worms are also stained with a fluorescent protein that specifically targets dead animals. Run the code below to load the fluorescence microscopy readout, and adapt the code from 2.1 to study the population of worms in the BBBC010 image. Can you spot the dead worms?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load fluorescence readout\n",
    "bbbc010_fluo = imageio.imread('data/Part_3/BBBC010/fluorescence.tif')\n",
    "bbbc010_fluo = sk.img_as_float(bbbc010_fluo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract intensity features from fluorescence readout\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the distribution of each feature as an histogram\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Shape-based measurements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the intensity of pixels within them, the shape of objects in images can hold relevant information. Shape is generally challenging to quantify with a single number, but can be reasonably characterized with collections of descriptors as we shall see."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Region properties\n",
    "\n",
    "The most classical and direct way of characterizing shape is to rely on collection of measurements referred to as *region properties*. Run the code below to extract a subset of scikit-image's own collection of region properties on the BBBC020 image. Their definition is provided at https://scikit-image.org/docs/stable/api/skimage.measure.html#skimage.measure.regionprops . Do you understand which aspect of the shape they capture?\n",
    "\n",
    "Below, you will also find a piece of code that colors instances in the image according to the value of a shape property of your choice. Play around with it to get an intuition of what these different metrics reveal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract shape features\n",
    "bbbc020_shape_feats = pd.DataFrame()\n",
    "\n",
    "for idx in bbbc020_instances.keys():\n",
    "    mask = bbbc020_instances[idx]\n",
    "    \n",
    "    all_props = skm.regionprops_table(mask, properties = ('area', 'perimeter', 'axis_major_length', \n",
    "                                      'axis_minor_length', 'eccentricity', 'solidity', 'extent'))\n",
    "    all_props['instance_id'] = idx\n",
    "    \n",
    "    if bbbc020_shape_feats.empty:\n",
    "        bbbc020_shape_feats = pd.DataFrame(all_props)\n",
    "    else:\n",
    "        bbbc020_shape_feats = pd.concat([bbbc020_shape_feats,pd.DataFrame(all_props)], axis=0)\n",
    "\n",
    "bbbc020_shape_feats.set_index('instance_id', inplace = True)\n",
    "\n",
    "display(bbbc020_shape_feats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Label objects according to a property of your choice\n",
    "feature = 'area'\n",
    "\n",
    "feature_mask = np.zeros(bbbc020.shape)\n",
    "for idx in bbbc020_instances.keys():\n",
    "    mask = bbbc020_instances[idx]\n",
    "    \n",
    "    feature_val = bbbc020_shape_feats.loc[idx, feature]\n",
    "    feature_mask[np.where(mask)] = feature_val\n",
    "\n",
    "plt.imshow(bbbc020, cmap='gray')\n",
    "im = plt.imshow(feature_mask, cmap='inferno', alpha=.4)\n",
    "plt.colorbar(im, fraction=0.05, pad=0.01)\n",
    "plt.axis('off')\n",
    "plt.title('Color-coded by '+feature, fontsize=9)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Fourier descriptors\n",
    "\n",
    "As an alternative to individual independent measurements as exactred from region propertiers, one can summarize a shape by the coefficients of its decomposition into an orthonormal basis. This is the idea behind Fourier descriptors (https://en.wikipedia.org/wiki/Fourier_analysis), implemented in the pyefd library (https://github.com/hbldh/pyefd)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyefd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the technical details of Fourier analysis is beyond the scope of this course, we can get a visual intuition of how that works. Run the code below to 1) decompose the shape of one of the objects in the BBBC020 image in the Fourier basis, 2) visualize its N-term approximation. By varying N, can you get a sense of how information is distributed across Fourier coefficients?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "N = 1 ##### Try to modify this #####\n",
    "\n",
    "for idx in bbbc020_instances.keys():\n",
    "    mask = bbbc020_instances[idx]\n",
    "    \n",
    "    # Retreive instance contour\n",
    "    outline = skm.find_contours(mask,1)\n",
    "    outline = outline[0]\n",
    "    \n",
    "    # Decompose into the Fourier basis\n",
    "    (A0,C0) = pyefd.calculate_dc_coefficients(outline)\n",
    "    fourier_coefficients = pyefd.elliptic_fourier_descriptors(outline, order = N)\n",
    "    \n",
    "    # Reconstruct\n",
    "    n_term_approximation = pyefd.reconstruct_contour(fourier_coefficients, locus = (A0, C0), num_points = len(outline))\n",
    "    \n",
    "    # Visualize the result\n",
    "    fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "    \n",
    "    axes[0].scatter(outline[:,1], outline[:,0], s = 5, c = 'red')\n",
    "    axes[0].axis('equal')\n",
    "    axes[0].axis('off')\n",
    "    axes[0].set_title(\"Original shape (\"+idx+\")\")\n",
    "\n",
    "    axes[1].scatter(n_term_approximation[:,1], n_term_approximation[:,0], s = 5, c = 'red')\n",
    "    axes[1].axis('equal')\n",
    "    axes[1].axis('off')\n",
    "    axes[1].set_title(str(N)+\"-term approximation (\"+idx+\")\")\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Classify worms\n",
    "\n",
    "When *C. elegans* worms die, they rigidify into a straight rod. Even in the absence of the extra fluorescence microscopy readout, we could therefore still try to assess whether the *C. elegans* worms in the BBBC010 image are alive or not. Adapt 3.1 to use measurements from the list of region properties that you think are appropriate to capture the \"straight\" shape phenotype, and study the population of worms in the BBBC010 image using only the brightfield data. Do you manage to reach the same conclusion as in 2.2?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract shape features\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the distribution of each feature as an histogram\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Export"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A collection of measurements describing an object as it appears on image data is referred to as a *feature vector*. Feature vectors are the basic ingredient of data analysis: they allow to use the visual information contained in images in statistical studies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Group measurements\n",
    "\n",
    "Run the lines below to assemble all the intensity (2.2) and shape (3.3) measurements you have extracted on each object of the BBBC010 image into one big matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_matrix = pd.concat([bbbc010_intensity_feats, bbbc010_shape_feats], axis=1)\n",
    "\n",
    "display(feature_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Export\n",
    "\n",
    "Run the line below to save the *feature matrix* you compiled into a spreadsheet. This file now quantitatively summarizes the content of the fluorescence and brightfield image and can be used for futher analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_matrix.to_csv('data/Part_3/bbbc010_features.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}