Notebook 1: Does a protein cluster?

Notebook 1: Does a protein cluster?#

# /// script
# requires-python = ">=3.10"
# dependencies = [
#     "matplotlib",
#     "numpy",
#     "scikit-image",
#     "scipy",
#     "tifffile",
#     "imagecodecs",
#     "pandas",
#     "seaborn",
#     "bobiac_tools @ git+https://github.com/bobiac/bobiac-tools.git"
# ]
# ///

Overview#

In this notebook, we will learn how to assess whether individual proteins are clustering.

Dataset#

Our data consists of fluorescence images with five channels:

Channel	Staining	Staining Description
0	Hoechst	Nuclear stain
1	Phalloidin	F-actin stain; whole cell marker
2	Protein A	Spot pattern
3	Protein B	Spot pattern
4	Protein C	Spot pattern

The data for this exercise is the same as used for the classic segmentation, deep learning segmentation and spot detection exercises and can be downloaded here. From the previous segmentation exercises you should have a collection of labelled masks for the cells, and coordinates for spots. If not you can download them from here.

0. Importing libraries#

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import skimage
import tifffile
from bobiac_tools import overlay_labels, ripleys_l
from scipy.spatial import KDTree

2. Load image, mask and spot table#

First, we load all the data that we will need for this notebook.

path_image = "../../../_static/images/quant/05_spatial_stats/images/F01_202_5ch.tif"
path_mask = "../../../_static/images/quant/05_spatial_stats/masks/F01_202_cell_labels.tif"
path_spots = "../../../_static/images/quant/05_spatial_stats/spots/F01_204_5ch_points.csv"

image = tifffile.imread(path_image)
mask = tifffile.imread(path_mask)
df_spots = pd.read_csv(path_spots)
df_spots = df_spots.rename({'axis-0':'channel', 'axis-1':'y','axis-2':'x'},axis=1)

mask = skimage.segmentation.clear_border(mask)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[3], line 7
path_spots = "../../../_static/images/quant/05_spatial_stats/spots/F01_204_5ch_points.csv"

image = tifffile.imread(path_image)
mask = tifffile.imread(path_mask)
----> 7 df_spots = pd.read_csv(path_spots)
df_spots = df_spots.rename({'axis-0':'channel', 'axis-1':'y','axis-2':'x'},axis=1)

mask = skimage.segmentation.clear_border(mask)

File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:873, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, skip_blank_lines, parse_dates, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, low_memory, memory_map, float_precision, storage_options, dtype_backend)
kwds_defaults = _refine_defaults_read(
   dialect,
   delimiter,
   (...)    869     dtype_backend=dtype_backend,
)
kwds.update(kwds_defaults)
--> 873 return _read(filepath_or_buffer, kwds)

File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:300, in _read(filepath_or_buffer, kwds)
_validate_names(kwds.get("names", None))
# Create the parser.
--> 300 parser = TextFileReader(filepath_or_buffer, **kwds)
if chunksize or iterator:
   return parser

File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1645, in TextFileReader.__init__(self, f, engine, **kwds)
   self.options["has_index_names"] = kwds["has_index_names"]
self.handles: IOHandles | None = None
-> 1645 self._engine = self._make_engine(f, self.engine)

File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1904, in TextFileReader._make_engine(self, f, engine)
   if "b" not in mode:
       mode += "b"
-> 1904 self.handles = get_handle(
   f,
   mode,
   encoding=self.options.get("encoding", None),
   compression=self.options.get("compression", None),
   memory_map=self.options.get("memory_map", False),
   is_text=is_text,
   errors=self.options.get("encoding_errors", "strict"),
   storage_options=self.options.get("storage_options", None),
)
assert self.handles is not None
f = self.handles.handle

File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/pandas/io/common.py:930, in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
elif isinstance(handle, str):
   # Check whether the filename is to be opened in binary mode.
   # Binary mode does not support 'encoding' and 'newline'.
   if ioargs.encoding and "b" not in ioargs.mode:
       # Encoding
--> 930         handle = open(
           handle,
           ioargs.mode,
           encoding=ioargs.encoding,
           errors=errors,
           newline="",
       )
   else:
       # Binary mode
       handle = open(handle, ioargs.mode)

FileNotFoundError: [Errno 2] No such file or directory: '../../../_static/images/quant/05_spatial_stats/spots/F01_204_5ch_points.csv'

Visually inspect the images, masks and spots#

To be sure that the data we are working with is corrects, we visualise it here.

