3_qc_passed_image_display.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "from tqdm import tqdm\n",
    "import napari\n",
    "import tifffile\n",
    "import numpy as np\n",
    "from cellpose import models\n",
    "import skimage\n",
    "from skimage import measure, exposure\n",
    "from scipy.ndimage import binary_erosion, binary_dilation, label\n",
    "import pyclesperanto_prototype as cle\n",
    "from apoc import ObjectSegmenter, PixelClassifier\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import count_particles_in_nuclei\n",
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the dataset that you want to analyze below (\"microglia\" or \"astrocyte\")\n",
    "dataset = \"microglia\"\n",
    "\n",
    "if dataset == \"microglia\":\n",
    "    directory_path = Path(\"./microglia_data/\")\n",
    "    print(f\"The following dataset will be analyzed: {str(directory_path)}\")\n",
    "    \n",
    "elif dataset == \"astrocyte\":\n",
    "    directory_path = Path(\"./astrocyte_data/\")\n",
    "    print(f\"The following dataset will be analyzed: {str(directory_path)}\")\n",
    "else:\n",
    "    print(\"Please define an existing dataset\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = []\n",
    "\n",
    "# Iterate through the lsm files in the directory\n",
    "for file_path in directory_path.glob(\"*.lsm\"):\n",
    "    images.append(str(file_path))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract all the files that did pass QC and display them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the results you want to explore\n",
    "csv_path = \"./results/microglia_results_cellpdia30_sigma1_dilrad4_dnad_obj_seg_v1_gliaero6_gliathr20_dnadero2.csv\" \n",
    "mouse_id_csv_path = \"./mouse_ids_Iba1.csv\"\n",
    "\n",
    "# Read the CSV file containing the analysis results\n",
    "df = pd.read_csv(csv_path)\n",
    "\n",
    "# Read the CSV file containing the sample_id data\n",
    "df_mouse_id = pd.read_csv(mouse_id_csv_path, delimiter=\";\", encoding=\"UTF-8\")\n",
    "\n",
    "# Convert the index into a column\n",
    "df.reset_index(inplace=True)\n",
    "\n",
    "# Extract 'tissue_location'\n",
    "df['tissue_location'] = df['filename'].str.split('_40X_').str[-1]\n",
    "\n",
    "# Extract 'staining_id'\n",
    "df['staining_id'] = df['filename'].str.extract('(\\d+)_40X')\n",
    "\n",
    "# Make sure staining_id are of the same data type before merging df and df_mouse_id\n",
    "df['staining_id'] = pd.to_numeric(df['staining_id'], errors='coerce')\n",
    "df_mouse_id['staining_id'] = pd.to_numeric(df_mouse_id['staining_id'], errors='coerce')\n",
    "\n",
    "# Merge both processed_results_df and mouse_id dataframes on staining_id\n",
    "merged_df = pd.merge(df, df_mouse_id, on=\"staining_id\")\n",
    "\n",
    "# Calculate mean area of the image occupied by glia+ signal\n",
    "glia_mask_area_mean = df['%_glia+_signal'].mean() \n",
    "\n",
    "# Calculate mean area of the image occupied by dna_damage_+ signal\n",
    "dna_damage_mask_area_mean = df['%_dna_damage_signal'].mean()\n",
    "\n",
    "# Define a function to determine staining quality\n",
    "def determine_stain_quality(value, mean_value):\n",
    "    if value < (mean_value + mean_value*3):\n",
    "        return \"optimal\"\n",
    "    else:\n",
    "        return \"suboptimal\"\n",
    "\n",
    "# Check stain quality for glia and create another column storing optimal or suboptimal if qc_passed or not    \n",
    "merged_df['glia_stain_quality_auto'] = merged_df['%_glia+_signal'].apply(lambda x: determine_stain_quality(x, glia_mask_area_mean))\n",
    "\n",
    "# Check stain quality for dna_damage and create another column storing optimal or suboptimal if qc_passed or not \n",
    "merged_df['dna_damage_stain_quality_auto'] = merged_df['%_dna_damage_signal'].apply(lambda x: determine_stain_quality(x, dna_damage_mask_area_mean))\n",
    "\n",
    "# Check for both stain qualities and store True qc_passed if both are optimal\n",
    "merged_df['staining_qc_passed'] = (merged_df['glia_stain_quality_auto'] == 'optimal') & (merged_df['dna_damage_stain_quality_auto'] == 'optimal')\n",
    "\n",
    "# Group the DataFrame by 'staining_id' and check if all 'staining_qc_passed' values are True, otherwise set them all to False\n",
    "merged_df['staining_qc_passed'] = merged_df.groupby('staining_id')['staining_qc_passed'].transform('all')\n",
    "\n",
    "# Now, if all 'staining_qc_passed' values for the same 'staining_id' were True, the column will remain True; otherwise, it will be False\n",
    "\n",
    "# Remove data from images with a poor quality stain (auto QC)\n",
    "auto_filtered_df = merged_df[merged_df['staining_qc_passed'] == True]\n",
    "\n",
    "# Extract the images not passing staining_qc\n",
    "qc_passed_df = merged_df[merged_df['staining_qc_passed'] == True]\n",
    "\n",
    "print(f\"{qc_passed_df.shape[0]} stains have passed QC and will be displayed\")\n",
    "\n",
    "qc_passed_list = qc_passed_df['index'].tolist()\n",
    "\n",
    "qc_passed_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract analysis settings from results.csv name and apply them to the images that passed QC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract analysis parameters from the CSV path\n",
    "extracted_values = re.findall(r'\\d+', csv_path)\n",
    "\n",
    "# Dynamically assign the extracted values to variables\n",
    "if len(extracted_values) >= 7:\n",
    "    cellpose_nuclei_diameter = int(extracted_values[0])\n",
    "    gaussian_sigma = int(extracted_values[1])\n",
    "    dilation_radius_nuclei = int(extracted_values[2])\n",
    "    dna_damage_segmenter_version = int(extracted_values[3])\n",
    "    glia_nuclei_colocalization_erosion = int(extracted_values[4])\n",
    "    dna_damage_erosion = int(extracted_values[6])\n",
    "    if \"_glia_sem_seg_v\" in str(csv_path):\n",
    "        glia_segmenter = True\n",
    "    else:\n",
    "        glia_segmenter = False\n",
    "\n",
    "if glia_segmenter:\n",
    "    glia_segmenter_version = int(extracted_values[5])\n",
    "    glia_channel_threshold = None\n",
    "    # Dinamically adjust plot titles\n",
    "    title = f\"cellpdia{cellpose_nuclei_diameter}_sigma{gaussian_sigma}_dilrad{dilation_radius_nuclei}_dnad_obj_seg_v{dna_damage_segmenter_version}_gliaero{glia_nuclei_colocalization_erosion}_glia_sem_seg_v{glia_segmenter_version}_dnadero{dna_damage_erosion}\"\n",
    "else:\n",
    "    glia_channel_threshold = int(extracted_values[5])\n",
    "    glia_segmenter_version = None\n",
    "    # Dinamically adjust plot titles\n",
    "    title = f\"cellpdia{cellpose_nuclei_diameter}_sigma{gaussian_sigma}_dilrad{dilation_radius_nuclei}_dnad_obj_seg_v{dna_damage_segmenter_version}_gliaero{glia_nuclei_colocalization_erosion}_gliathr{glia_channel_threshold}_dnadero{dna_damage_erosion}\"\n",
    "    \n",
    "\n",
    "# Print the assigned analysis parameters\n",
    "print(f\"Cellpose nuclei diameter: {cellpose_nuclei_diameter}\")\n",
    "print(f\"Gaussian sigma: {gaussian_sigma}\")\n",
    "print(f\"Dilation radius nuclei: {dilation_radius_nuclei}\")\n",
    "print(f\"Dna damage segmenter version: {dna_damage_segmenter_version}\")\n",
    "print(f\"Glia erosion: {glia_nuclei_colocalization_erosion}\")\n",
    "print(f\"Glia threshold: {glia_channel_threshold}\")\n",
    "print(f\"Glia semantic segmentation version: {glia_segmenter_version}\")\n",
    "print(f\"DNA damage foci erosion: {dna_damage_erosion}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display in-notebook (using matplotlib) the images that passed QC using the above extracted analysis settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames_indexerror = []\n",
    "\n",
    "model = models.Cellpose(gpu=True, model_type='nuclei')\n",
    "\n",
    "for index in qc_passed_list:\n",
    "\n",
    "        image = images[index]\n",
    "\n",
    "        # Extract filename\n",
    "        file_path = Path(image)\n",
    "        filename = file_path.stem\n",
    "\n",
    "        # Read the image file\n",
    "        img = tifffile.imread(image)\n",
    "        \n",
    "        try:\n",
    "\n",
    "            # Extract the stack containing the nuclei (0), dna_damage (1) and microglia channel (2)\n",
    "            nuclei_stack = img[:, 0, :, :]\n",
    "            dna_damage_stack = img[:, 1, :, :]\n",
    "            microglia_stack = img[:, 2, :, :]\n",
    "\n",
    "            # Perform maximum intensity projections\n",
    "            nuclei_mip = np.max(nuclei_stack, axis=0)\n",
    "            dna_damage_mip = np.max(dna_damage_stack, axis=0)\n",
    "            microglia_mip = np.max(microglia_stack, axis=0)\n",
    "\n",
    "            # Create a copy of nuclei_mip\n",
    "            input_img = nuclei_mip.copy()\n",
    "\n",
    "            # Might need to perform a Gaussian-blur before\n",
    "            post_gaussian_img = skimage.filters.gaussian(input_img, sigma=gaussian_sigma)\n",
    "\n",
    "            # Apply Contrast Stretching to improve Cellpose detection of overly bright nuclei\n",
    "            p2, p98 = np.percentile(post_gaussian_img, (2, 98))\n",
    "            img_rescale = exposure.rescale_intensity(post_gaussian_img, in_range=(p2, p98))\n",
    "\n",
    "            # Predict nuclei nuclei_masks using cellpose\n",
    "            nuclei_masks, flows, styles, diams = model.eval(\n",
    "                img_rescale,\n",
    "                diameter=cellpose_nuclei_diameter,\n",
    "                channels=[0, 0],\n",
    "                net_avg=False,\n",
    "            )\n",
    "\n",
    "            # Dilate nuclei to make sure the dna_damage_foci objects sit inside\n",
    "            dilated_nuclei_masks = cle.dilate_labels(\n",
    "                nuclei_masks, radius=dilation_radius_nuclei\n",
    "            )\n",
    "            # Erode dilated_nuclei_masks to obtain separate objects upon binarization\n",
    "            eroded_nuclei_masks = cle.erode_labels(dilated_nuclei_masks, radius=1)\n",
    "\n",
    "            # Set a threshold value for the pixels in microglia channel\n",
    "            threshold = glia_channel_threshold  # Based on the microglia marker intensity across images\n",
    "\n",
    "            # Create a new array with the same shape as nuclei_masks, initialized with False\n",
    "            result_array = np.full_like(nuclei_masks, False, dtype=bool)\n",
    "\n",
    "            if glia_segmenter:\n",
    "                # Predict glia mask using the semantic_segmenter instead of thresholding\n",
    "                segmenter = PixelClassifier(\n",
    "                    opencl_filename=f\"./semantic_segmenters/microglia_segmenter_v{glia_segmenter_version}.cl\"\n",
    "                )\n",
    "                above_threshold_indices = segmenter.predict(image=microglia_mip)\n",
    "                # Extract ndarray from OCLArray\n",
    "                above_threshold_indices = cle.pull(above_threshold_indices)\n",
    "                # Transform into a boolean array\n",
    "                above_threshold_indices = above_threshold_indices > 1\n",
    "            else:\n",
    "                # Find indices where values in values_array are above the threshold\n",
    "                above_threshold_indices = microglia_mip > threshold\n",
    "\n",
    "            # Update the corresponding positions in the result_array based on the mask_array\n",
    "            nuclei_masks_bool = nuclei_masks.astype(\n",
    "                bool\n",
    "            )  # Make a boolean copy of nuclei_masks to be able to use logical operators\n",
    "            result_array[nuclei_masks_bool & above_threshold_indices] = True\n",
    "\n",
    "            # Convert the boolean array to a binary array\n",
    "            binary_result_array = result_array.astype(int)\n",
    "\n",
    "            # Now, result_array contains True where both conditions are satisfied, and False otherwise\n",
    "            # viewer.add_labels(binary_result_array, name=f\"{filename}_nuclei+glia_coloc\")\n",
    "\n",
    "            # Erode binary_result_array to get rid of small nuclei pixels colocalizing with glia branches\n",
    "            # Set the structuring element for erosion\n",
    "            structuring_element = np.ones(\n",
    "                (glia_nuclei_colocalization_erosion, glia_nuclei_colocalization_erosion),\n",
    "                dtype=bool,\n",
    "            )  # You can adjust the size and shape\n",
    "\n",
    "            # Perform binary erosion\n",
    "            eroded_array = binary_erosion(\n",
    "                binary_result_array, structure=structuring_element\n",
    "            )\n",
    "\n",
    "            # Now I want to recover just the nuclei_masks that are sitting on top of binary_results_array\n",
    "            labels, num_labels = measure.label(nuclei_masks, return_num=True)\n",
    "\n",
    "            # Create an array of indices corresponding to the True values in binary_result_array\n",
    "            true_indices = np.where(eroded_array)\n",
    "\n",
    "            # Initialize an array to store labels for each processed region\n",
    "            processed_region_labels = np.zeros_like(nuclei_masks, dtype=int)\n",
    "\n",
    "            # Label index for processed regions\n",
    "            label_index = 1\n",
    "\n",
    "            # Iterate over each connected component\n",
    "            for label_ in range(1, num_labels + 1):\n",
    "                # Extract the region corresponding to the current label\n",
    "                region = labels == label_\n",
    "\n",
    "                # Check if any True value in binary_result_array is present in the region\n",
    "                if np.any(region[true_indices]):\n",
    "                    # Assign a unique label to the processed region\n",
    "                    processed_region_labels[region] = label_index\n",
    "                    label_index += 1\n",
    "\n",
    "            # viewer.add_labels(processed_region_labels, name=f\"{filename}_glia+_nuclei_mask\")\n",
    "\n",
    "            # Dilate processed_regions to make sure the dna_damage_foci objects sit inside\n",
    "            dilated_glia_pos_nuclei = cle.dilate_labels(processed_region_labels, radius=dilation_radius_nuclei)\n",
    "            # Erode processed_regions to obtain separate objects upon binarization\n",
    "            eroded_glia_pos_nuclei = cle.erode_labels(dilated_glia_pos_nuclei, radius=1)\n",
    "\n",
    "            # Apply object segmenter from APOC\n",
    "            segmenter = ObjectSegmenter(\n",
    "                opencl_filename=f\"./object_segmenters/dna_damage_object_segmenter_v{dna_damage_segmenter_version}.cl\"\n",
    "            )\n",
    "            dna_damage_masks = segmenter.predict(image=dna_damage_mip)\n",
    "\n",
    "            # Erode dna_damage_masks to get rid of small artifacts\n",
    "            # Set the structuring element for erosion\n",
    "            structuring_element_damage = np.ones(\n",
    "                (dna_damage_erosion, dna_damage_erosion), dtype=bool\n",
    "            )  # You can adjust the size and shape\n",
    "\n",
    "            # Perform binary erosion\n",
    "            eroded_dna_damage = binary_erosion(\n",
    "                dna_damage_masks, structure=structuring_element_damage\n",
    "            )\n",
    "\n",
    "            # Perform binary dilation to fill holes\n",
    "            dilated_dna_damage = binary_dilation(\n",
    "                eroded_dna_damage, structure=structuring_element_damage\n",
    "            )\n",
    "\n",
    "            # Label the DNA damage particles\n",
    "            labeled_dna_damage, num_dna_damage = label(dilated_dna_damage)\n",
    "\n",
    "            # Label the nuclei in eroded_glia_pos_nuclei and eroded_nuclei_masks\n",
    "            labeled_glia_nuclei, num_glia_nuclei = label(eroded_glia_pos_nuclei)\n",
    "            labeled_nuclei, num_nuclei = label(eroded_nuclei_masks)\n",
    "\n",
    "            # Count particles in glia nuclei\n",
    "            particles_in_glia_nuclei = count_particles_in_nuclei(\n",
    "                labeled_dna_damage, labeled_glia_nuclei, num_glia_nuclei\n",
    "            )\n",
    "\n",
    "            # Count particles in other nuclei\n",
    "            particles_in_nuclei = count_particles_in_nuclei(\n",
    "                labeled_dna_damage, labeled_nuclei, num_nuclei\n",
    "            )\n",
    "\n",
    "            # The arrays particles_in_glia_nuclei and particles_in_nuclei now contain the counts of DNA damage particles\n",
    "            # in each nucleus of dilated_glia_pos_nuclei and dilated_nuclei_masks respectively\n",
    "            \n",
    "            # Visualize the original image and segmentation nuclei_masks\n",
    "            print(filename)\n",
    "            print(particles_in_glia_nuclei)\n",
    "            \n",
    "            plt.figure(figsize=(40, 40))\n",
    "            plt.subplot(1, 7, 1)\n",
    "            plt.imshow(nuclei_mip, cmap='gray')\n",
    "            plt.title('Input Nuclei MIP Image')\n",
    "            plt.axis(\"off\")\n",
    "\n",
    "            plt.subplot(1, 7, 2)\n",
    "            plt.imshow(eroded_nuclei_masks, cmap='viridis')\n",
    "            plt.title('Segmentation nuclei')\n",
    "            plt.axis(\"off\")\n",
    "\n",
    "            plt.subplot(1, 7, 3)\n",
    "            plt.imshow(microglia_mip, cmap='gray')\n",
    "            plt.title('Input Glia MIP Image')\n",
    "            plt.axis(\"off\")\n",
    "\n",
    "            plt.subplot(1, 7, 4)\n",
    "            plt.imshow(above_threshold_indices, cmap='viridis')\n",
    "            plt.title('Segmentation Glia')\n",
    "            plt.axis(\"off\")\n",
    "\n",
    "            plt.subplot(1, 7, 5)\n",
    "            plt.imshow(eroded_glia_pos_nuclei, cmap='viridis')\n",
    "            plt.title('Segmentation glia+ nuclei')\n",
    "            plt.axis(\"off\")\n",
    "\n",
    "            plt.subplot(1, 7, 6)\n",
    "            plt.imshow(dna_damage_mip, cmap='gray')\n",
    "            plt.title('Input DNA damage')\n",
    "            plt.axis(\"off\")\n",
    "\n",
    "            plt.subplot(1, 7, 7)\n",
    "            plt.imshow(dilated_dna_damage, cmap='gray')\n",
    "            plt.title('DNA damage mask')\n",
    "            plt.axis(\"off\")\n",
    "\n",
    "            plt.show()\n",
    "            \n",
    "        except IndexError:\n",
    "            print(f\"Filename: {filename} displays an IndexError\")\n",
    "            filenames_indexerror.append(filename)\n",
    "            pass"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "microglia",
   "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.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}