clustering-objects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# QuPath for Python Programmers 🐍\n",
    "## Graph clustering detections\n",
    "**Alan O'Callaghan, Léo Leplat, Peter Bankhead, Fiona Inglis, Laura Nicolás Sáenz**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Creating a gateway"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "First, open `Patient_test_1.ome.tiff` from the example project, and ensure you have an annotation selected (preferably a small one). Then, create a gateway as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qubalab.qupath import qupath_gateway\n",
    "from qubalab.images.qupath_server import QuPathServer\n",
    "\n",
    "gateway = qupath_gateway.create_gateway(auth_token=None, port=25333)\n",
    "qupath_server = QuPathServer(gateway) # read pixels from qupath via gateway"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "gateway.clearDetections()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Starting state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "Here, we're using a [CC0 fluorescence image](https://ftp.ebi.ac.uk/biostudies/fire/S-BIAD/463/S-BIAD463/Files/my_submission/Validation_raw/DCIS/Patient_test_2.ome.tiff) with cells detected using **InstanSeg**.\n",
    "\n",
    "This is only a 5-channel image, but this type of approach may be very useful\n",
    "for highly multiplexed images.\n",
    "\n",
    "QuPath snapshot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.imshow(qupath_gateway.create_snapshot())\n",
    "plt.axis(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Fetch InstanSeg model\n",
    "\n",
    "First, let's fetch the fluorescence InstanSeg model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from instanseg import InstanSeg\n",
    "instanseg = InstanSeg(\"fluorescence_nuclei_and_cells\", image_reader= \"tiffslide\", verbosity=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requesting pixels from an annotation\n",
    "\n",
    "Next, let's use the bounding box of the selected annotation to request pixels.\n",
    "\n",
    "This is easy with a rectangular annotation, but for odd shapes you would need make the input or filter the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "roi = gateway.getSelectedObject().getROI()\n",
    "x = roi.getBoundsX()\n",
    "y = roi.getBoundsY()\n",
    "width = roi.getBoundsWidth()\n",
    "height = roi.getBoundsHeight()\n",
    "tile = qupath_server.read_region(downsample=1, x=x, y=y, width=width, height=height)\n",
    "\n",
    "from qubalab.display.plot import plotImage\n",
    "_, ax = plt.subplots()\n",
    "plotImage(ax, tile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Running InstanSeg\n",
    "\n",
    "Now let's run InstanSeg! Here the InstanSeg Python package is handling tiling and merging for us, and provides us with two channels - one for nuclei and one for cell boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "labeled_output = instanseg.eval_medium_image(image = tile, return_image_tensor=False)\n",
    "\n",
    "from instanseg.utils.utils import show_images\n",
    "show_images(labeled_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Moving from InstanSeg to QuPath\n",
    "\n",
    "Let's just use the cell channel, for simplicity. We need to convert the objects to ensure they match the location in the original image, since InstanSeg here doesn't know that it's just working on a patch of a larger image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from instanseg.utils.utils import labels_to_features\n",
    "# channel 1 is cells, 0 is nuclei\n",
    "instanseg_objects = labels_to_features(labeled_output[:,1,:].numpy())\n",
    "\n",
    "import numpy as np\n",
    "def correct_coords(obj):\n",
    "    obj[\"geometry\"][\"coordinates\"] = [[tuple(np.add(coord, (x, y))) for coord in coord_set] for coord_set in obj[\"geometry\"][\"coordinates\"]]\n",
    "    obj[\"properties\"][\"object_type\"] = \"detection\"\n",
    "    return obj\n",
    "\n",
    "instanseg_objects = [correct_coords(obj) for obj in instanseg_objects]\n",
    "qupath_gateway.add_objects(instanseg_objects)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## InstanSeg output in QuPath\n",
    "\n",
    "Then, we can add the objects back to QuPath and take a snapshot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "qupath_gateway.refresh_qupath()\n",
    "\n",
    "plt.imshow(qupath_gateway.create_snapshot())\n",
    "plt.axis(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Detections\n",
    "We can get a list of all detections in the current image,\n",
    "and make \"intensity measurements\" on them.\n",
    "\n",
    "This is maybe a good example of why writing *too much* Java code in Python with QuBaLab is a bad idea..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "detections = gateway.getDetectionObjects()\n",
    "builder = gateway.jvm.qupath.ext.instanseg.core.DetectionMeasurer.builder()\n",
    "measurer = builder.downsample(1.0).build()\n",
    "measurer.makeMeasurements(\n",
    "    gateway.getCurrentImageData(),\n",
    "    gateway.getDetectionObjects())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "Now, let's get a list of the measurements we made."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "names = detections[0].getMeasurementList().getNames()\n",
    "\n",
    "# filter measurements for simplicity - you could also add measurements here!\n",
    "# we're just using mean intensity measurements per channel\n",
    "names = [x for x in names if x.endswith(\"Mean\")]\n",
    "names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Measurement collection\n",
    "\n",
    "We can make a pandas DataFrame of the detection measurements. This can be somewhat slow for many objects.\n",
    "\n",
    "You could also use this to assess feature importance, perform feature\n",
    "engineering..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(columns = names)\n",
    "for det in detections:\n",
    "    df = df._append(\n",
    "        {name: det.getMeasurements().get(name) for name in names},\n",
    "        ignore_index=True)\n",
    "\n",
    "# standardise the columns to avoid scale and shift effects\n",
    "normalized_df = (df - df.mean()) / df.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Running UMAP\n",
    "\n",
    "We could now run a dimensionality reduction algorithm, which might help us to visualise classifiers or other analyses in QuPath.\n",
    "\n",
    "We can also add these features back into QuPath."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import umap\n",
    "\n",
    "# run a dimensionality reduction algorithm on the measurements\n",
    "embedding = umap.UMAP().fit_transform(normalized_df)\n",
    "\n",
    "# assign back to measurement list\n",
    "for i in range(embedding.shape[0]):\n",
    "    detections[i].getMeasurementList().put(\"UMAP1\", float(embedding[i][0]))\n",
    "    detections[i].getMeasurementList().put(\"UMAP2\", float(embedding[i][1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "We could then create interactive plots in QuPath using these measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Preparing to cluster cells\n",
    "\n",
    "We can identify the KNN graph of the cells, and turn this into an adjacency matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.neighbors import kneighbors_graph\n",
    "import igraph\n",
    "\n",
    "# find KNN graph\n",
    "A = kneighbors_graph(normalized_df, 50)\n",
    "\n",
    "# convert matrix to adjacency graph\n",
    "g = igraph.Graph.Adjacency((A > 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "For large datasets you would use an approximate method, and a suitable K value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Running clustering\n",
    "\n",
    "Then, we can cluster the adjacency. We could also have done K-means or something similar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import leidenalg\n",
    "\n",
    "# partition the KNN graph based on graph modularity\n",
    "partition = leidenalg.find_partition(g, leidenalg.ModularityVertexPartition)\n",
    "\n",
    "# assign the partitions as classes to the original QuPath objects\n",
    "for i in range(embedding.shape[0]):\n",
    "    detections[i].setClassification(f\"Cluster {partition.membership[i]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Clustered cells in QuPath\n",
    "\n",
    "If we refresh QuPath and show its state, we can see that our clustering looks...\n",
    "somewhat successful?\n",
    "\n",
    "Validating the results is another (much more difficult) story that would\n",
    "require careful annotation and assessment of cluster stability.\n",
    "\n",
    "Visualising the features that drive the clustering would also be important, and may be easier in Python using for example scanpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "qupath_gateway.refresh_qupath()\n",
    "plt.imshow(qupath_gateway.create_snapshot())\n",
    "plt.axis(False)\n",
    "plt.show()"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}