overlay_labels(image[2:, :, :])

../../../_images/d77e376639dfad93eafcb79c5caafda25d6c8c24ca99b93d5ac304e0c0273657.png

(<Figure size 800x800 with 1 Axes>, <Axes: >)

overlay_labels(label_mask=mask, coordinates=df_spots)

../../../_images/1431eec3a3f28e8d094e47291fc187732b4d1d0b83dfc8d222bef3d32640e714.png

(<Figure size 800x800 with 1 Axes>,
 <Axes: title={'center': 'Segmentation overlay'}>)

Using the focus attribute in overlay we can zoom in to view single objects.

overlay_labels(
    image=image[2:, :, :], label_mask=mask, coordinates=df_spots, focus_object=11
)

../../../_images/cd3fbda63a0a442c7b5e9446f6c3e37a685995ec78f75421800a09cfbc226596.png

(<Figure size 800x586.081 with 1 Axes>,
 <Axes: title={'center': 'Segmentation overlay'}>)

3. Make An Analysis Plan#

Question#

Are the spots clustered?

Approach#

To answer this question, we need to first define a reference for comparison. A common comparison is complete spatial randomness (CSR) (i.e., a fully random distribution of points). We can then define a metric that measures if a collection of points are clustered, and then use this metric to compare our dataset collection of points to a collection of random points.

Metric for Clustering#

Here, we use the average distance of each point to its nearest neighbour as the metric for clustering.

Workflow#

1. Generate a reference for comparison#

Randomly distributed set of points inside each cell’s cytoplasm. We will call this the random dataset, which we will compare to our protein dataset.

2. Measure the nearest neighbor distance#

For each spot, identify the closest spot within a cell and measure the distance.

3. Calculate the average nearest neighbour distance#

For each cell, calculate the average nearest neighbor distance for a spot.

4. Calculate the Clark-Evans Index#

The Clark-Evans Index is a ratio between the average nearest neighbour distances for the protein dataset and the random dataset. This index represents how random or clustered the dot protein distribution is. - Ratio = 1: The distribution is random - Ratio > 1: The distribution is clusered

5. Calculate the p-value of the Clark-Evans Index#

Calculate the p-value using a comparison of multiple random patterns.

To start, we focus on Protein A (Channel 2).

4. Assemble the Protein Dataset and the Random Dataset#

Protein Dataset: Map spots to cells and count spots per cell#

We need to know how many points are within each cell to place the same number of random points in the next step. If we would generate more or less random points, we would artificially measure clustering or dispersion.

First we generate a dataframe for the cells by using regionprops_table.

df_cell = pd.DataFrame(
    skimage.measure.regionprops_table(mask, properties=["area", "label"])
)
print(df_cell)

      area  label
 5634.0      1
 4728.0      2
 7677.0      3
 6824.0      4
 9638.0      5
..     ...    ...
6577.0     80
8947.0     81
2622.0     82
1669.0     84
4207.0     86

[76 rows x 2 columns]

Then we count the spots per cell (like in Exercise B: Nested Objects).

# map the cell label to the spots
df_spots["label_cell"] = mask[df_spots["y"].astype(int), df_spots["x"].astype(int)]

# delete spots in the background
df_spots = df_spots.loc[df_spots["label_cell"] != 0]

# count number of spots per cell label
mapping_spots = (
    df_spots.loc[df_spots["channel"] == 2]["label_cell"].value_counts().reset_index()
)

# add the spot count to the cell dataframe
df_cell = df_cell.merge(
    mapping_spots, left_on="label", right_on="label_cell", how="inner"
).drop(columns="label_cell")
print(df_spots)
print(df_cell)

       channel     y    x  label_cell
          0    14  672           1
          0    31  700           1
          0    64  672           1
          0     2  686           1
          0     1  678           1
...        ...   ...  ...         ...
      2  1000  768          86
      2  1020  709          86
      2   980  799          86
      2  1004  779          86
      2  1029  764          86

[10657 rows x 4 columns]
      area  label  count
 5634.0      1     37
 4728.0      2     34
 7677.0      3     36
 6824.0      4     44
 9638.0      5     56
..     ...    ...    ...
6577.0     80     41
8947.0     81     37
2622.0     82     18
1669.0     84     12
4207.0     86     35

[76 rows x 3 columns]

overlay_labels(label_mask=mask, df=df_cell, id_col="label", measurement_col="count")

../../../_images/262ff13188eb53e65aedd5ce11d1447e4e001ab4f638aa1a9c3cac22843338cf.png

