Skip to content

Spatial stats #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
AnthonyChristidis wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Spatial stats #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
87a4fe2
Added complete spatial statistics pipeline and documentation notebooks
AnthonyChristidis May 10, 2026
51c6572
Resolved merge conflict in Results.py
AnthonyChristidis May 10, 2026
f780412
Fixed expression bar plot formatting and updated neighborhood gating …
AnthonyChristidis May 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions environment.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,10 @@ dependencies:
- lap
- poppy
- zernike
- squidpy
- anndata
- esda
- libpysal
- -f "https://download.pytorch.org/whl/torch_stable.html"
- torch==1.8.1+cu111
- torchvision==0.9.1+cu111
Expand Down
326 changes: 326 additions & 0 deletions notebooks/SpatialStats_SingleTimepoint.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,326 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# oyLabImaging Spatial Statistics Pipeline\n",
"## Single-Timepoint Example\n",
"\n",
"This notebook demonstrates how to extract, calculate, and visualize spatiotemporal interactions from segmented microscopy data. It covers:\n",
"1. Global Spatial Autocorrelation (Moran's I, Geary's C)\n",
"2. Local Spatial Hotspot Detection (Local Moran's I)\n",
"3. Spatial Expression Mapping\n",
"4. Categorical Neighborhood Enrichment\n",
"5. Interactive Napari Visualization (Stacked Metric Layers)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2\n",
"%gui qt\n",
"%matplotlib inline\n",
"\n",
"import sys\n",
"import os\n",
"import dill\n",
"import numpy as np\n",
"import pandas as pd\n",
"import anndata as ad\n",
"\n",
"# Point to our modified package directory\n",
"sys.path.insert(0, \"../oyLabImaging\")\n",
"from oyLabImaging import Metadata\n",
"from oyLabImaging.Processing.Results import results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1. Data Loading & Environment Setup"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Define the path to the dataset\n",
"data_path = '/bigstore/pirlo/Images2025/Sanne/2025.11.23_VME221_IFNdecoy_B18Rtransient/pSTAT1_cntrls_1/'\n",
"print(\"Loading results.pickle...\")\n",
"\n",
"# Load the master results object\n",
"with open(os.path.join(data_path, \"results.pickle\"), \"rb\") as f:\n",
" R = dill.load(f)\n",
"R.pth = data_path\n",
"\n",
"# Load the individual Position files (PosLbls) to get the single-cell data\n",
"for pos_name in R.PosNames:\n",
" pkl_file = os.path.join(data_path, \"PosLbls\", f\"{pos_name}.pkl\")\n",
" if os.path.exists(pkl_file):\n",
" with open(pkl_file, \"rb\") as f:\n",
" P = dill.load(f)\n",
" P.pth = data_path\n",
" R.PosLbls[pos_name] = P\n",
"\n",
"# Grab the positions we want to compare\n",
"test_positions = ['B6-Site_0', 'B5-Site_0']\n",
"\n",
"print(f\"✓ Loaded successfully\")\n",
"print(f\" Positions we are testing: {test_positions}\")\n",
"print(f\" Timepoints found: {len(R.frames)}\")\n",
"print(f\" Exact Channels found: {list(R.channels)}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. Discrete Cell Classification (Quadrant Gating)\n",
"Neighborhood Enrichment is designed for discrete cell types. We mimic flow-cytometry \"Quadrant Gating\" by finding the top 10% of expressors for each marker and categorizing every cell."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for pos in test_positions:\n",
" # 1. Find the global 90th percentile for the markers \n",
" all_red = np.concatenate([R.PosLbls[pos].mean('Red')[t] for t in range(len(R.frames)) if R.PosLbls[pos].num[t] > 0])\n",
" all_farred = np.concatenate([R.PosLbls[pos].mean('FarRed')[t] for t in range(len(R.frames)) if R.PosLbls[pos].num[t] > 0])\n",
" \n",
" red_90th = np.percentile(all_red, 90)\n",
" farred_90th = np.percentile(all_farred, 90)\n",
"\n",
" # 2. Classify each cell based on these biological thresholds\n",
" for t in range(len(R.frames)):\n",
" if R.PosLbls[pos].num[t] > 0:\n",
" r_vals = R.PosLbls[pos].mean('Red')[t]\n",
" fr_vals = R.PosLbls[pos].mean('FarRed')[t]\n",
" \n",
" states = []\n",
" for r, fr in zip(r_vals, fr_vals):\n",
" if r > red_90th and fr > farred_90th:\n",
" states.append('Double+') # High Red, High FarRed\n",
" elif r > red_90th:\n",
" states.append('Red+') # High Red only\n",
" elif fr > farred_90th:\n",
" states.append('FarRed+') # High FarRed only\n",
" else:\n",
" states.append('Low') # Background/Bystander cells\n",
" \n",
" R.PosLbls[pos].framelabels[t].regionprops['Cell_State'] = states"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. Run Spatial Statistics\n",
"Calculating Global and Local autocorrelation, Bivariate interactions, and Neighborhood Enrichment."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"target_channels = ['Red', 'FarRed']\n",
"biv_pairs = [('Red', 'FarRed')]\n",
"\n",
"print(\"\\n--- Running Spatial Stats Integration ---\")\n",
"R.calculate_spatial_stats(\n",
" Position=test_positions, \n",
" metrics=[\n",
" 'morans_i', \n",
" 'gearys_c', \n",
" 'neighborhood_enrichment', \n",
" 'bivariate_moran', \n",
" 'local_morans_i', \n",
" 'local_bivariate_moran'\n",
" ],\n",
" channels=target_channels,\n",
" bivariate_pairs=biv_pairs,\n",
" cluster_key='Cell_State', # Uses the discrete Cell_State column we generated\n",
" n_neighs=10, # Uses 10 neighbors for this specific dataset \n",
" export_h5ad=True, # Exports AnnData for external use\n",
" save=False \n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4. Diagnostics: View Raw Outputs & Counts"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"\\n--- Diagnostic: Cell Counts per Frame ---\")\n",
"counts_dict = {'Frame': R.frames}\n",
"for pos in test_positions:\n",
" counts_dict[pos] = R.PosLbls[pos].num\n",
"df_counts = pd.DataFrame(counts_dict)\n",
"print(df_counts) \n",
"\n",
"print(\"\\n--- Viewing Raw Local Statistics (First 5 Cells) ---\")\n",
"df = R.PosLbls[test_positions[0]].framelabels[0].regionprops\n",
"print(df[['mean_Red', 'local_moran_I_Red', 'local_moran_q_Red', 'local_moran_p_Red']].head())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5. Visualization: Generating Matplotlib Plots\n",
"Showcasing all functionality: Global metrics, grouped bar chart summaries, static X/Y spatial maps, and categorical enrichment."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Define exact colors to match the biology using the exact channel names\n",
"my_colors = {\n",
" 'Red': '#e41a1c', # Univariate: Red\n",
" 'FarRed': '#800080', # Univariate: Purple\n",
" 'Red vs FarRed': '#ff7f00', # Bivariate: Orange\n",
" ('Red+', 'Red+'): '#e41a1c', # Nhood: Red - Red\n",
" ('Red+', 'FarRed+'): '#ff7f00', # Nhood: Red - FarRed\n",
" ('Double+', 'Double+'): '#4daf4a', # Nhood: Double+ - Double+ (Green)\n",
" ('Double+', 'Low'): '#377eb8', # Nhood: Double+ - Low (Blue)\n",
"}\n",
"\n",
"# --- A. GLOBAL METRICS ---\n",
"print(\"\\n--- Plotting Univariate: Moran's I ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='morans_i', channels=target_channels, custom_colors=my_colors)\n",
"\n",
"print(\"\\n--- Plotting Univariate: Geary's C ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='gearys_c', channels=target_channels, custom_colors=my_colors)\n",
"\n",
"print(\"\\n--- Plotting Bivariate Moran's I ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='bivariate_moran', custom_colors=my_colors)\n",
"\n",
"\n",
"# --- B. SUMMARY BAR CHARTS ---\n",
"# Because there is only 1 timepoint, 'summary' forces a clean grouped bar chart.\n",
"print(\"\\n--- Plotting Expression Summary (Mean Intensity) ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='expression', channels=target_channels, plot_type='summary', custom_colors=my_colors)\n",
"\n",
"print(\"\\n--- Plotting Local Moran's I Summary (% Hotspots) ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='local_morans_i', channels=target_channels, plot_type='summary', custom_colors=my_colors)\n",
"\n",
"print(\"\\n--- Plotting Local Bivariate Moran's I Summary (% Hotspots) ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='local_bivariate_moran', plot_type='summary', custom_colors=my_colors)\n",
"\n",
"\n",
"# --- C. SPATIAL SNAPSHOT MAPS ---\n",
"print(\"\\n--- Plotting Expression (Spatial Map Snapshot) ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='expression', channels=target_channels, plot_type='spatial_map', frame_idx=0)\n",
"\n",
"print(\"\\n--- Plotting Local Univariate (Spatial Map Snapshot) ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='local_morans_i', channels=target_channels, plot_type='spatial_map', frame_idx=0)\n",
"\n",
"print(\"\\n--- Plotting Local Bivariate (Spatial Map Snapshot) ---\")\n",
"R.plot_spatial_stats(Position=test_positions, metric='local_bivariate_moran', plot_type='spatial_map', frame_idx=0)\n",
"\n",
"\n",
"# --- D. NEIGHBORHOOD ENRICHMENT ---\n",
"print(\"\\n--- Plotting Neighborhood Enrichment Bar Chart ---\")\n",
"my_pairs = [\n",
" ('Red+', 'Red+'), \n",
" ('Red+', 'FarRed+'), \n",
" ('Double+', 'Double+'), \n",
" ('Double+', 'Low')\n",
"]\n",
"R.plot_spatial_stats(Position=test_positions, metric='neighborhood_enrichment', nhood_pairs=my_pairs, custom_colors=my_colors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 6. Interactive Napari Visualization: Stacked Metric Layers\n",
"This loads the underlying TIFF image and stacks 3 toggleable point layers (Expression, Local Moran's I, and Bivariate Moran's I) in the same viewer!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"\\n--- Opening Napari to view EXPRESSION & HOTSPOTS ---\")\n",
"\n",
"# Call 1: Clears the viewer, and loads the Red expression map\n",
"viewer = R.show_spatial_map_napari(\n",
" pos=test_positions[0], \n",
" Channel='Red',\n",
" metric='expression', \n",
" frame_idx=0, \n",
" size=10,\n",
" load_images=True, # Loads the TIF because it's only 1 frame!\n",
" clear_viewer=True \n",
")\n",
"\n",
"# Call 2: Stacks the Red Hotspots on top of the same viewer\n",
"R.show_spatial_map_napari(\n",
" pos=test_positions[0], \n",
" Channel='Red',\n",
" metric='local_morans_i', \n",
" frame_idx=0, \n",
" size=10,\n",
" load_images=False \n",
")\n",
"\n",
"# Call 3: Stacks the Red vs FarRed Bivariate Hotspots on top\n",
"R.show_spatial_map_napari(\n",
" pos=test_positions[0], \n",
" Channel='Red vs FarRed',\n",
" metric='local_bivariate_moran', \n",
" frame_idx=0, \n",
" size=10,\n",
" load_images=False \n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Loading