initial commit

32c190f1 · Guillaume Cerutti · 32c190f1 · 32c190f1 · 32c190f1
Commit 32c190f1 authored Apr 3, 2023 by Guillaume Cerutti
--- a/.gitignore
+++ b/.gitignore
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+
+# Packages
+*.egg
+*.egg-info
+.eggs
+.Python
+*.pth
+dist/
+build/
+env/
+downloads/
+eggs/
+parts/
+bin/
+var/
+sdist/
+develop-eggs/
+.installed.cfg
+lib/
+lib64/
+
+# editors
+.idea/
+
+# Vim files
+*.swp
+*.*~
+
+# Mr Developer
+.mr.developer.cfg
+.project
+.pydevproject
+.settings
+
+
+# C extensions
+*.so
+*.dll
+*.dylib
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+
+# Compiled Object files
+*.os
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+.amlog
+.sconsign.dblite
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+
+# PyBuilder
+target/
+
+# jupyter notebooks
+.ipynb_checkpoints/
+
+# Microsoft office temporary files
+~$*.docx
+~$*.pptx
+~$*.xlsx
+
+# svn
+.svn
+# coverage
+.coverage
+
+
+.DS_Store
+
+# user custom filters
--- a/environment.yml
+++ b/environment.yml
+name: nuclei-quantification
+
+channels:
+  - mosaic
+  - conda-forge
+  - morpheme
+  - defaults
+
+dependencies:
+  - python =3.9
+  - jupyter
+  - numpy
+  - scipy
+  - pandas
+  - matplotlib
+  - pyvista
+  - ipyvtklink
+  - timagetk =3
+  - visu_core
--- a/notebook/test_nuclei_quantification.ipynb
+++ b/notebook/test_nuclei_quantification.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "6704c6e8",
+   "metadata": {},
+   "source": [
+    "# Nuclei Fie/Cas9 Quantification\n",
+    "\n",
+    "### Necessary imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7002eb02",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import scipy.ndimage as nd\n",
+    "import pandas as pd\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import pyvista as pv\n",
+    "pv.set_jupyter_backend('ipyvtklink')\n",
+    "\n",
+    "from timagetk import SpatialImage\n",
+    "from timagetk.io import imread\n",
+    "from timagetk.algorithms.peak_detection import detect_nuclei\n",
+    "from timagetk.algorithms.signal_quantification import quantify_nuclei_signal_intensity\n",
+    "\n",
+    "from visu_core.vtk.utils.image_tools import vtk_image_data_from_image\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c196d5ba",
+   "metadata": {},
+   "source": [
+    "### Data specification"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ab076632",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dirname = \"../data/\"\n",
+    "\n",
+    "filename = \"230315_GK_86_9_T1_5\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1023b95",
+   "metadata": {},
+   "source": [
+    "### Parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "86777931",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "detection_threshold = 100,\n",
+    "detection_radius_range = (0.8, 3)\n",
+    "\n",
+    "quantification_radius = 2.5"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "165d542f",
+   "metadata": {},
+   "source": [
+    "### Image loading"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8bfe671b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "image_filename = f\"{dirname}/{filename}.tif\"\n",
+    "img = imread(image_filename, channel_names=['FieGFP', 'Calco', 'Cas9RFP'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2234e546",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "image_datas = {}\n",
+    "\n",
+    "for channel_name in img.channel_names:\n",
+    "    image_data = vtk_image_data_from_image(img[channel_name].get_array(), img.voxelsize)\n",
+    "    image_data = pv.UniformGrid(image_data)\n",
+    "    image_data.point_data.GetArray(0).SetName(channel_name)\n",
+    "    image_datas[channel_name] = image_data\n",
+    "\n",
+    "plotter = pv.Plotter(notebook=True)\n",
+    "plotter.set_background('w')\n",
+    "\n",
+    "plotter.add_volume(image_datas['Cas9RFP'], cmap='Greys', clim=(0, 2000))\n",
+    "\n",
+    "plotter.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6195742e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "figure = plt.figure(figsize=(20, 10))\n",
+    "\n",
+    "extent = (\n",
+    "    -img.voxelsize[2]/2,\n",
+    "    img.extent[2]+img.voxelsize[2]/2,\n",
+    "    img.extent[1]+img.voxelsize[1]/2,\n",
+    "    -img.voxelsize[1]/2,\n",
+    ")\n",
+    "\n",
+    "figure.add_subplot(1, 2 ,1)\n",
+    "figure.gca().imshow(\n",
+    "    np.max([img['FieGFP'].get_array(), img['Cas9RFP'].get_array()], axis=0).max(axis=0),\n",
+    "    extent=extent\n",
+    ")\n",
+    "\n",
+    "figure.add_subplot(1, 2 ,2)\n",
+    "figure.gca().imshow(\n",
+    "    np.sum([img['FieGFP'].get_array(), img['Cas9RFP'].get_array()], axis=0).max(axis=0),\n",
+    "    extent=extent\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7ab91cea",
+   "metadata": {},
+   "source": [
+    "### Nuclei detection on fused channels"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0d13fd28",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nuclei_img = SpatialImage(\n",
+    "    np.sum([img['FieGFP'].get_array(), img['Cas9RFP'].get_array()], axis=0),\n",
+    "    voxelsize=img.voxelsize\n",
+    ")\n",
+    "\n",
+    "nuclei_points = detect_nuclei(\n",
+    "    nuclei_img,\n",
+    "    threshold=detection_threshold,\n",
+    "    radius_range=detection_radius_range\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "66ee5adb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "figure = plt.figure(figsize=(10, 10))\n",
+    "\n",
+    "extent = (\n",
+    "    -img.voxelsize[2]/2,\n",
+    "    img.extent[2]+img.voxelsize[2]/2,\n",
+    "    img.extent[1]+img.voxelsize[1]/2,\n",
+    "    -img.voxelsize[1]/2,\n",
+    ")\n",
+    "\n",
+    "figure.gca().imshow(\n",
+    "    nuclei_img.get_array().max(axis=0),\n",
+    "    extent=extent\n",
+    ")\n",
+    "\n",
+    "figure.gca().scatter(\n",
+    "    nuclei_points[:, 0],\n",
+    "    nuclei_points[:, 1],\n",
+    "    color='r'\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6e171d20",
+   "metadata": {},
+   "source": [
+    "### Nuclei signal quantification"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7c5f3bb1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nuclei_signals = {\n",
+    "    signal_name: quantify_nuclei_signal_intensity(\n",
+    "        img[signal_name],\n",
+    "        nuclei_points,\n",
+    "        nuclei_sigma=quantification_radius\n",
+    "    )\n",
+    "    for signal_name in ['FieGFP', 'Cas9RFP']\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3d7b5c4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "figure = plt.figure(figsize=(20, 10))\n",
+    "\n",
+    "extent = (\n",
+    "    -img.voxelsize[2]/2,\n",
+    "    img.extent[2]+img.voxelsize[2]/2,\n",
+    "    img.extent[1]+img.voxelsize[1]/2,\n",
+    "    -img.voxelsize[1]/2,\n",
+    ")\n",
+    "\n",
+    "for i, signal_name in enumerate(['FieGFP', 'Cas9RFP']):\n",
+    "    figure.add_subplot(1, 2, i+1)\n",
+    "    figure.gca().imshow(\n",
+    "        img[signal_name].get_array().max(axis=0),\n",
+    "        extent=extent,\n",
+    "    )\n",
+    "\n",
+    "    figure.gca().scatter(\n",
+    "        nuclei_points[:, 0],\n",
+    "        nuclei_points[:, 1],\n",
+    "        c=nuclei_signals[signal_name],\n",
+    "        cmap='Reds'\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e7b756f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nuclei_mask = nuclei_points[:, 1]>200\n",
+    "#nuclei_mask = np.ones_like(nuclei_points[:, 1]).astype(bool)\n",
+    "\n",
+    "figure = plt.figure(figsize=(10, 10))\n",
+    "\n",
+    "figure.gca().scatter(\n",
+    "    nuclei_signals['Cas9RFP'][nuclei_mask],\n",
+    "    nuclei_signals['FieGFP'][nuclei_mask],\n",
+    "    c=nuclei_points[:, 1][nuclei_mask],\n",
+    "    cmap=\"Reds\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "95a97827",
+   "metadata": {},
+   "source": [
+    "### Data export"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0bea4d68",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nuclei_data = {\n",
+    "    'id': range(len(nuclei_points)),\n",
+    "}\n",
+    "nuclei_data.update({\n",
+    "    f'center_{dim}': nuclei_points[:, i]\n",
+    "    for i, dim in enumerate('xyz')\n",
+    "})\n",
+    "nuclei_data.update(nuclei_signals)\n",
+    "\n",
+    "nuclei_df = pd.DataFrame(nuclei_data)\n",
+    "nuclei_df.head()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "487153ec",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nuclei_filename = f\"{dirname}/{filename}_nuclei_data.csv\"\n",
+    "nuclei_df.to_csv(nuclei_filename, index=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e7b650b8",
+   "metadata": {},
+   "source": [
+    "### Individual nuclei slices"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "95d74f9c",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "nucleus_window = 30\n",
+    "\n",
+    "nucleus_extent = (\n",
+    "    (-nucleus_window-0.5)*img.voxelsize[2],\n",
+    "    (nucleus_window+0.5)*img.voxelsize[2],\n",
+    "    (nucleus_window+0.5)*img.voxelsize[1],\n",
+    "    (-nucleus_window-0.5)*img.voxelsize[1],\n",
+    ")\n",
+    "\n",
+    "for nucleus_id in range(len(nuclei_points)):\n",
+    "    nucleus_coords = (nuclei_points[nucleus_id][::-1]/img.voxelsize).astype(int)\n",
+    "\n",
+    "    figure = plt.figure(figsize=(20, 10))\n",
+    "    \n",
+    "    for i, signal_name in enumerate(['FieGFP', 'Cas9RFP']):\n",
+    "        figure.add_subplot(1, 2, i+1)\n",
+    "\n",
+    "        nucleus_img = img[signal_name].get_array()[nucleus_coords[0], \n",
+    "                          nucleus_coords[1]-40:nucleus_coords[1]+40,\n",
+    "                          nucleus_coords[2]-40:nucleus_coords[2]+40]\n",
+    "\n",
+    "        figure.gca().imshow(\n",
+    "            nucleus_img,\n",
+    "            extent=nucleus_extent,\n",
+    "        )\n",
+    "        \n",
+    "        figure.gca().set_title(signal_name, size=18)\n",
+    "    figure.suptitle(f\"Nucleus {nucleus_id}\", size=24)\n",
+    "    figure.tight_layout()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e32bc393",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.9.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
+%% Cell type:markdown id:6704c6e8 tags:
+
+# Nuclei Fie/Cas9 Quantification
+
+### Necessary imports
+
+%% Cell type:code id:7002eb02 tags:
+
+``` python
+import numpy as np
+import scipy.ndimage as nd
+import pandas as pd
+
+import matplotlib.pyplot as plt
+import pyvista as pv
+pv.set_jupyter_backend('ipyvtklink')
+
+from timagetk import SpatialImage
+from timagetk.io import imread
+from timagetk.algorithms.peak_detection import detect_nuclei
+from timagetk.algorithms.signal_quantification import quantify_nuclei_signal_intensity
+
+from visu_core.vtk.utils.image_tools import vtk_image_data_from_image
+```
+
+%% Cell type:markdown id:c196d5ba tags:
+
+### Data specification
+
+%% Cell type:code id:ab076632 tags:
+
+``` python
+dirname = "../data/"
+
+filename = "230315_GK_86_9_T1_5"
+```
+
+%% Cell type:markdown id:a1023b95 tags:
+
+### Parameters
+
+%% Cell type:code id:86777931 tags:
+
+``` python
+detection_threshold = 100,
+detection_radius_range = (0.8, 3)
+
+quantification_radius = 2.5
+```
+
+%% Cell type:markdown id:165d542f tags:
+
+### Image loading
+
+%% Cell type:code id:8bfe671b tags:
+
+``` python
+image_filename = f"{dirname}/{filename}.tif"
+img = imread(image_filename, channel_names=['FieGFP', 'Calco', 'Cas9RFP'])
+```
+
+%% Cell type:code id:2234e546 tags:
+
+``` python
+image_datas = {}
+
+for channel_name in img.channel_names:
+    image_data = vtk_image_data_from_image(img[channel_name].get_array(), img.voxelsize)
+    image_data = pv.UniformGrid(image_data)
+    image_data.point_data.GetArray(0).SetName(channel_name)
+    image_datas[channel_name] = image_data
+
+plotter = pv.Plotter(notebook=True)
+plotter.set_background('w')
+
+plotter.add_volume(image_datas['Cas9RFP'], cmap='Greys', clim=(0, 2000))
+
+plotter.show()
+```
+
+%% Cell type:code id:6195742e tags:
+
+``` python
+figure = plt.figure(figsize=(20, 10))
+
+extent = (
+    -img.voxelsize[2]/2,
+    img.extent[2]+img.voxelsize[2]/2,
+    img.extent[1]+img.voxelsize[1]/2,
+    -img.voxelsize[1]/2,
+)
+
+figure.add_subplot(1, 2 ,1)
+figure.gca().imshow(
+    np.max([img['FieGFP'].get_array(), img['Cas9RFP'].get_array()], axis=0).max(axis=0),
+    extent=extent
+)
+
+figure.add_subplot(1, 2 ,2)
+figure.gca().imshow(
+    np.sum([img['FieGFP'].get_array(), img['Cas9RFP'].get_array()], axis=0).max(axis=0),
+    extent=extent
+)
+```
+
+%% Cell type:markdown id:7ab91cea tags:
+
+### Nuclei detection on fused channels
+
+%% Cell type:code id:0d13fd28 tags:
+
+``` python
+nuclei_img = SpatialImage(
+    np.sum([img['FieGFP'].get_array(), img['Cas9RFP'].get_array()], axis=0),
+    voxelsize=img.voxelsize
+)
+
+nuclei_points = detect_nuclei(
+    nuclei_img,
+    threshold=detection_threshold,
+    radius_range=detection_radius_range
+)
+```
+
+%% Cell type:code id:66ee5adb tags:
+
+``` python
+figure = plt.figure(figsize=(10, 10))
+
+extent = (
+    -img.voxelsize[2]/2,
+    img.extent[2]+img.voxelsize[2]/2,
+    img.extent[1]+img.voxelsize[1]/2,
+    -img.voxelsize[1]/2,
+)
+
+figure.gca().imshow(
+    nuclei_img.get_array().max(axis=0),
+    extent=extent
+)
+
+figure.gca().scatter(
+    nuclei_points[:, 0],
+    nuclei_points[:, 1],
+    color='r'
+)
+```
+
+%% Cell type:markdown id:6e171d20 tags:
+
+### Nuclei signal quantification
+
+%% Cell type:code id:7c5f3bb1 tags:
+
+``` python
+nuclei_signals = {
+    signal_name: quantify_nuclei_signal_intensity(
+        img[signal_name],
+        nuclei_points,
+        nuclei_sigma=quantification_radius
+    )
+    for signal_name in ['FieGFP', 'Cas9RFP']
+}
+```
+
+%% Cell type:code id:a3d7b5c4 tags:
+
+``` python
+import matplotlib.pyplot as plt
+
+figure = plt.figure(figsize=(20, 10))
+
+extent = (
+    -img.voxelsize[2]/2,
+    img.extent[2]+img.voxelsize[2]/2,
+    img.extent[1]+img.voxelsize[1]/2,
+    -img.voxelsize[1]/2,
+)
+
+for i, signal_name in enumerate(['FieGFP', 'Cas9RFP']):
+    figure.add_subplot(1, 2, i+1)
+    figure.gca().imshow(
+        img[signal_name].get_array().max(axis=0),
+        extent=extent,
+    )
+
+    figure.gca().scatter(
+        nuclei_points[:, 0],
+        nuclei_points[:, 1],
+        c=nuclei_signals[signal_name],
+        cmap='Reds'
+    )
+```
+
+%% Cell type:code id:6e7b756f tags:
+
+``` python
+nuclei_mask = nuclei_points[:, 1]>200
+#nuclei_mask = np.ones_like(nuclei_points[:, 1]).astype(bool)
+
+figure = plt.figure(figsize=(10, 10))
+
+figure.gca().scatter(
+    nuclei_signals['Cas9RFP'][nuclei_mask],
+    nuclei_signals['FieGFP'][nuclei_mask],
+    c=nuclei_points[:, 1][nuclei_mask],
+    cmap="Reds",
+)
+```
+
+%% Cell type:markdown id:95a97827 tags:
+
+### Data export
+
+%% Cell type:code id:0bea4d68 tags:
+
+``` python
+nuclei_data = {
+    'id': range(len(nuclei_points)),
+}
+nuclei_data.update({
+    f'center_{dim}': nuclei_points[:, i]
+    for i, dim in enumerate('xyz')
+})
+nuclei_data.update(nuclei_signals)
+
+nuclei_df = pd.DataFrame(nuclei_data)
+nuclei_df.head()
+```
+
+%% Cell type:code id:487153ec tags:
+
+``` python
+nuclei_filename = f"{dirname}/{filename}_nuclei_data.csv"
+nuclei_df.to_csv(nuclei_filename, index=False)
+```
+
+%% Cell type:markdown id:e7b650b8 tags:
+
+### Individual nuclei slices
+
+%% Cell type:code id:95d74f9c tags:
+
+``` python
+nucleus_window = 30
+
+nucleus_extent = (
+    (-nucleus_window-0.5)*img.voxelsize[2],
+    (nucleus_window+0.5)*img.voxelsize[2],
+    (nucleus_window+0.5)*img.voxelsize[1],
+    (-nucleus_window-0.5)*img.voxelsize[1],
+)
+
+for nucleus_id in range(len(nuclei_points)):
+    nucleus_coords = (nuclei_points[nucleus_id][::-1]/img.voxelsize).astype(int)
+
+    figure = plt.figure(figsize=(20, 10))
+
+    for i, signal_name in enumerate(['FieGFP', 'Cas9RFP']):
+        figure.add_subplot(1, 2, i+1)
+
+        nucleus_img = img[signal_name].get_array()[nucleus_coords[0],
+                          nucleus_coords[1]-40:nucleus_coords[1]+40,
+                          nucleus_coords[2]-40:nucleus_coords[2]+40]
+
+        figure.gca().imshow(
+            nucleus_img,
+            extent=nucleus_extent,
+        )
+
+        figure.gca().set_title(signal_name, size=18)
+    figure.suptitle(f"Nucleus {nucleus_id}", size=24)
+    figure.tight_layout()
+```
+
+%% Cell type:code id:e32bc393 tags:
+
+``` python
+```