(<Figure size 800x800 with 2 Axes>, <Axes: title={'center': 'count'}>)

Random Dataset: Generate random points#

Let’s generate our first set of random points. To make it more transparent, let’s focus on a single cell first.

cell_id = 11
spots = df_spots.loc[(df_spots["channel"] == 0) & (df_spots["label_cell"] == cell_id)][
    ["y", "x"]
].to_numpy()

overlay_labels(label_mask=mask, coordinates=spots, focus_object=cell_id)

../../../_images/dabd0c6303f604975c2fd63a15845948e49dbea9fad8124d81337a6edf781b87.png

(<Figure size 800x586.081 with 1 Axes>,
 <Axes: title={'center': 'Segmentation overlay'}>)

We can use a pre-written function to generate the random points. In short, what it does is:

It calculates the bounding box of the cell (the box that touches the edge of the cell at every side).
It draws random points inside that box and tests if they fall within the cell.
- If not, the point is discarded.
- If it is within the cell it is kept. This sequence is repeated until the cell contains the specified number of points.

This procedure is necessary since cells have irregular shapes.

# or use:
# from bobiac_tools import random_points_in_mask


def random_points_in_mask(mask: np.ndarray, cell_label: int, n: int) -> np.ndarray:
    """
    Generate n random points uniformly distributed within a labeled cell region.

    Points are placed continuously using rejection sampling:
    candidate points are drawn uniformly from the bounding box of the cell and
    accepted only if they fall within the cell mask. This correctly handles
    irregular cell shapes and serves as a null model for Complete Spatial
    Randomness (CSR).

    Parameters
    ----------
    mask : 2D numpy array
        Label mask where each cell is identified by a unique integer ID
        and background is 0.
    cell_label : int
        ID of the cell to sample within.
    n : int
        Number of random points to generate.

    Returns
    -------
    points : (n, 2) numpy array
        Random point coordinates in (row, col) / (y, x) order.
    """
    ys, xs = np.where(mask == cell_label)
    y_min, y_max = ys.min(), ys.max()
    x_min, x_max = xs.min(), xs.max()
    points = []
    while len(points) < n:
        y = np.random.uniform(y_min, y_max)
        x = np.random.uniform(x_min, x_max)
        if mask[int(round(y)), int(round(x))] == cell_label:
            points.append([y, x])
    return np.array(points)

Look up how many points are in our cell of interest.

n_points = df_cell.loc[df_cell["label"] == cell_id]["count"].values[0]
print(n_points)

Create the same number of random points inside that cell.

rnd_points = random_points_in_mask(mask, cell_label=cell_id, n=n_points)
rnd_points.shape

(74, 2)

overlay_labels(label_mask=mask, coordinates=[spots, rnd_points], focus_object=cell_id)

../../../_images/4a1984fe6f807fa13795a8c73fce8ca1b4fd974ba0151287cc454498a0b1cd9a.png

(<Figure size 800x586.081 with 1 Axes>,
 <Axes: title={'center': 'Segmentation overlay'}>)

5. Calculate Nearest Neighbor Distances#

Identifying Neighbors: KD Trees#

Now that we have two sets of points, our protein dataset and the random dataset as a control, we can find the nearest neighbour of every point and calculate the distance. To find neighbors, we will use a KD tree, a spatial partitioning algorithm that identities neighbours very effectively.

The syntax for the KD tree function is a bit strange at first:

KDtree(points_A).query(points_B, k=k)

KDTree(points_A) builds the KD tree data structure.
.query(points_B, k=k) computes the distance to the nearest neighbors. It measures the distance from every point in points_B to the k nearest neighbours in points_A. If the two point sets are the same (like in our case), the first distance will be 0, since it is the distance from each point to itself.

# the query below returns an array of distances and an array of indices. We don't need the indicis so we assign it to the variable `_`.
dist_spots, _ = KDTree(spots).query(spots, k=2)
print(dist_spots)

