Detecting spots with Spotiflow

# Libraries for plotting, reading images and tables
import matplotlib.pyplot as plt
import tifffile as tiff
import numpy as np
import pandas as pd

# We import the spotiflow library
from spotiflow.model import Spotiflow

hela = tiff.imread("data/hela_cells.tif")
img_lyso = hela[:, :, 0]

_ = plt.figure(figsize=(8, 6))
_ = plt.imshow(img_lyso, cmap="gray", vmax=np.percentile(img_lyso, 99.8))
_ = plt.axis("off")

model_spots = Spotiflow.from_pretrained("general")

points_lyso, details = model_spots.predict(img_lyso)

print(f"found {len(points_lyso)} spots")
print(points_lyso[:5])

_ = plt.figure(figsize=(8, 6))
_ = plt.imshow(img_lyso, cmap="gray", vmax=np.percentile(img_lyso, 99.8))
_ = plt.scatter(points_lyso[:, 1], points_lyso[:, 0], s=15, facecolors="none", edgecolors="orange")
_ = plt.axis("off")

# Zoom in on a 150 x 150 pixel region
y0, x0, size = 130, 50, 150

_ = plt.figure(figsize=(7, 7))
_ = plt.imshow(img_lyso, cmap="gray", vmax=np.percentile(img_lyso, 99.8))
_ = plt.scatter(points_lyso[:, 1], points_lyso[:, 0], s=80, facecolors="none", edgecolors="orange")
_ = plt.xlim(x0, x0 + size)
_ = plt.ylim(y0 + size, y0)   # y the other way round, so the image is not flipped
_ = plt.axis("off")

fig, axs = plt.subplots(1, 2, figsize=(12, 6))
_ = axs[0].imshow(img_lyso, cmap="gray", vmax=np.percentile(img_lyso, 99.8))
_ = axs[0].set_title("image")
_ = axs[1].imshow(details.heatmap, cmap="magma")
_ = axs[1].scatter(points_lyso[:, 1], points_lyso[:, 0], s=80, facecolors="none", edgecolors="cyan")
_ = axs[1].set_title("predicted heatmap and spots")
for ax in axs:
    _ = ax.set_xlim(x0, x0 + size)
    _ = ax.set_ylim(y0 + size, y0)
    _ = ax.axis("off")

print(f"lowest confidence of a detected spot: {details.prob.min():.2f}")

img = tiff.imread("data/MAX_2h_IR_Position002.tif")
img_foci = img[0, :, :]
img_nuclei = img[1, :, :]

_ = plt.figure(figsize=(7, 7))
_ = plt.imshow(img_foci, cmap="gray", vmax=50)
_ = plt.axis("off")

points, details = model_spots.predict(img_foci)

print(f"found {len(points)} spots")

_ = plt.figure(figsize=(8, 8))
_ = plt.imshow(img_foci, cmap="gray", vmax=50)
_ = plt.scatter(points[:, 1], points[:, 0], s=20, facecolors="none", edgecolors="orange")
_ = plt.axis("off")

# Zoom in on a 200 x 200 pixel region around one nucleus
y0, x0, size = 740, 20, 200

_ = plt.figure(figsize=(7, 7))
_ = plt.imshow(img_foci, cmap="gray", vmax=np.percentile(img_foci, 99.8))
_ = plt.scatter(points[:, 1], points[:, 0], s=80, facecolors="none", edgecolors="orange")
_ = plt.xlim(x0, x0 + size)
_ = plt.ylim(y0 + size, y0)   # y the other way round, so the image is not flipped
_ = plt.axis("off")

from cellpose import models

model_nuc = models.CellposeModel(model_type="nuclei")
masks, _, _ = model_nuc.eval(img_nuclei, channels=[0, 0], diameter=80)

print(f"{masks.max()} nuclei")

# The coordinates are floats, so round them to whole pixels
spot_y = np.round(points[:, 0]).astype(int)
spot_x = np.round(points[:, 1]).astype(int)

# Look up the nucleus label under every spot (0 means: not in a nucleus)
spot_label = masks[spot_y, spot_x]

print(f"{np.sum(spot_label > 0)} of {len(points)} spots are inside a nucleus")

# Count how often each label occurs: that is the number of foci per nucleus
counts = np.bincount(spot_label, minlength=masks.max() + 1)
foci_per_nucleus = counts[1:]   # drop the background (label 0)

print(f"foci per nucleus: {foci_per_nucleus}")
print(f"median: {np.median(foci_per_nucleus):.0f}")

_ = plt.hist(foci_per_nucleus, bins=15)
_ = plt.xlabel("Foci per nucleus")
_ = plt.ylabel("Number of nuclei")

def count_foci_per_nucleus(path):
    """Return the number of foci in each nucleus of a two-channel image."""

    img = tiff.imread(path)

    # Channel 0: the foci
    points, _ = model_spots.predict(img[0, :, :], verbose=False)

    # Channel 1: the nuclei
    masks, _, _ = model_nuc.eval(img[1, :, :], channels=[0, 0], diameter=80)

    # Which nucleus is each spot in, and how many spots per nucleus?
    spot_y = np.round(points[:, 0]).astype(int)
    spot_x = np.round(points[:, 1]).astype(int)
    counts = np.bincount(masks[spot_y, spot_x], minlength=masks.max() + 1)

    return counts[1:]

files = [
    "data/MAX_2h_control_Position001.tif",
    "data/MAX_2h_control_Position002.tif",
    "data/MAX_2h_IR_Position001.tif",
    "data/MAX_2h_IR_Position002.tif",
]

results = []
for path in files:
    condition = "IR" if "IR" in path else "control"
    for n_foci in count_foci_per_nucleus(path):
        results.append({"file": path, "condition": condition, "foci": n_foci})

results = pd.DataFrame(results)
print(f"{len(results)} nuclei in total")
results.head()

results.groupby("condition")["foci"].agg(["count", "mean", "median"])

# Put the control first, then the irradiated cells
conditions = ["control", "IR"]
data = [results.loc[results["condition"] == c, "foci"] for c in conditions]

_ = plt.boxplot(data, tick_labels=conditions)
_ = plt.ylabel("Foci per nucleus")