diff --git a/cuisto/config.py b/cuisto/config.py index 2d62a83..e0b7eb2 100644 --- a/cuisto/config.py +++ b/cuisto/config.py @@ -4,6 +4,7 @@ """ +import os import tomllib import warnings @@ -62,9 +63,10 @@ def __init__(self, config_file): def get_blacklist(self): """Wraps cuisto.utils.get_blacklist.""" - self.atlas["blacklist"] = utils.get_blacklist( - self.files["blacklist"], self.bg_atlas - ) + if os.path.isfile(self.files["blacklist"]): + self.atlas["blacklist"] = utils.get_blacklist( + self.files["blacklist"], self.bg_atlas + ) def get_leaves_list(self): """Wraps utils.get_leaves_list.""" diff --git a/docs/images/cuisto-showcase.svg b/docs/images/cuisto-showcase.svg index 8977cc7..a45f162 100644 --- a/docs/images/cuisto-showcase.svg +++ b/docs/images/cuisto-showcase.svg @@ -2,14 +2,14 @@ + x="-7.641221" + y="107.91206">Projections heatmap diff --git a/examples/batch_process_animals.py b/examples/batch_process_animals.py index f8dca82..1752f3e 100644 --- a/examples/batch_process_animals.py +++ b/examples/batch_process_animals.py @@ -37,8 +37,8 @@ # full path to where are the TSV files and {animal}_detections folders # dl_example = False # animals = ("animalid0", "animalid1") -# input_dir = /path/to/your/data -# config_file = /path/to/your/config.toml +# input_dir = "~/.cuisto/example" +# config_file = "~/.cuisto/config_multi.toml" # Download example data into cuisto default folder. Remove this block if using your data if dl_example: diff --git a/examples/cells_distributions.py b/examples/cells_distributions.py new file mode 100644 index 0000000..0b0d21d --- /dev/null +++ b/examples/cells_distributions.py @@ -0,0 +1,97 @@ +# This script shows how to load data exported from QuPath, compute metrics and display +# them, according to the configuration file. This is meant for a single-animal. +# +# There are some conventions that need to be met in the QuPath project so that the +# measurements are usable with `cuisto`: +# + Objects' classifications must be derived, eg. be in the form "something: else". The +# primary classification ("something") will be refered to "object_type" and the +# secondary classification ("else") to "channel" in the configuration file. +# + Only one "object_type" can be processed at once, but supports any numbers of +# channels. +# + Annotations (brain regions) must have properly formatted measurements. For punctual +# objects, it would be the count. Run the "add_regions_count.groovy" script to add +# them. The measurements names must be in the form "something: else name", for +# instance, "something: else Count". "name" is refered to "base_measurement" in the +# configuration file. +# +# The data was generated from QuPath with stardist cell detection followed by a pixel +# classifier "Classify" function on toy data. +# +# If you cloned the repository, this example should be ready to run. +# +# Find this script as a notebook in the documentation : +# https://teamncmc.github.io/cuisto/demo_notebooks/cells_distributions.html + +from pathlib import Path + +import matplotlib.pyplot as plt +import pandas as pd + +import cuisto + +# Full path to your configuration file, edited according to your need beforehand +wdir = Path(__file__).parent.parent / "resources" +config_file = wdir / "demo_config_cells.toml" + +# - Files +# animal identifier +animal = "animalid0" +# set the full path to the annotations tsv file from QuPath +annotations_file = wdir / "cells_measurements_annotations.tsv" +# set the full path to the detections tsv file from QuPath +detections_file = wdir / "cells_measurements_detections.tsv" + +# get configuration +cfg = cuisto.config.Config(config_file) +# update configuration file paths (so that this example can self-run) +cfg.files["blacklist"] = wdir / "demo_atlas_blacklist_brain.toml" +cfg.files["fusion"] = wdir / "demo_atlas_fusion_brain.toml" +cfg.files["infos"] = wdir / "demo_info_cells.toml" +cfg.get_blacklist() + +# read data +df_annotations = pd.read_csv(annotations_file, index_col="Object ID", sep="\t") +df_detections = pd.read_csv(detections_file, index_col="Object ID", sep="\t") + +# remove annotations that are not brain regions +df_annotations = df_annotations[df_annotations["Classification"] != "Region*"] +df_annotations = df_annotations[df_annotations["ROI"] != "Rectangle"] + +# convert atlas coordinates from mm to microns +df_detections[["Atlas_X", "Atlas_Y", "Atlas_Z"]] = df_detections[ + ["Atlas_X", "Atlas_Y", "Atlas_Z"] +].multiply(1000) + +# have a look +print(df_annotations.head()) +print(df_detections.head()) + +# get distributions per regions, spatial distributions and coordinates +df_regions, dfs_distributions, df_coordinates = cuisto.process.process_animal( + animal, df_annotations, df_detections, cfg, compute_distributions=True +) + +# have a look +print(df_regions.head()) +print(df_coordinates.head()) + +# plot distributions per regions +figs_regions = cuisto.display.plot_regions(df_regions, cfg) +# specify which regions to plot +# figs_regions = cuisto.display.plot_regions(df_regions, cfg, names_list=["GRN", "IRN", "MDRNv"]) + +# save as svg +# figs_regions[0].savefig(r"C:\Users\glegoc\Downloads\regions_count.svg") +# figs_regions[1].savefig(r"C:\Users\glegoc\Downloads\regions_density.svg") + +# plot 1D distributions +fig_distrib = cuisto.display.plot_1D_distributions( + dfs_distributions, cfg, df_coordinates=df_coordinates +) +# If there were several `animal` in the measurement file, it would be displayed as +# mean +/- sem instead. + +# plot heatmap (all types of cells pooled) +fig_heatmap = cuisto.display.plot_2D_distributions(df_coordinates, cfg) + +plt.show() diff --git a/examples/density_map.py b/examples/density_map.py new file mode 100644 index 0000000..2bb017f --- /dev/null +++ b/examples/density_map.py @@ -0,0 +1,167 @@ +# Draw 2D heatmaps as density isolines. +# +# This script does not actually use `cuisto` and relies only on +# brainglobe-heatmap (https://brainglobe.info/documentation/brainglobe-heatmap/index.html) +# to extract brain structures outlines. +# +# Only the detections measurements with atlas coordinates exported from QuPath are used. +# +# You need to select the range of data to be used, the regions outlines will be +# extracted at the centroid of that range. Therefore, a range that is too large will be +# misleading and irrelevant. +# +# If you cloned the repository, this example should be ready to run. +# +# Find this script as a notebook in the documentation : +# https://teamncmc.github.io/cuisto/demo_notebooks/density_map.html + +from pathlib import Path + +import brainglobe_heatmap as bgh +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns + +# path to the exported measurements from QuPath +wdir = Path(__file__).parent.parent / "resources" +filename = wdir / "cells_measurements_detections.tsv" + +# Settings +# atlas to use +atlas_name = "allen_mouse_10um" +# brain regions whose outlines will be plotted +regions = ["root", "CB", "MY", "GRN", "IRN"] +# range to include, in Allen coordinates, in microns +ap_lims = [9800, 10000] # lims : [0, 13200] for coronal +ml_lims = [5600, 5800] # lims : [0, 11400] for sagittal +dv_lims = [3900, 4100] # lims : [0, 8000] for top +# number of isolines +nlevels = 5 +# color mapping between classification and matplotlib color +palette = {"Cells: marker-": "#d8782f", "Cells: marker+": "#8ccb73"} + +df = pd.read_csv(filename, sep="\t") +print(df.head()) + +# Here we can filter out classifications we don't wan't to display. +# select objects +# df = df[df["Classification"] == "example: classification"] + +# get outline coordinates in coronal (=frontal) orientation +coords_coronal = bgh.get_structures_slice_coords( + regions, + orientation="frontal", + atlas_name=atlas_name, + position=(np.mean(ap_lims), 0, 0), +) +# get outline coordinates in sagittal orientation +coords_sagittal = bgh.get_structures_slice_coords( + regions, + orientation="sagittal", + atlas_name=atlas_name, + position=(0, 0, np.mean(ml_lims)), +) +# get outline coordinates in top (=horizontal) orientation +coords_top = bgh.get_structures_slice_coords( + regions, + orientation="horizontal", + atlas_name=atlas_name, + position=(0, np.mean(dv_lims), 0), +) + +# Coronal projection +# select objects within the rostro-caudal range +df_coronal = df[ + (df["Atlas_X"] >= ap_lims[0] / 1000) & (df["Atlas_X"] <= ap_lims[1] / 1000) +] + +plt.figure() + +for struct_name, contours in coords_coronal.items(): + for cont in contours: + plt.fill(cont[:, 0] / 1000, cont[:, 1] / 1000, lw=1, fc="none", ec="k") + +# see https://seaborn.pydata.org/generated/seaborn.kdeplot.html to customize +ax = sns.kdeplot( + df_coronal, + x="Atlas_Z", + y="Atlas_Y", + hue="Classification", + levels=nlevels, + common_norm=False, + palette=palette, +) +ax.invert_yaxis() +sns.despine(left=True, bottom=True) +plt.axis("equal") +plt.xlabel(None) +plt.ylabel(None) +plt.xticks([]) +plt.yticks([]) +plt.plot([2, 3], [8, 8], "k", linewidth=3) +plt.text(2, 7.9, "1 mm") + +# Sagittal projection +# select objects within the medio-lateral range +df_sagittal = df[ + (df["Atlas_Z"] >= ml_lims[0] / 1000) & (df["Atlas_Z"] <= ml_lims[1] / 1000) +] + +plt.figure() + +for struct_name, contours in coords_sagittal.items(): + for cont in contours: + plt.fill(cont[:, 0] / 1000, cont[:, 1] / 1000, lw=1, fc="none", ec="k") + +# see https://seaborn.pydata.org/generated/seaborn.kdeplot.html to customize +ax = sns.kdeplot( + df_sagittal, + x="Atlas_X", + y="Atlas_Y", + hue="Classification", + levels=nlevels, + common_norm=False, + palette=palette, +) +ax.invert_yaxis() +sns.despine(left=True, bottom=True) +plt.axis("equal") +plt.xlabel(None) +plt.ylabel(None) +plt.xticks([]) +plt.yticks([]) +plt.plot([2, 3], [7.1, 7.1], "k", linewidth=3) +plt.text(2, 7, "1 mm") + +# Top projection +# select objects within the dorso-ventral range +df_top = df[(df["Atlas_Y"] >= dv_lims[0] / 1000) & (df["Atlas_Y"] <= dv_lims[1] / 1000)] + +plt.figure() + +for struct_name, contours in coords_top.items(): + for cont in contours: + plt.fill(-cont[:, 0] / 1000, cont[:, 1] / 1000, lw=1, fc="none", ec="k") + +# see https://seaborn.pydata.org/generated/seaborn.kdeplot.html to customize +ax = sns.kdeplot( + df_top, + x="Atlas_Z", + y="Atlas_X", + hue="Classification", + levels=nlevels, + common_norm=False, + palette=palette, +) +ax.invert_yaxis() +sns.despine(left=True, bottom=True) +plt.axis("equal") +plt.xlabel(None) +plt.ylabel(None) +plt.xticks([]) +plt.yticks([]) +plt.plot([0.5, 1.5], [0.5, 0.5], "k", linewidth=3) +plt.text(0.5, 0.4, "1 mm") + +plt.show() diff --git a/examples/fibers_coverage.py b/examples/fibers_coverage.py new file mode 100644 index 0000000..d3c11d4 --- /dev/null +++ b/examples/fibers_coverage.py @@ -0,0 +1,78 @@ +# Plot regions coverage percentage in the spinal cord. +# +# This showcases that any brainglobe atlases should be supported. +# +# Here we're going to quantify the percentage of area of each spinal cord regions +# innervated by axons. +# +# The "area µm^2" measurement for each annotations can be created in QuPath with a pixel +# classifier, using the Measure button. +# +# We're going to consider that the "area µm^2" measurement generated by the pixel +# classifier is an object count. +# `cuisto` computes a density, which is the count in each region divided by its aera. +# Therefore, in this case, it will be actually the fraction of area covered by fibers +# in a given color. +# +# The data was generated using QuPath with a pixel classifier on toy data. +# +# If you cloned the repository, this example should be ready to run. +# +# Find this script as a notebook in the documentation : +# https://teamncmc.github.io/cuisto/demo_notebooks/fibers_coverage.html + +from pathlib import Path + +import matplotlib.pyplot as plt +import pandas as pd + +import cuisto + +# Full path to your configuration file, edited according to your need beforehand +wdir = Path(__file__).parent.parent / "resources" +config_file = wdir / "demo_config_fibers.toml" + +# - Files +# not important if only one animal +animal = "animalid1-SC" +# set the full path to the annotations tsv file from QuPath +annotations_file = wdir / "fibers_measurements_annotations.tsv" + +# get configuration +cfg = cuisto.config.Config(config_file) +# update configuration file paths (so that this example can self-run) +cfg.files["blacklist"] = wdir / "demo_atlas_blacklist_cord.toml" +cfg.files["fusion"] = wdir / "demo_atlas_fusion_cord.toml" +cfg.get_blacklist() + +# read data +df_annotations = pd.read_csv(annotations_file, index_col="Object ID", sep="\t") +df_detections = pd.DataFrame() # empty DataFrame + +# remove annotations that are not brain regions +df_annotations = df_annotations[df_annotations["Classification"] != "Region*"] +df_annotations = df_annotations[df_annotations["ROI"] != "Rectangle"] + +# have a look +print(df_annotations.head()) + +# get distributions per regions, spatial distributions and coordinates +df_regions, dfs_distributions, df_coordinates = cuisto.process.process_animal( + animal, df_annotations, df_detections, cfg, compute_distributions=False +) + +# convert the "density µm^-2" column, which is actually the coverage fraction, to a percentage +df_regions["density µm^-2"] = df_regions["density µm^-2"] * 100 + +# have a look +print(df_regions.head()) + +# plot distributions per regions +fig_regions = cuisto.display.plot_regions(df_regions, cfg) +# specify which regions to plot +# fig_regions = hq.display.plot_regions(df_regions, cfg, names_list=["Rh9", "Sr9", "8Sp"]) + +# save as svg +# fig_regions[0].savefig(r"C:\Users\glegoc\Downloads\nice_figure.svg") + +plt.show() diff --git a/examples/fibers_length_multi.py b/examples/fibers_length_multi.py new file mode 100644 index 0000000..bb85436 --- /dev/null +++ b/examples/fibers_length_multi.py @@ -0,0 +1,49 @@ +# # Fibers length in multi animals +# This example uses synthetic data to showcase how `cuisto` can be used in a pipeline. +# +# Annotations measurements should be exported from QuPath, following the required +# directory structure. +# +# Alternatively, you can merge all your CSV files yourself, one per animal, adding an +# animal ID to each table. Those can be processed with the +# `cuisto.process.process_animal()` function, in a loop, collecting the results at each +# iteration and finally concatenating the results. Finally, those can be used with +# `display` module. +# +# If you cloned the repository, this example should be ready to run. +# +# Find this script as a notebook in the documentation : +# https://teamncmc.github.io/cuisto/demo_notebooks/fibers_length_multi.html + +from pathlib import Path + +import matplotlib.pyplot as plt + +import cuisto + +# Full path to your configuration file, edited according to your need beforehand +wdir = Path(__file__).parent.parent / "resources" +config_file = wdir / "demo_config_multi.toml" + +# Files +datadir = wdir / "multi" +animals = ["mouse0", "mouse1"] + +# get configuration +cfg = cuisto.Config(config_file) +# update configuration file paths (so that this example can self-run) +cfg.files["blacklist"] = wdir / "demo_atlas_blacklist_brain.toml" +cfg.files["fusion"] = wdir / "demo_atlas_fusion_brain.toml" +cfg.get_blacklist() + +# get distributions per regions +df_regions, _, _ = cuisto.process.process_animals( + datadir, animals, cfg, compute_distributions=False +) + +# have a look +print(df_regions.head(10)) + +figs_regions = cuisto.display.plot_regions(df_regions, cfg) + +plt.show() diff --git a/examples/segment_images.py b/examples/segment_images.py index ff2c9bd..1817ddc 100644 --- a/examples/segment_images.py +++ b/examples/segment_images.py @@ -25,27 +25,68 @@ To exclude objects near the edges of an ROI, specify the path to masks stored as images with the same names as probabilities images (without their suffix). +The first block downloads toy data to showcase usage. You can drag & drop the prediction +map into QuPath, run this script and drag & drop the resulting geojson files to see the +results. Remove that block to use with your own data. + """ +import requests from tqdm import tqdm from cuisto import segmentation +### Configure the example to fetch data online. See after this block for the actual +# segmentation option + +# --- Set up the example +# Configure so that toy data is fetched online. You can safely remove this +# part when running on your own data +dl_example = True +example_url = "https://github.com/TeamNCMC/cuisto/raw/main/resources/example-seg.tar.gz" + +# Download example data into cuisto default folder. Remove this block if using your data +if dl_example: + import tarfile + import sys + from pathlib import Path + + default_destination = Path.home() / ".cuisto" + print(f"Downloading example data to {default_destination}...") + if not default_destination.exists(): + default_destination.mkdir() + + response = requests.get(example_url) + if response.ok: + tarname = default_destination / "example-seg.tar.gz" + with open(tarname, "wb") as fid: + fid.write(response.content) + with tarfile.open(tarname) as tar: + tar.extractall(tarname.parent, filter="data") + + else: + msg = ( + "Download failed. Download manually here : " + "https://github.com/TeamNCMC/cuisto/tree/main/resources" + ) + print(msg) + sys.exit() +### Remove the block above to run on your own data # --- Parameters -IMAGES_DIR = "/path/to/images" +IMAGES_DIR = default_destination / "example-seg" / "probabilities" """Full path to the images to segment.""" MASKS_DIR = "path/to/corresponding/masks" """Full path to the masks, to exclude objects near the brain edges (set to None or empty string to disable this feature).""" MASKS_EXT = "tiff" """Masks files extension.""" -SEGTYPE = "cells" +SEGTYPE = "fibers" """Type of segmentation, must match one the hardcoded keywords to associate it to 'fibers', 'points' or polygon'. See `cuisto.segmentation` doc.""" IMG_SUFFIX = "_Probabilities.tiff" """Images suffix, including extension. Masks must be the same name without the suffix.""" -ORIGINAL_PIXELSIZE = 0.4500 +ORIGINAL_PIXELSIZE = 0.5476 * 2 """Original images pixel size in microns. This is in case the pixel classifier uses a lower resolution, yielding smaller probability maps, so output objects coordinates need to be rescaled to the full size images. The pixel size is written in the "Image" @@ -55,22 +96,22 @@ { "name": "cy5", "target_channel": 0, - "proba_threshold": 0.85, - "qp_class": "Fibers: Cy5", + "proba_threshold": 0.6, + "qp_class": "Fibers: marker1", "qp_color": [164, 250, 120], }, { "name": "dsred", "target_channel": 1, - "proba_threshold": 0.65, - "qp_class": "Fibers: DsRed", + "proba_threshold": 0.6, + "qp_class": "Fibers: marker2", "qp_color": [224, 153, 18], }, { "name": "egfp", "target_channel": 2, - "proba_threshold": 0.85, - "qp_class": "Fibers: EGFP", + "proba_threshold": 0.6, + "qp_class": "Fibers: marker3", "qp_color": [135, 11, 191], }, ] diff --git a/pyproject.toml b/pyproject.toml index d2d558e..36b89d4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "cuisto" -version = "2025.5.16" +version = "2025.5.28" authors = [{ name = "Guillaume Le Goc", email = "g.legoc@posteo.org" }] description = "Quantification of objects in histological slices" readme = "README.md" diff --git a/resources/example-seg.tar.gz b/resources/example-seg.tar.gz new file mode 100644 index 0000000..1d5a40f Binary files /dev/null and b/resources/example-seg.tar.gz differ diff --git a/resources/example.tar.gz b/resources/example.tar.gz index ad58bb1..9f87c19 100644 Binary files a/resources/example.tar.gz and b/resources/example.tar.gz differ diff --git a/scripts/qupath-utils/segmentation/exportFibersAtlasCoordinates.groovy b/scripts/qupath-utils/segmentation/exportFibersAtlasCoordinates.groovy index 3bde47d..f15bcff 100644 --- a/scripts/qupath-utils/segmentation/exportFibersAtlasCoordinates.groovy +++ b/scripts/qupath-utils/segmentation/exportFibersAtlasCoordinates.groovy @@ -20,10 +20,18 @@ print('Exporting fibers coordinates...') +def pixelToAtlasTransform = + AtlasTools + .getAtlasToPixelTransform(getCurrentImageData()) + .inverse() // pixel to atlas = inverse of atlas to pixel + // Parameters def folderPrefix = 'id_' // output folder name, "segmentation" will be appended def segTag = 'fibers' // type of segmentation def atlas = 'CCFv3' // stored in the file for reference +def swapXZFlag = true // invert X and Z coordinates, should be true if CCFv3 Fiji +def mirrorLeftRight = true // invert left and right, should be true if CCfv3 Fiji +def midline = 5.70 // mediolateral midline in mm (5.70 for CCFv3) // Build file name def projectParentDir = new File(buildFilePath(PROJECT_BASE_DIR)).getParent()