[[ 0.         11.40175425]
 [ 0.          5.09901951]
 [ 0.          8.24621125]
 [ 0.          9.8488578 ]
 [ 0.          6.        ]
 [ 0.          8.06225775]
 [ 0.          5.83095189]
 [ 0.          6.32455532]
 [ 0.          4.47213595]
 [ 0.         13.41640786]
 [ 0.          3.16227766]
 [ 0.          9.89949494]
 [ 0.         10.81665383]
 [ 0.         14.86606875]
 [ 0.          6.08276253]
 [ 0.          7.        ]
 [ 0.          8.60232527]
 [ 0.          5.09901951]
 [ 0.          2.82842712]
 [ 0.          9.8488578 ]
 [ 0.          8.94427191]
 [ 0.          3.16227766]
 [ 0.          9.48683298]
 [ 0.         13.15294644]
 [ 0.         15.5241747 ]
 [ 0.         11.18033989]
 [ 0.          8.24621125]
 [ 0.          9.8488578 ]
 [ 0.         11.40175425]
 [ 0.          7.28010989]
 [ 0.          6.08276253]
 [ 0.          8.24621125]
 [ 0.         17.        ]
 [ 0.          4.47213595]
 [ 0.          3.60555128]
 [ 0.          3.16227766]
 [ 0.          5.83095189]
 [ 0.          6.08276253]
 [ 0.          5.        ]
 [ 0.          2.82842712]
 [ 0.          4.47213595]
 [ 0.          5.09901951]
 [ 0.         10.81665383]
 [ 0.          5.09901951]
 [ 0.          8.54400375]
 [ 0.          8.        ]
 [ 0.          8.06225775]
 [ 0.          3.60555128]
 [ 0.         15.03329638]
 [ 0.          8.48528137]
 [ 0.         11.66190379]
 [ 0.          4.12310563]
 [ 0.          2.82842712]
 [ 0.          3.16227766]
 [ 0.          5.        ]
 [ 0.          8.06225775]
 [ 0.          4.47213595]
 [ 0.          5.09901951]
 [ 0.          5.        ]
 [ 0.          2.82842712]
 [ 0.         15.        ]
 [ 0.          2.82842712]
 [ 0.         13.41640786]
 [ 0.          8.06225775]
 [ 0.         11.18033989]
 [ 0.          4.12310563]
 [ 0.          5.09901951]
 [ 0.          3.        ]
 [ 0.          9.21954446]
 [ 0.          8.60232527]
 [ 0.         13.60147051]
 [ 0.          9.8488578 ]
 [ 0.          5.        ]
 [ 0.          5.        ]
 [ 0.          3.16227766]
 [ 0.          3.16227766]
 [ 0.          4.47213595]
 [ 0.          8.06225775]
 [ 0.          6.32455532]
 [ 0.         11.18033989]
 [ 0.         10.63014581]
 [ 0.          5.        ]
 [ 0.          8.54400375]
 [ 0.          9.48683298]
 [ 0.          2.82842712]
 [ 0.          3.60555128]
 [ 0.          3.        ]]

We do the same for the random spots.

# the query below returns an array of distances and an array of indices. We don't need the indicis so we assign it to the variable `_`.
dist_rnd, _ = KDTree(rnd_points).query(rnd_points, k=2)
print(dist_rnd)

[[ 0.         15.3958515 ]
 [ 0.          7.20580281]
 [ 0.          3.96582018]
 [ 0.         19.44095057]
 [ 0.          4.13343581]
 [ 0.          5.08668562]
 [ 0.          4.13343581]
 [ 0.          6.16258696]
 [ 0.          7.29186544]
 [ 0.         12.63603643]
 [ 0.          3.33469311]
 [ 0.         12.03641756]
 [ 0.          2.67812992]
 [ 0.         11.60068589]
 [ 0.          7.29186544]
 [ 0.          4.53051891]
 [ 0.         12.55383452]
 [ 0.         10.98919797]
 [ 0.         12.24902428]
 [ 0.         10.98919797]
 [ 0.          2.29882827]
 [ 0.          5.02544218]
 [ 0.          6.63401915]
 [ 0.          7.97604165]
 [ 0.          6.15636203]
 [ 0.         13.07662977]
 [ 0.          2.67812992]
 [ 0.          7.7644232 ]
 [ 0.         11.86730545]
 [ 0.          5.68846943]
 [ 0.          5.36514523]
 [ 0.         10.39214696]
 [ 0.          5.150342  ]
 [ 0.          5.02544218]
 [ 0.          7.7546019 ]
 [ 0.          8.54076906]
 [ 0.          3.96582018]
 [ 0.          3.57326492]
 [ 0.          6.46116824]
 [ 0.          3.33469311]
 [ 0.          8.37671678]
 [ 0.          2.29882827]
 [ 0.         12.41018317]
 [ 0.         16.86317781]
 [ 0.         13.22748147]
 [ 0.          7.97604165]
 [ 0.         11.60068589]
 [ 0.          2.72435881]
 [ 0.          2.72435881]
 [ 0.          6.45302973]
 [ 0.         11.57169912]
 [ 0.         12.83321148]
 [ 0.         14.93836663]
 [ 0.         11.88885135]
 [ 0.         14.21393896]
 [ 0.          5.68846943]
 [ 0.          4.85144427]
 [ 0.          5.150342  ]
 [ 0.          7.7546019 ]
 [ 0.          7.56202746]
 [ 0.          4.53051891]
 [ 0.          8.43900709]
 [ 0.          6.46116824]
 [ 0.         12.83321148]
 [ 0.          9.90728851]
 [ 0.         19.22655019]
 [ 0.         14.93836663]
 [ 0.          3.57326492]
 [ 0.         15.19319633]
 [ 0.         14.95995653]
 [ 0.          4.85144427]
 [ 0.          5.08668562]
 [ 0.         15.05848189]
 [ 0.         17.10888677]]

6. Calculate the Clark-Evans Index for the protein and random datasets#

Now, let’s calculate the mean nearest neighbour distances for both spot patterns and compare them using the Clark-Evans Index.

spot_dist_avg = dist_spots[:, 1].mean()
rnd_dist_avg = dist_rnd[:, 1].mean()
print("Mean NN distance spots:", spot_dist_avg)
print("Mean NN distance random:", rnd_dist_avg)
print("Clark-Evans index (ratio):", spot_dist_avg / rnd_dist_avg)

Mean NN distance spots: 7.360255526882802
Mean NN distance random: 8.590688160680662
Clark-Evans index (ratio): 0.8567713539609649

7. Repeat Steps 4-6 with Multiple Random Patterns#

Concept#

A single random pattern as a control dataset is not a reliable comparison because a random pattern could be clustered by chance. To account for this possibility, we repeat this procedure with additional random patterns.

n_repeats = 10

list_rnd_dist = []  # this list holds the average NN distances of the random points for each run

for _i in range(n_repeats):
    rnd_points = random_points_in_mask(mask, cell_label=cell_id, n=n_points)
    dist_rnd, _ = KDTree(rnd_points).query(rnd_points, k=2)
    rnd_dist_avg = dist_rnd[:, 1].mean()
    list_rnd_dist.append(rnd_dist_avg)

rnd_dist_avg = np.array(list_rnd_dist).mean()
clark_evans = spot_dist_avg / rnd_dist_avg

print("Mean NN distance spots:", spot_dist_avg)
print(f"Mean NN distance random ({n_repeats} runs):", rnd_dist_avg)
print(f"Clark-Evans index ({n_repeats} runs):", clark_evans)

Mean NN distance spots: 7.360255526882802
Mean NN distance random (10 runs): 8.132006740422481
Clark-Evans index (10 runs): 0.9050970765059173

Can we trust this metric now? A little more, but it would be great to have a p-value. We can calculate the p-value by asking ‘In what fraction of the runs was the observed NN distance smaller than random?’.

runs_rnd_smaller = np.array(list_rnd_dist) <= spot_dist_avg
print(runs_rnd_smaller)

[False False False False False False False False False False]

The fraction of Trues in this list is our p-value.

Clustered spots → NND_observed is small → most simulations produce a larger NND → few simulations have NND_sim ≤ NND_obs → p close to 0
Random spots → NND_observed is typical → roughly half the simulations fall below it → p ≈ 0.5

pval = runs_rnd_smaller.mean()
print("p-value for CE index:", pval)

p-value for CE index: 0.0

8. Repeat Measurements for All Cells#

So far, we have focussed on one single cell. Of course, we need to analyse all cells.

✍️ Exercise: Summarise the per-cell clustering analysis as a function#

To make processing individual cells simpler, we can summarise our above analysis in a function called nn_dist. The function should have the following features:

Parameters:

spots : (n, 2) numpy array Spot coordinates in (row, col) / (y, x) order.
mask : 2D numpy array Label mask where each cell is identified by a unique integer ID.
cell_id : int ID of the cell in mask that spots belong to.
n_repeats : int Number of CSR simulations to run.

Returns

spot_dist_avg : float Mean nearest-neighbor distance of the observed spots.
rnd_dist_avg : float Mean of the per-simulation mean nearest-neighbor distances under CSR.
clark_evans : float Clark-Evans index: spot_dist_avg / rnd_dist_avg.
p_value : float Fraction of CSR simulations with mean NND <= spot_dist_avg.

# or use:
# from bobiac_tools import nn_dist


def nn_dist(
    spots: np.ndarray, mask: np.ndarray, cell_id: int, n_repeats: int
) -> tuple[float, float, float, float]:
    """
    Test whether spots in a cell are randomly distributed using nearest-neighbor
    distances and a Monte Carlo simulation of Complete Spatial Randomness (CSR).

    For each simulation run, n_repeats random point patterns are generated inside
    the cell mask and their mean nearest-neighbor distance is computed. The
    Clark-Evans index (observed / random NND) and a Monte Carlo p-value are
    derived from this null distribution.

    Parameters
    ----------
    spots : (n, 2) numpy array
        Spot coordinates in (row, col) / (y, x) order.
    mask : 2D numpy array
        Label mask where each cell is identified by a unique integer ID.
    cell_id : int
        ID of the cell in mask that spots belong to.
    n_repeats : int
        Number of CSR simulations to run.

    Returns
    -------
    spot_dist_avg : float
        Mean nearest-neighbor distance of the observed spots.
    rnd_dist_avg : float
        Mean of the per-simulation mean nearest-neighbor distances under CSR.
    clark_evans : float
        Clark-Evans index: spot_dist_avg / rnd_dist_avg.
        < 1 indicates clustering, > 1 indicates dispersion.
    p_value : float
        Fraction of CSR simulations with mean NND <= spot_dist_avg.
        Small values (< 0.05) indicate significant clustering.
    """
    n_points = spots.shape[0]

    dist_spots, _ = KDTree(spots).query(spots, k=2)
    spot_dist_avg = dist_spots[
        :, 1
    ].mean()  # k=2: column 0 is self, column 1 is nearest neighbor

    list_rnd_dist = []
    for _i in range(n_repeats):
        rnd_points = random_points_in_mask(mask, cell_label=cell_id, n=n_points)
        dist_rnd, _ = KDTree(rnd_points).query(rnd_points, k=2)
        rnd_dist_avg = dist_rnd[:, 1].mean()
        list_rnd_dist.append(rnd_dist_avg)

    rnd_dist_avg = np.array(list_rnd_dist).mean()
    clark_evans = spot_dist_avg / rnd_dist_avg
    pval = (np.array(list_rnd_dist) <= spot_dist_avg).mean()

    return (
        spot_dist_avg,
        rnd_dist_avg,
        clark_evans,
        pval,
    )


nn_dist(spots, mask, cell_id, 100)

(np.float64(7.360255526882802),
 np.float64(7.468141977074545),
 np.float64(0.985553776224002),
 np.float64(0.39))

channel = 2
n_repeats = 10

results = []
for _, row in df_cell.iterrows():
    # get the spots for channel 0 and the current cell
    spots_in_cell = df_spots[
        (df_spots["label_cell"] == row["label"]) & (df_spots["channel"] == channel)
    ][["y", "x"]].values
    # calculate the NN distance, Clark-Evans index and p-value for that cell
    results.append(nn_dist(spots_in_cell, mask, row["label"], n_repeats=n_repeats))

df_cell[["nnd_observed", "nnd_random", "clark_evans", "p_value"]] = results
df_cell

	area	label	count	nnd_observed	nnd_random	clark_evans	p_value
0	5634.0	1	37	6.756012	6.924861	0.975617	0.4
1	4728.0	2	34	6.870816	7.063149	0.972770	0.5
2	7677.0	3	36	7.618762	8.190042	0.930247	0.0
3	6824.0	4	44	6.343310	6.879484	0.922062	0.2
4	9638.0	5	56	6.074178	7.106470	0.854739	0.0
...	...	...	...	...	...	...	...
71	6577.0	80	41	6.366875	7.422158	0.857820	0.0
72	8947.0	81	37	6.781684	8.195194	0.827520	0.1
73	2622.0	82	18	7.013325	6.818760	1.028534	0.4
74	1669.0	84	12	7.230385	7.857231	0.920221	0.3
75	4207.0	86	35	6.692059	6.049204	1.106271	0.8

76 rows × 7 columns

df_cell.loc[df_cell["label"] == cell_id]

	area	label	count	nnd_observed	nnd_random	clark_evans	p_value
9	17676.0	11	74	6.978486	8.249791	0.845898	0.1

sns.histplot(
    df_cell,
    x="clark_evans",
)
plt.axvline(df_cell["clark_evans"].mean(), c="red", lw=10)

<matplotlib.lines.Line2D at 0x115379810>

../../../_images/6299a728c039d3baa57ddcce554bd7d6837926dbbc1525502be1476cc27928d5.png

overlay_labels(
    image=image[2, :, :],
    label_mask=mask,
    df=df_cell,
    id_col="label",
    measurement_col="clark_evans",
)

../../../_images/f3391a7c3d743d9e611d0c6a9c4659deab7788217651bfdfec5183fe59458845.png

(<Figure size 800x800 with 2 Axes>, <Axes: title={'center': 'clark_evans'}>)

9. Perform the Same Calculation on the Other Channels#

Let’s use the tools we developed above to see if the spots in the other channels.

Protein B

channel = 1
n_repeats = 10

results = []
for _, row in df_cell.iterrows():
    # get the spots for channel 0 and the current cell
    spots_in_cell = df_spots[
        (df_spots["label_cell"] == row["label"]) & (df_spots["channel"] == channel)
    ][["y", "x"]].values
    # calculate the NN distance, Clark-Evans index and p-value for that cell
    results.append(nn_dist(spots_in_cell, mask, row["label"], n_repeats=n_repeats))

df_cell[["nnd_observed", "nnd_random", "clark_evans", "p_value"]] = results
df_cell

	area	label	count	nnd_observed	nnd_random	clark_evans	p_value
0	5634.0	1	37	4.167969	4.564722	0.913083	0.1
1	4728.0	2	34	3.905908	4.780596	0.817034	0.0
2	7677.0	3	36	4.235827	5.130219	0.825662	0.0
3	6824.0	4	44	4.595475	4.768810	0.963652	0.3
4	9638.0	5	56	4.891864	5.448673	0.897808	0.1
...	...	...	...	...	...	...	...
71	6577.0	80	41	4.507801	4.889029	0.922024	0.0
72	8947.0	81	37	4.271835	5.095960	0.838279	0.0
73	2622.0	82	18	3.185283	3.613965	0.881382	0.1
74	1669.0	84	12	2.656640	3.111253	0.853881	0.0
75	4207.0	86	35	3.621850	4.064461	0.891102	0.2

76 rows × 7 columns

sns.histplot(
    df_cell,
    x="clark_evans",
)
plt.axvline(df_cell["clark_evans"].mean(), c="red", lw=10)

<matplotlib.lines.Line2D at 0x1154aa850>

../../../_images/777d2c060f355891b26d41cc10ea102e33339b8064e020cc722207ccefd6aaab.png

overlay_labels(
    image=image[3, :, :],
    label_mask=mask,
    df=df_cell,
    id_col="label",
    measurement_col="clark_evans",
)

../../../_images/a22c87c051ace1be90e909de2927027cbcffda2bfad3dfbae6a6541d7011b31a.png

(<Figure size 800x800 with 2 Axes>, <Axes: title={'center': 'clark_evans'}>)

Protein C

channel = 2
n_repeats = 10

results = []
for _, row in df_cell.iterrows():
    # get the spots for channel 0 and the current cell
    spots_in_cell = df_spots[
        (df_spots["label_cell"] == row["label"]) & (df_spots["channel"] == channel)
    ][["y", "x"]].values
    # calculate the NN distance, Clark-Evans index and p-value for that cell
    results.append(nn_dist(spots_in_cell, mask, row["label"], n_repeats=n_repeats))

df_cell[["nnd_observed", "nnd_random", "clark_evans", "p_value"]] = results
df_cell

	area	label	count	nnd_observed	nnd_random	clark_evans	p_value
0	5634.0	1	37	6.756012	6.553187	1.030951	0.8
1	4728.0	2	34	6.870816	7.297388	0.941545	0.4
2	7677.0	3	36	7.618762	8.277627	0.920404	0.2
3	6824.0	4	44	6.343310	6.747335	0.940121	0.2
4	9638.0	5	56	6.074178	7.017423	0.865585	0.0
...	...	...	...	...	...	...	...
71	6577.0	80	41	6.366875	7.179537	0.886809	0.2
72	8947.0	81	37	6.781684	8.287768	0.818276	0.1
73	2622.0	82	18	7.013325	7.077613	0.990917	0.5
74	1669.0	84	12	7.230385	6.531853	1.106943	0.6
75	4207.0	86	35	6.692059	6.389435	1.047363	0.6

76 rows × 7 columns

sns.histplot(
    df_cell,
    x="clark_evans",
)
plt.axvline(df_cell["clark_evans"].mean(), c="red", lw=10)

<matplotlib.lines.Line2D at 0x11580b390>

../../../_images/afaaa092976de5d74f2d6728584f664087b68135a246458d23255261201831f1.png

overlay_labels(
    image=image[4, :, :],
    label_mask=mask,
    df=df_cell,
    id_col="label",
    measurement_col="clark_evans",
)

../../../_images/afa38bd40637557d3b6761e69119ffa6351055c5a63d4aae4fac8cfb6d9c0259.png

(<Figure size 800x800 with 2 Axes>, <Axes: title={'center': 'clark_evans'}>)

It seems like protein A is not clustered, protein B is clustering and protein C is a little ambiguous.

10. Measure the Length Scale of Clustering: Ripleys L Function#

We know that Protein B is clustered compared to random. However, we don’t know what scale the clustering is occurring on.

A common method to investigate the length scales of spatial patterns is Ripley’s L function. The idea is to calculate the density of spots around each spot for different radii and compare that to a random point pattern. The density of a random pattern is independent of the radius, while clustered points will exhibit higher density at shorter radii than longer ones. This way we can infer the average cluster size and distance between clusters.

Here, we only demonstrate the concept on the spot pattern of protein B in one cell.

spots = df_spots.loc[(df_spots["channel"] == 1) & (df_spots["label_cell"] == cell_id)][
    ["y", "x"]
].to_numpy()
overlay_labels(label_mask=mask, coordinates=spots, focus_object=cell_id)

../../../_images/490eb7bad9e67cbceb922a9d851348e435f78dc0660c372dc26e8bf5d6b8b765.png

(<Figure size 800x586.081 with 1 Axes>,
 <Axes: title={'center': 'Segmentation overlay'}>)

r_values, L_obs, L_sims = ripleys_l(spots, mask, cell_id=cell_id, n_repeats=200)

plt.plot(r_values, L_obs, color="black", label="observed")
plt.fill_between(
    r_values,
    np.percentile(L_sims, 2.5, axis=0),
    np.percentile(L_sims, 97.5, axis=0),
    alpha=0.3,
    label="95% CSR envelope",
)

plt.axhline(0, color="grey", linestyle="--")
plt.xlabel("r (pixels)")
plt.ylabel("L(r) - r")
plt.legend()

<matplotlib.legend.Legend at 0x1159e6e90>

../../../_images/d7a25598292c2c67e046930514bbf88ac0dcac152f879d7440afbb39941ee689.png

We observe clustering with scale 2-10 pixels, where the observed curve is above the 95% CSR envelope. Moreover, we see a corresponding dispersion from 15 pixels onward. This is because our random null distribution has the same number of points as the observed distribution.

Conclusion#

Here we have learned how to assess whether a spot pattern is clustered compared to random. Note the following:

The metric we used to measure clustering (nearest neighbor distance) was a choice. Other metrics, such as number of spots in a certain radius, distance to k nearest neighbors, and many others, are valid. Can you think of a spot pattern where the nearest neighbor distance is not an appropriate metric?
We chose a completely random spatial pattern as a null distribution. This was a choice. An equally valid option would have been to translate, rotate, etc., the spot pattern.
The formalism with which we defined spatial relationships (KDTree / Euclidean distance nearest neighbor search) was a choice. An equally valid approach would have been to divide space using a Voronoi diagram, or to use the distance between cell membranes.

All these choices need to be made with your scientific question in mind.

Notebook 1: Does a protein cluster?

Contents

Notebook 1: Does a protein cluster?#

Overview#

Dataset#

0. Importing libraries#

2. Load image, mask and spot table#

Visually inspect the images, masks and spots#

3. Make An Analysis Plan#

Question#

Approach#

Metric for Clustering#

Workflow#

1. Generate a reference for comparison#

2. Measure the nearest neighbor distance#

3. Calculate the average nearest neighbour distance#

4. Calculate the Clark-Evans Index#

5. Calculate the p-value of the Clark-Evans Index#

4. Assemble the Protein Dataset and the Random Dataset#

Protein Dataset: Map spots to cells and count spots per cell#

Random Dataset: Generate random points#

5. Calculate Nearest Neighbor Distances#

Identifying Neighbors: KD Trees#

6. Calculate the Clark-Evans Index for the protein and random datasets#

7. Repeat Steps 4-6 with Multiple Random Patterns#

Concept#

8. Repeat Measurements for All Cells#

✍️ Exercise: Summarise the per-cell clustering analysis as a function#

9. Perform the Same Calculation on the Other Channels#

10. Measure the Length Scale of Clustering: Ripleys L Function#

Conclusion#