Benchmarking stGP using Excitatory Neuron in Human Aging Brains
This tutorial walks through the full stGP pipeline on the human aging brain MERFISH dataset (Jeffries et al., Nature 2025), focusing on Excitatory Neurons (ext) as the target cell type. Twelve tissue sections span donors aged 15–87 years. The laminar strcuture of excitatory neurons provide an opportunity for evaluating the spatial structure after decomposing the spatial and temporal effects.
1. Setup
[1]:
%matplotlib inline
import os
import pickle
import sys
import warnings
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
from IPython.display import display
PROJECT_DIR = Path.cwd()
os.chdir(PROJECT_DIR)
sys.path.insert(0, str(PROJECT_DIR))
sys.path.insert(0, str(PROJECT_DIR.parent))
from plots import (
plot_alpha_trajectory,
plot_representative_continuous_tiles,
plot_w_heatmap_vertical,
representative_slices_by_age,
set_nature_style,
)
set_nature_style()
warnings.filterwarnings("ignore", category=FutureWarning)
RAW_DATA_DIR = Path("/import/home2/share/byual/HumanBrainMERFISH+sc_Nature2025_Jeffries")
DATA_QC = Path("data/qc/human_merfish_qc.h5ad")
DATA_PROC = Path("data/processed")
RESULTS_DIR = Path("Results/stgp")
FIGURES_DIR = Path("Figure/ext")
FIGURES_DIR.mkdir(parents=True, exist_ok=True)
CELLTYPE = "ext" # target cell type
adata_ext = sc.read_h5ad(DATA_PROC / f"{CELLTYPE}.h5ad")
print(f"Working directory: {PROJECT_DIR}")
Working directory: /import/home4/byual/stGP-0707/RealData_HumanBrainMERFISH
2. Fitting stGP
stGP models gene expression in each tissue slice as a sum of p latent programs. Each program is characterised by:
W (gene loadings,
p × G): non-negative gene weights defining the programH (cell scores,
N × p): overall activity of each program in each cellb (spatial field,
N × p): spatially smooth residual componentα (age effect,
p × S): how program amplitude varies across slices/ages:math:`sigma_{age}^2`, :math:`tau_{spa}^2` (variance components,
p × 2): strength of temporal and spatial effects of each program
The model rank p is selected automatically by greedy forward selection.
[2]:
import time
from stgp.estimation import fit_pfactor_auto
from stgp.kernels import (
bandwidth_select_spatial, bandwidth_select_temporal,
build_K_age, build_K_spa_list_from_stacked
)
from stgp.preprocessing import standardize_coords_list, log1p_norm_centered_list
OUT_DIR = RESULTS_DIR / CELLTYPE
OUT_DIR.mkdir(parents=True, exist_ok=True)
PKL_PATH = OUT_DIR / "stgp_result.pkl"
age_arr = pd.to_numeric(adata_ext.obs["age"], errors="coerce").to_numpy(float)
groups = adata_ext.obs["id_region"].astype(str).to_numpy()
uniq, inv = np.unique(groups, return_inverse=True)
idx_per_group = [np.sort(np.where(inv == t)[0]) for t in range(len(uniq))]
adata_temp = adata_ext.copy()
sc.pp.normalize_total(adata_temp, target_sum=1e3)
adata_prep = sc.pp.log1p(adata_temp, copy = True)
Y_list = [adata_prep.X[ix].toarray() for ix in idx_per_group]
Y_list, _ = log1p_norm_centered_list(Y_list, target_sum = 1000)
nlist = np.array([len(ix) for ix in idx_per_group])
ages = np.array([age_arr[ix[0]] for ix in idx_per_group])
sort_ord = np.argsort(ages); ages = ages[sort_ord]
slices = uniq.copy(); slices = slices[sort_ord]
nlist = nlist[sort_ord]
Y_list = [Y_list[i] for i in sort_ord]
[3]:
# ── Build GP kernels ─────────────────────────────────────────────────────
coords_list = standardize_coords_list([adata_ext.obsm["spatial"][ix] for ix in idx_per_group])
coords_list = [coords_list[i] for i in sort_ord]
gamma_spa = bandwidth_select_spatial(coords_list, frac=0.01, rho=0.6)
gamma_age = bandwidth_select_temporal(ages, rho=np.exp(-1.5))
print(f" gamma_spa = {gamma_spa:.4f} | gamma_age = {gamma_age:.4f}")
K_age = build_K_age(ages, gamma_age, kernel="rbf", standardize=True)
K_spa_list = build_K_spa_list_from_stacked(
np.vstack(coords_list), nlist, gamma_spa, standardize=False, jitter=1e-6
)
gamma_spa = 0.1473 | gamma_age = 1.0997
[4]:
t0 = time.perf_counter()
res = fit_pfactor_auto(
Y_list=Y_list, Nlist=nlist, K_age=K_age, Kspa_list=K_spa_list,
p_max=10, k=15,
inner_rank1_tol=1e-4, rel_improve_total_tol=0.002, backfit_tol=1e-4, prune_energy_frac = 0.005,
random_state=0, verbose=1,
)
print(f"Runtime: {time.perf_counter() - t0:.1f}s | programs selected: {res['W'].shape[0]}")
res["gamma_age"] = gamma_age; res["gamma_spa"] = gamma_spa
with open(PKL_PATH, "wb") as f:
pickle.dump(res, f)
print(f"Saved: {PKL_PATH}")
[sweep=001] dW_rel=1.218e-01 dTheta_rel=1.983e-02 time=1.873e+01
[sweep=002] dW_rel=6.676e-02 dTheta_rel=8.814e-03 time=1.305e+01
[sweep=003] dW_rel=2.279e-02 dTheta_rel=3.723e-03 time=4.392e+00
[sweep=004] dW_rel=1.621e-02 dTheta_rel=2.770e-03 time=5.251e+00
[sweep=005] dW_rel=1.241e-02 dTheta_rel=2.481e-03 time=3.715e+00
[sweep=006] dW_rel=9.979e-03 dTheta_rel=2.230e-03 time=3.517e+00
[sweep=007] dW_rel=6.024e-02 dTheta_rel=1.084e-02 time=4.760e+00
[sweep=008] dW_rel=2.301e-02 dTheta_rel=1.035e-02 time=4.719e+00
[sweep=009] dW_rel=7.135e-02 dTheta_rel=1.357e-02 time=5.090e+00
[sweep=010] dW_rel=3.012e-02 dTheta_rel=1.148e-02 time=5.086e+00
[sweep=011] dW_rel=5.342e-02 dTheta_rel=8.292e-03 time=6.332e+00
[sweep=012] dW_rel=2.314e-02 dTheta_rel=8.913e-03 time=4.721e+00
[sweep=013] dW_rel=2.005e-02 dTheta_rel=7.392e-03 time=4.580e+00
[sweep=014] dW_rel=7.595e-02 dTheta_rel=6.371e-03 time=5.787e+00
[sweep=015] dW_rel=1.148e-01 dTheta_rel=1.352e-02 time=6.291e+00
[sweep=016] dW_rel=3.844e-01 dTheta_rel=4.893e-02 time=1.202e+01
[sweep=017] dW_rel=2.464e-02 dTheta_rel=1.177e-02 time=6.969e+00
[sweep=018] dW_rel=2.072e-02 dTheta_rel=1.006e-02 time=5.821e+00
[sweep=019] dW_rel=5.216e-02 dTheta_rel=9.616e-03 time=5.774e+00
[sweep=020] dW_rel=1.504e-01 dTheta_rel=1.497e-02 time=7.798e+00
[sweep=021] dW_rel=1.504e-02 dTheta_rel=9.105e-03 time=5.902e+00
[sweep=022] dW_rel=1.365e-02 dTheta_rel=7.652e-03 time=5.393e+00
[sweep=023] dW_rel=1.226e-02 dTheta_rel=6.611e-03 time=4.988e+00
[sweep=024] dW_rel=1.097e-02 dTheta_rel=5.676e-03 time=4.703e+00
[sweep=025] dW_rel=9.814e-03 dTheta_rel=4.970e-03 time=4.439e+00
[sweep=026] dW_rel=8.805e-03 dTheta_rel=4.571e-03 time=4.388e+00
[sweep=027] dW_rel=7.923e-03 dTheta_rel=4.068e-03 time=4.352e+00
[sweep=028] dW_rel=7.158e-03 dTheta_rel=3.713e-03 time=4.251e+00
[sweep=029] dW_rel=3.839e-02 dTheta_rel=3.312e-03 time=4.238e+00
[sweep=030] dW_rel=1.419e-02 dTheta_rel=3.542e-03 time=4.625e+00
[sweep=031] dW_rel=4.370e-03 dTheta_rel=3.756e-03 time=4.265e+00
[sweep=032] dW_rel=4.049e-03 dTheta_rel=3.373e-03 time=4.246e+00
[sweep=033] dW_rel=3.775e-03 dTheta_rel=3.024e-03 time=4.066e+00
[sweep=034] dW_rel=3.536e-03 dTheta_rel=2.665e-03 time=3.810e+00
[sweep=035] dW_rel=3.323e-03 dTheta_rel=2.399e-03 time=3.794e+00
[sweep=036] dW_rel=3.132e-03 dTheta_rel=2.274e-03 time=3.671e+00
[sweep=037] dW_rel=2.958e-03 dTheta_rel=2.097e-03 time=3.414e+00
[sweep=038] dW_rel=9.045e-02 dTheta_rel=2.311e-03 time=4.553e+00
[sweep=039] dW_rel=3.684e-03 dTheta_rel=1.962e-03 time=4.290e+00
[sweep=040] dW_rel=2.609e-03 dTheta_rel=1.462e-03 time=4.069e+00
[sweep=041] dW_rel=2.460e-03 dTheta_rel=1.394e-03 time=3.542e+00
[sweep=042] dW_rel=2.334e-03 dTheta_rel=1.318e-03 time=3.098e+00
[sweep=043] dW_rel=2.225e-03 dTheta_rel=1.370e-03 time=2.896e+00
[sweep=044] dW_rel=2.128e-03 dTheta_rel=1.318e-03 time=2.743e+00
[sweep=045] dW_rel=2.038e-03 dTheta_rel=1.308e-03 time=3.074e+00
[sweep=046] dW_rel=1.954e-03 dTheta_rel=1.284e-03 time=3.114e+00
[sweep=047] dW_rel=1.880e-03 dTheta_rel=1.235e-03 time=3.141e+00
[sweep=048] dW_rel=1.808e-03 dTheta_rel=1.163e-03 time=2.726e+00
[sweep=049] dW_rel=1.742e-03 dTheta_rel=1.239e-03 time=2.961e+00
[sweep=050] dW_rel=1.679e-03 dTheta_rel=1.041e-03 time=2.746e+00
[sweep=051] dW_rel=1.621e-03 dTheta_rel=1.182e-03 time=2.761e+00
[sweep=052] dW_rel=1.567e-03 dTheta_rel=1.073e-03 time=2.903e+00
[sweep=053] dW_rel=1.516e-03 dTheta_rel=1.010e-03 time=2.739e+00
[sweep=054] dW_rel=1.467e-03 dTheta_rel=9.738e-04 time=3.146e+00
[sweep=055] dW_rel=1.421e-03 dTheta_rel=1.058e-03 time=3.094e+00
[sweep=056] dW_rel=1.379e-03 dTheta_rel=9.958e-04 time=2.765e+00
[sweep=057] dW_rel=1.336e-03 dTheta_rel=1.012e-03 time=2.618e+00
[sweep=058] dW_rel=1.298e-03 dTheta_rel=8.135e-04 time=2.627e+00
[sweep=059] dW_rel=1.260e-03 dTheta_rel=8.382e-04 time=2.550e+00
[sweep=060] dW_rel=1.225e-03 dTheta_rel=9.158e-04 time=2.945e+00
[sweep=061] dW_rel=1.192e-03 dTheta_rel=8.612e-04 time=2.705e+00
[sweep=062] dW_rel=1.158e-03 dTheta_rel=7.832e-04 time=2.593e+00
[sweep=063] dW_rel=1.129e-03 dTheta_rel=8.074e-04 time=2.677e+00
[sweep=064] dW_rel=1.098e-03 dTheta_rel=7.496e-04 time=2.427e+00
[sweep=065] dW_rel=1.071e-03 dTheta_rel=7.867e-04 time=2.722e+00
[sweep=066] dW_rel=1.046e-03 dTheta_rel=7.951e-04 time=2.709e+00
[sweep=067] dW_rel=1.019e-03 dTheta_rel=6.357e-04 time=2.337e+00
[sweep=068] dW_rel=9.945e-04 dTheta_rel=7.402e-04 time=2.558e+00
[sweep=069] dW_rel=9.707e-04 dTheta_rel=7.351e-04 time=2.568e+00
[sweep=070] dW_rel=9.484e-04 dTheta_rel=6.862e-04 time=2.655e+00
[sweep=071] dW_rel=9.272e-04 dTheta_rel=7.563e-04 time=2.476e+00
[sweep=072] dW_rel=9.050e-04 dTheta_rel=6.131e-04 time=2.447e+00
[sweep=073] dW_rel=8.856e-04 dTheta_rel=6.641e-04 time=2.493e+00
[sweep=074] dW_rel=8.650e-04 dTheta_rel=6.834e-04 time=2.512e+00
[sweep=075] dW_rel=8.485e-04 dTheta_rel=6.747e-04 time=2.706e+00
[sweep=076] dW_rel=8.292e-04 dTheta_rel=4.801e-04 time=2.303e+00
[sweep=077] dW_rel=8.134e-04 dTheta_rel=6.730e-04 time=2.588e+00
[sweep=078] dW_rel=7.942e-04 dTheta_rel=5.725e-04 time=2.380e+00
[sweep=079] dW_rel=7.791e-04 dTheta_rel=7.191e-04 time=2.518e+00
[sweep=080] dW_rel=7.621e-04 dTheta_rel=4.509e-04 time=2.325e+00
[sweep=081] dW_rel=7.481e-04 dTheta_rel=5.570e-04 time=2.377e+00
[sweep=082] dW_rel=7.344e-04 dTheta_rel=5.176e-04 time=2.311e+00
[sweep=083] dW_rel=7.194e-04 dTheta_rel=6.636e-04 time=2.623e+00
[sweep=084] dW_rel=7.058e-04 dTheta_rel=4.437e-04 time=2.331e+00
[sweep=085] dW_rel=6.928e-04 dTheta_rel=6.150e-04 time=2.555e+00
[sweep=086] dW_rel=6.802e-04 dTheta_rel=5.009e-04 time=2.471e+00
[sweep=087] dW_rel=6.676e-04 dTheta_rel=5.193e-04 time=2.409e+00
[sweep=088] dW_rel=6.543e-04 dTheta_rel=3.888e-04 time=2.348e+00
[sweep=089] dW_rel=6.440e-04 dTheta_rel=4.915e-04 time=2.480e+00
[sweep=090] dW_rel=6.324e-04 dTheta_rel=5.130e-04 time=2.487e+00
[sweep=091] dW_rel=6.215e-04 dTheta_rel=5.642e-04 time=2.557e+00
[sweep=092] dW_rel=6.094e-04 dTheta_rel=3.800e-04 time=2.401e+00
[sweep=093] dW_rel=6.000e-04 dTheta_rel=5.203e-04 time=2.426e+00
[sweep=094] dW_rel=5.898e-04 dTheta_rel=5.332e-04 time=2.628e+00
[sweep=095] dW_rel=5.804e-04 dTheta_rel=4.591e-04 time=2.414e+00
[sweep=096] dW_rel=5.703e-04 dTheta_rel=3.432e-04 time=2.271e+00
[sweep=097] dW_rel=5.622e-04 dTheta_rel=3.860e-04 time=2.414e+00
[sweep=098] dW_rel=5.528e-04 dTheta_rel=4.820e-04 time=2.416e+00
[sweep=099] dW_rel=5.428e-04 dTheta_rel=4.828e-04 time=2.475e+00
[sweep=100] dW_rel=5.343e-04 dTheta_rel=3.658e-04 time=2.247e+00
[sweep=101] dW_rel=5.259e-04 dTheta_rel=3.697e-04 time=2.167e+00
[sweep=102] dW_rel=5.190e-04 dTheta_rel=4.885e-04 time=2.705e+00
[sweep=103] dW_rel=5.104e-04 dTheta_rel=3.256e-04 time=2.288e+00
[sweep=104] dW_rel=5.037e-04 dTheta_rel=5.953e-04 time=2.779e+00
[sweep=105] dW_rel=4.939e-04 dTheta_rel=4.050e-04 time=2.490e+00
[sweep=106] dW_rel=4.872e-04 dTheta_rel=3.387e-04 time=2.297e+00
[sweep=107] dW_rel=4.804e-04 dTheta_rel=4.519e-04 time=2.481e+00
[sweep=108] dW_rel=4.726e-04 dTheta_rel=4.185e-04 time=2.555e+00
[sweep=109] dW_rel=4.649e-04 dTheta_rel=2.764e-04 time=2.112e+00
[sweep=110] dW_rel=4.597e-04 dTheta_rel=3.905e-04 time=2.472e+00
[sweep=111] dW_rel=4.527e-04 dTheta_rel=3.127e-04 time=2.191e+00
[sweep=112] dW_rel=4.476e-04 dTheta_rel=4.385e-04 time=2.427e+00
[sweep=113] dW_rel=4.378e-04 dTheta_rel=2.666e-04 time=2.113e+00
[sweep=114] dW_rel=4.361e-04 dTheta_rel=4.639e-04 time=2.555e+00
[sweep=115] dW_rel=4.279e-04 dTheta_rel=3.562e-04 time=2.413e+00
[sweep=116] dW_rel=4.220e-04 dTheta_rel=2.351e-04 time=2.113e+00
[sweep=117] dW_rel=4.172e-04 dTheta_rel=4.005e-04 time=2.415e+00
[sweep=118] dW_rel=4.107e-04 dTheta_rel=3.698e-04 time=2.498e+00
[sweep=119] dW_rel=4.054e-04 dTheta_rel=3.294e-04 time=2.205e+00
[sweep=120] split factors 1&2 (cos=0.900)
[sweep=120] dW_rel=5.652e-01 dTheta_rel=2.216e-01 time=7.794e+01
[sweep=121] dW_rel=5.605e-01 dTheta_rel=1.511e-01 time=2.836e+01
[sweep=122] dW_rel=1.895e-01 dTheta_rel=7.754e-02 time=1.355e+01
[sweep=123] dW_rel=9.604e-02 dTheta_rel=1.595e-02 time=8.839e+00
[sweep=124] dW_rel=3.321e-02 dTheta_rel=6.604e-03 time=6.724e+00
[sweep=125] dW_rel=1.847e-02 dTheta_rel=3.499e-03 time=5.211e+00
[sweep=126] dW_rel=1.330e-02 dTheta_rel=2.806e-03 time=5.427e+00
[sweep=127] dW_rel=1.049e-02 dTheta_rel=2.520e-03 time=5.026e+00
[sweep=128] dW_rel=8.571e-03 dTheta_rel=2.200e-03 time=4.955e+00
[sweep=129] dW_rel=7.152e-03 dTheta_rel=1.879e-03 time=4.790e+00
[sweep=130] dW_rel=6.064e-03 dTheta_rel=1.618e-03 time=4.442e+00
[sweep=131] dW_rel=5.202e-03 dTheta_rel=1.399e-03 time=4.320e+00
[sweep=132] dW_rel=4.501e-03 dTheta_rel=1.188e-03 time=4.217e+00
[sweep=133] dW_rel=3.926e-03 dTheta_rel=1.019e-03 time=4.080e+00
[sweep=134] dW_rel=3.444e-03 dTheta_rel=9.106e-04 time=4.053e+00
[sweep=135] dW_rel=3.035e-03 dTheta_rel=7.897e-04 time=4.311e+00
[sweep=136] dW_rel=2.686e-03 dTheta_rel=7.068e-04 time=3.932e+00
[sweep=137] dW_rel=2.383e-03 dTheta_rel=6.076e-04 time=3.723e+00
[sweep=138] dW_rel=2.123e-03 dTheta_rel=5.489e-04 time=3.708e+00
[sweep=139] dW_rel=1.894e-03 dTheta_rel=4.765e-04 time=3.501e+00
[sweep=140] dW_rel=1.695e-03 dTheta_rel=4.640e-04 time=3.789e+00
[sweep=141] dW_rel=1.520e-03 dTheta_rel=3.795e-04 time=3.136e+00
[sweep=142] dW_rel=1.366e-03 dTheta_rel=3.598e-04 time=3.145e+00
[sweep=143] dW_rel=1.229e-03 dTheta_rel=3.059e-04 time=3.015e+00
[sweep=144] dW_rel=1.107e-03 dTheta_rel=2.841e-04 time=2.893e+00
[sweep=145] dW_rel=9.985e-04 dTheta_rel=2.596e-04 time=2.922e+00
[sweep=146] dW_rel=8.997e-04 dTheta_rel=2.326e-04 time=2.831e+00
[sweep=147] dW_rel=8.164e-04 dTheta_rel=2.241e-04 time=2.889e+00
[sweep=148] dW_rel=7.370e-04 dTheta_rel=1.626e-04 time=2.548e+00
[sweep=149] dW_rel=6.672e-04 dTheta_rel=1.733e-04 time=2.684e+00
[sweep=150] dW_rel=6.028e-04 dTheta_rel=1.385e-04 time=2.476e+00
[sweep=151] dW_rel=5.479e-04 dTheta_rel=1.475e-04 time=2.550e+00
[sweep=152] dW_rel=4.964e-04 dTheta_rel=1.387e-04 time=2.586e+00
[sweep=153] dW_rel=4.515e-04 dTheta_rel=1.189e-04 time=2.581e+00
[sweep=154] dW_rel=4.101e-04 dTheta_rel=1.120e-04 time=2.661e+00
[sweep=155] dW_rel=3.721e-04 dTheta_rel=8.214e-05 time=2.250e+00
[sweep=156] dW_rel=3.365e-04 dTheta_rel=9.021e-05 time=2.354e+00
[sweep=157] dW_rel=3.075e-04 dTheta_rel=8.937e-05 time=2.326e+00
[sweep=158] dW_rel=2.812e-04 dTheta_rel=7.797e-05 time=2.425e+00
[sweep=159] dW_rel=2.543e-04 dTheta_rel=5.501e-05 time=2.013e+00
[sweep=160] dW_rel=2.303e-04 dTheta_rel=7.846e-05 time=2.696e+00
[sweep=161] dW_rel=2.095e-04 dTheta_rel=7.636e-05 time=2.391e+00
[sweep=162] dW_rel=1.946e-04 dTheta_rel=7.112e-05 time=2.317e+00
[sweep=163] dW_rel=1.753e-04 dTheta_rel=3.885e-05 time=2.071e+00
[sweep=164] dW_rel=1.591e-04 dTheta_rel=4.328e-05 time=2.251e+00
[sweep=165] dW_rel=1.463e-04 dTheta_rel=5.235e-05 time=2.364e+00
[sweep=166] dW_rel=1.324e-04 dTheta_rel=1.530e-05 time=1.827e+00
[sweep=167] dW_rel=1.230e-04 dTheta_rel=3.756e-05 time=2.089e+00
[sweep=168] dW_rel=1.128e-04 dTheta_rel=3.351e-05 time=2.084e+00
[sweep=169] dW_rel=1.003e-04 dTheta_rel=3.773e-05 time=2.057e+00
[auto_rank prune] p=6 -> 4 (dropped factors [4, 5])
[sweep=001] dW_rel=1.245e-01 dTheta_rel=1.387e-02 time=6.477e+00
[sweep=002] dW_rel=7.621e-02 dTheta_rel=5.089e-03 time=7.603e+00
[sweep=003] dW_rel=1.009e-02 dTheta_rel=3.829e-03 time=4.611e+00
[sweep=004] dW_rel=6.470e-03 dTheta_rel=3.201e-03 time=4.369e+00
[sweep=005] dW_rel=4.908e-03 dTheta_rel=2.254e-03 time=3.824e+00
[sweep=006] dW_rel=4.011e-03 dTheta_rel=1.614e-03 time=3.680e+00
[sweep=007] dW_rel=3.425e-03 dTheta_rel=1.297e-03 time=3.535e+00
[sweep=008] dW_rel=2.996e-03 dTheta_rel=1.078e-03 time=3.277e+00
[sweep=009] dW_rel=2.660e-03 dTheta_rel=8.631e-04 time=3.011e+00
[sweep=010] dW_rel=2.382e-03 dTheta_rel=7.650e-04 time=2.831e+00
[sweep=011] dW_rel=2.143e-03 dTheta_rel=6.405e-04 time=2.622e+00
[sweep=012] dW_rel=1.937e-03 dTheta_rel=6.050e-04 time=2.310e+00
[sweep=013] dW_rel=1.754e-03 dTheta_rel=5.650e-04 time=2.683e+00
[sweep=014] dW_rel=1.596e-03 dTheta_rel=4.686e-04 time=2.381e+00
[sweep=015] dW_rel=1.453e-03 dTheta_rel=4.563e-04 time=2.315e+00
[sweep=016] dW_rel=1.324e-03 dTheta_rel=3.831e-04 time=2.160e+00
[sweep=017] dW_rel=1.206e-03 dTheta_rel=3.610e-04 time=2.238e+00
[sweep=018] dW_rel=1.106e-03 dTheta_rel=3.655e-04 time=1.677e+02
[sweep=019] dW_rel=1.011e-03 dTheta_rel=2.582e-04 time=5.777e+01
[sweep=020] dW_rel=9.239e-04 dTheta_rel=2.773e-04 time=3.450e+00
[sweep=021] dW_rel=8.508e-04 dTheta_rel=2.714e-04 time=4.327e+00
[sweep=022] dW_rel=7.807e-04 dTheta_rel=2.375e-04 time=3.882e+00
[sweep=023] dW_rel=7.143e-04 dTheta_rel=2.226e-04 time=4.503e+00
[sweep=024] dW_rel=6.560e-04 dTheta_rel=1.749e-04 time=3.748e+00
[sweep=025] dW_rel=6.030e-04 dTheta_rel=1.558e-04 time=3.224e+00
[sweep=026] dW_rel=5.560e-04 dTheta_rel=2.150e-04 time=2.844e+00
[sweep=027] dW_rel=5.136e-04 dTheta_rel=1.606e-04 time=3.504e+00
[sweep=028] dW_rel=4.723e-04 dTheta_rel=1.747e-04 time=4.196e+00
[sweep=029] dW_rel=4.299e-04 dTheta_rel=7.535e-05 time=3.386e+00
[sweep=030] dW_rel=4.016e-04 dTheta_rel=1.434e-04 time=4.784e+00
[sweep=031] dW_rel=3.660e-04 dTheta_rel=1.009e-04 time=2.914e+00
[sweep=032] dW_rel=3.428e-04 dTheta_rel=1.087e-04 time=3.636e+00
[sweep=033] dW_rel=3.149e-04 dTheta_rel=1.077e-04 time=1.650e+02
[sweep=034] dW_rel=2.921e-04 dTheta_rel=1.044e-04 time=5.502e+01
[sweep=035] dW_rel=2.667e-04 dTheta_rel=6.739e-05 time=2.047e+00
[sweep=036] dW_rel=2.456e-04 dTheta_rel=6.696e-05 time=1.909e+00
[sweep=037] dW_rel=2.276e-04 dTheta_rel=1.059e-04 time=2.099e+00
[sweep=038] dW_rel=2.101e-04 dTheta_rel=5.481e-05 time=2.098e+00
[sweep=039] dW_rel=1.947e-04 dTheta_rel=9.319e-05 time=2.822e+00
[sweep=040] dW_rel=1.775e-04 dTheta_rel=5.144e-05 time=2.163e+00
[sweep=041] dW_rel=1.677e-04 dTheta_rel=7.898e-05 time=2.664e+00
[sweep=042] dW_rel=1.532e-04 dTheta_rel=3.685e-05 time=2.104e+00
[sweep=043] dW_rel=1.416e-04 dTheta_rel=7.375e-05 time=2.690e+00
[sweep=044] dW_rel=1.311e-04 dTheta_rel=1.755e-05 time=1.753e+00
[sweep=045] dW_rel=1.224e-04 dTheta_rel=4.527e-05 time=2.035e+00
[sweep=046] dW_rel=1.113e-04 dTheta_rel=1.772e-05 time=1.666e+00
[sweep=047] dW_rel=1.039e-04 dTheta_rel=3.846e-05 time=1.914e+00
Runtime: 1375.4s | programs selected: 4
Saved: Results/stgp/ext/stgp_result.pkl
[5]:
ADATA_PATH = OUT_DIR / "adata_with_scores.h5ad"
age_arr = pd.to_numeric(adata_ext.obs["age"], errors="coerce").to_numpy(float)
groups = adata_ext.obs["id_region"].astype(str).to_numpy()
uniq, inv = np.unique(groups, return_inverse=True)
idx_per_group = [np.sort(np.where(inv == t)[0]) for t in range(len(uniq))]
# Apply the same age-ascending sort used during model fitting
_ages_raw = np.array([age_arr[ix[0]] for ix in idx_per_group])
sort_ord = np.argsort(_ages_raw)
idx_sorted = [idx_per_group[i] for i in sort_ord] # cell indices in age order
slices_sorted = uniq[sort_ord] # id_region in age order
ages_sorted = _ages_raw[sort_ord] # ages ascending
adata = adata_ext.copy()
all_idx = np.concatenate(idx_sorted) # res["H"] rows follow this order
H_arr = np.empty_like(res["H"]); H_arr[all_idx] = res["H"]
b_arr = np.empty_like(res["b"]); b_arr[all_idx] = res["b"]
adata.obsm["X_stgp"] = H_arr.astype(np.float32)
adata.obsm["X_stgp_spatial"] = b_arr.astype(np.float32)
adata.uns["stgp"] = dict(
groups=slices_sorted.tolist(), ages=ages_sorted.tolist(),
gamma_age=float(res["gamma_age"]), gamma_spa=float(res["gamma_spa"]),
p_selected=res["W"].shape[0],
alpha=np.asarray(res["alpha"]).tolist(),
alpha_lower=np.asarray(res["alpha_lower"]).tolist(),
alpha_upper=np.asarray(res["alpha_upper"]).tolist(),
theta=np.asarray(res["theta"]).tolist(),
sigma2e=float(res.get("sigma2e", np.nan)),
)
adata.write_h5ad(str(ADATA_PATH), compression="gzip")
print(f"Saved: {ADATA_PATH}")
# Also write W.csv for enrichment
p_sel = res["W"].shape[0]
W_df = pd.DataFrame(res["W"],
index=[f"stGP{j+1}" for j in range(p_sel)],
columns=adata.var_names.astype(str))
W_df.to_csv(OUT_DIR / "W.csv")
Saved: Results/stgp/ext/adata_with_scores.h5ad
[6]:
# ── Reload fitted outputs ──────────────────────────────────────────────────
ADATA_PATH = OUT_DIR / "adata_with_scores.h5ad"
adata = sc.read_h5ad(str(ADATA_PATH))
W_df = pd.read_csv(OUT_DIR / "W.csv", index_col=0)
stgp_info = adata.uns["stgp"]
p_sel = stgp_info["p_selected"]
slices = np.array(stgp_info["groups"])
W_df.index = [f"stGP{i+1}" for i in range(len(W_df))]
print(f"Loaded: {adata.n_obs} cells | {p_sel} programs | {len(slices)} slices")
Loaded: 61162 cells | 4 programs | 12 slices
3. Model Outputs
3.1 Gene Loadings (W matrix)
Each row of W is a non-negative weight vector over the measured genes. Positive-weight genes define the molecular identity of each program; the magnitude reflects the gene’s contribution to that program.
[7]:
W_df = pd.read_csv(OUT_DIR / "W.csv", index_col=0)
W_df.index = [f"stGP{i+1}" for i in range(len(W_df))]
print("Top 10 genes per program:")
for prog, row in W_df.iterrows():
top = row[row > 0].sort_values(ascending=False).head(10)
print(f" {prog}: {', '.join(top.index.tolist())}")
Top 10 genes per program:
stGP1: CBLN2, CUX2, LAMP5, GRIK4, COL19A1, NEUROD1, ONECUT2, SYN3, EPHB1, C1QL3
stGP2: CBLN2, FEZF2, HS3ST4, GRIK4, COL19A1, MOG, PDZRN4, NEUROD6, TBR1, SORCS3
stGP3: RORB, NEUROD6, PLCH1, CNTN5, ZMAT4, CUX2, SORCS2, PVALB, FEZF2, FAM241B
stGP4: AP1G2, NOXA1, HSF4, PDIA2, CORO6, RGS11, NEIL1, NPM2, TMEM145, CHRD
3.2 Spatial Gene-Program Maps
The spatial field b captures the within-slice smooth variation of each program. We tile all tissue sections ordered by donor age so any age-related spatial patterns become visible.
[8]:
# CUX2: canonical L2/3 excitatory-neuron marker – used below as a reference layer
sc.pl.embedding(adata[adata.obs['age']==28].copy(), basis = 'spatial', color = 'CUX2', vmax = 10)
[9]:
from plots import plot_stgp_spatial_programs
scores_df = pd.DataFrame(
adata.obsm["X_stgp_spatial"],
index=adata.obs_names,
columns=[f"stGP{j+1}" for j in range(p_sel)],
)
figs = plot_stgp_spatial_programs(
stgp_adata=adata, scores=scores_df,
age_unit="years",
ncols=4, fg_dot_size=5.0,
)
for j, fig in enumerate(figs):
fig.savefig(FIGURES_DIR / f"spatial_stGP{j+1}.png", dpi=300, bbox_inches="tight")
plt.close(fig)
figs[0]
[9]:
3.3 Age Trajectories (α)
α(t) is the posterior mean age effect of each program — it quantifies how the program amplitude changes across the human lifespan (15–87 yr). The shaded band shows the 95% posterior credible interval.
[10]:
plot_alpha_trajectory(
adata.uns["stgp"],
2,
out_dir=FIGURES_DIR,
stem="alpha_trajectory_stGP3",
);
3.4 Variance Components
[11]:
theta = np.array(stgp_info['theta'])
prog_names = [f'stGP{j+1}' for j in range(p_sel)]
prog_colors = plt.cm.tab10.colors[:p_sel]
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), constrained_layout=True)
for j, (col, name) in enumerate(zip(prog_colors, prog_names)):
axes[0].bar(j, theta[j, 0], color=col, edgecolor='white', linewidth=0.6)
axes[1].bar(j, theta[j, 1], color=col, edgecolor='white', linewidth=0.6)
for ax in axes:
ax.set_xticks(range(p_sel))
ax.set_xticklabels(prog_names, rotation=30, ha='right')
ax.set_xlabel('Program')
axes[0].set_ylabel(r'$\sigma_{\mathrm{age}}^2$')
axes[0].set_title(r'Temporal variance component ($\sigma_{\mathrm{age}}^2$)')
axes[1].set_ylabel(r'$\tau_{\mathrm{spa}}^2$')
axes[1].set_title(r'Spatial variance component ($\tau_{\mathrm{spa}}^2$)')
fig.savefig(FIGURES_DIR / "variance_components.png", dpi=400, bbox_inches="tight")
3.5 Spatial Visualisation of a Single Slice
We inspect one tissue slice in detail, showing the spatial b field of each program (the smooth within-slice component).
[12]:
slice_ages = [(adata.obs.loc[adata.obs['id_region'] == sid, 'age'].iloc[0], sid)
for sid in adata.obs['id_region'].unique()]
slice_ages.sort()
_, example_slice = slice_ages[len(slice_ages) // 2]
sub = adata[adata.obs['id_region'].astype(str) == example_slice].copy()
age_val = sub.obs['age'].iloc[0]
fig, axes = plt.subplots(1, p_sel, figsize=(4.5 * p_sel, 4.5), constrained_layout=True)
b = sub.obsm['X_stgp_spatial']
xy = np.asarray(sub.obsm['spatial'])
for j, ax in enumerate(np.atleast_1d(axes)):
v99 = np.nanpercentile(np.abs(b[:, j]), 99)
sc_ref = ax.scatter(xy[:, 0], xy[:, 1], c=b[:, j],
cmap='RdBu_r', vmin=-v99, vmax=v99,
s=10, linewidths=0, rasterized=True)
ax.set_aspect('equal'); ax.axis('off')
ax.set_title(f'stGP{j+1}')
plt.colorbar(sc_ref, ax=ax, shrink=0.7, pad=0.01)
fig.savefig(FIGURES_DIR / f"spatial_b_{example_slice}.png", dpi=400, bbox_inches="tight")
4. Benchmarking Analysis
The benchmarking logic is consolidated in benchmarking_ext.py so this notebook stays readable and the exported figure/source-data layout is generated from one maintained implementation.
This section compares stGP against STAMP, MEFISTO, Popari, and SpatialPCA using three complementary views:
marker-gene/program correlations for CUX2, RORB, and HS3ST4;
spatial embedding panels and representative slice panels;
clustering recovery against marker-derived layer labels and high-resolution
celltype2layer labels.
All outputs are written under Figure/ext/benchmark.
[13]:
from benchmarking_ext import BASELINE_OBSM_KEYS, LAYER_MARKERS, LAYER_SAFE, METHODS, run_ext_benchmarking
OUT_DIR = RESULTS_DIR / CELLTYPE
ADATA_PATH = OUT_DIR / "adata_with_scores.h5ad"
BL_DIR = Path("Results/baselines")
baseline_adatas = {
"STAMP": sc.read_h5ad(BL_DIR / "stamp_k=3/ext/adata_with_scores.h5ad"),
"MEFISTO": sc.read_h5ad(BL_DIR / "mefisto/ext/adata_with_scores.h5ad"),
"Popari": sc.read_h5ad(BL_DIR / "popari/ext/res_popari.h5ad"),
"SpatialPCA": sc.read_h5ad(BL_DIR / "spatialpca/ext/adata_with_scores.h5ad"),
}
for method, obsm_key in BASELINE_OBSM_KEYS.items():
if obsm_key not in baseline_adatas[method].obsm:
raise KeyError(f"{method} is missing .obsm[{obsm_key!r}]")
if baseline_adatas[method].n_obs != adata.n_obs:
raise ValueError(f"{method} has {baseline_adatas[method].n_obs} cells; expected {adata.n_obs}")
input_summary = pd.DataFrame(
[
{"dataset": "stGP", "cells": adata.n_obs, "slices": len(slices), "embedding": "X_stgp_spatial"},
*[
{
"dataset": method,
"cells": baseline_adatas[method].n_obs,
"slices": baseline_adatas[method].obs["id_region"].nunique(),
"embedding": BASELINE_OBSM_KEYS[method],
}
for method in BASELINE_OBSM_KEYS
],
]
)
display(input_summary)
| dataset | cells | slices | embedding | |
|---|---|---|---|---|
| 0 | stGP | 61162 | 12 | X_stgp_spatial |
| 1 | STAMP | 61162 | 12 | X_stamp |
| 2 | MEFISTO | 61162 | 12 | X_mefisto |
| 3 | Popari | 61162 | 12 | X |
| 4 | SpatialPCA | 61162 | 12 | X_spatialpca |
[14]:
BENCHMARK_DIR = FIGURES_DIR / "benchmark"
benchmark_outputs = run_ext_benchmarking(
adata=adata,
adata_prep=adata_prep,
baseline_adatas=baseline_adatas,
benchmark_dir=BENCHMARK_DIR,
slices=slices,
methods=METHODS,
dpi=400,
)
output_counts = pd.Series(
{
"correlation_figures": len(benchmark_outputs["correlation_figures"]),
"spatial_figures": len(benchmark_outputs["spatial_figures"]),
"cluster_figures": len(benchmark_outputs["cluster_figures"]),
"metric_figures": len(benchmark_outputs["metric_figures"]),
"summary_figures": len(benchmark_outputs["summary_figures"]),
},
name="n_outputs",
)
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
/home/byual/.conda/envs/stGP/lib/python3.11/site-packages/sklearn/manifold/_spectral_embedding.py:324: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
warnings.warn(
[15]:
import os
from benchmarking_ext import (
BASELINE_OBSM_KEYS,
CELLTYPE2_LAYER_MAP,
LAYER_COLORS,
LAYER_MARKERS,
METHODS,
_representative_slices_by_age,
_slice_ages,
_slice_order,
)
from plots import (
plot_alpha_trajectory,
plot_embedding_marker_correlation_boxplots,
plot_human_cluster_embedding_composite,
plot_layer_cluster_tiles,
plot_marker_region_ari_nmi_combined,
plot_representative_continuous_tiles,
plot_representative_embedding_tiles,
)
from utils import (
as_1d_array,
best_program_by_correlation,
clean_layer_labels,
one_to_one_series_from_prediction_csv,
)
PROJECT_ROOT = Path(os.environ.get('STGP_REPRO_ROOT', PROJECT_DIR.parent))
HUMAN_DIR = PROJECT_ROOT / 'RealData_HumanBrainMERFISH'
BENCHMARK_RESULTS_DIR = HUMAN_DIR / 'Results' / 'benchmark' / CELLTYPE
stale_output_stems = [
'clust_MEFISTO', 'clust_Popari', 'clust_STAMP', 'clust_SpatialPCA',
'MEFISTO', 'Popari', 'STAMP', 'SpatialPCA',
'human_marker_region_ARI', 'human_marker_region_NMI',
'human_L4_embedding_vs_marker_correlation',
]
[16]:
adata = sc.read_h5ad(HUMAN_DIR / 'Results' / 'stgp' / CELLTYPE / 'adata_with_scores.h5ad')
adata_log = sc.pp.log1p(adata.copy(), copy=True)
baseline_paths = {
'STAMP': HUMAN_DIR / 'Results' / 'baselines' / 'stamp_k=3' / CELLTYPE / 'adata_with_scores.h5ad',
'MEFISTO': HUMAN_DIR / 'Results' / 'baselines' / 'mefisto' / CELLTYPE / 'adata_with_scores.h5ad',
'Popari': HUMAN_DIR / 'Results' / 'baselines' / 'popari' / CELLTYPE / 'res_popari.h5ad',
'SpatialPCA': HUMAN_DIR / 'Results' / 'baselines' / 'spatialpca' / CELLTYPE / 'adata_with_scores.h5ad',
}
baseline_adatas = {method: sc.read_h5ad(path) for method, path in baseline_paths.items()}
baseline_adatas = {
method: bl if bl.obs_names.equals(adata.obs_names) else bl[adata.obs_names].copy()
for method, bl in baseline_adatas.items()
}
slices_sorted = _slice_order(adata)
rep_slices = _representative_slices_by_age(adata, slices_sorted)
rep_slices = ['5657_rep1' if sid == '5657_rep0' else sid for sid in rep_slices]
age_by_slice = _slice_ages(adata, slices_sorted)
ids = adata.obs['id_region'].astype(str).to_numpy()
obs_index = pd.Index(adata.obs_names.astype(str))
# Fig. 3 uses the L4 marker RORB as the representative embedding benchmark.
layer = 'L4'
gene = LAYER_MARKERS[layer]
gene_idx = np.where(adata_log.var_names == gene)[0][0]
rorb_log = as_1d_array(adata_log.X[:, gene_idx])
bl_k_l4 = {
method: best_program_by_correlation(baseline_adatas[method], BASELINE_OBSM_KEYS[method], rorb_log)
for method in ['MEFISTO', 'Popari', 'STAMP', 'SpatialPCA']
}
baseline_embedding_values = {
method: np.asarray(baseline_adatas[method].obsm[BASELINE_OBSM_KEYS[method]])[:, k]
for method, k in bl_k_l4.items()
}
celltype2 = adata.obs['celltype2'].astype('string')
celltype_values = clean_layer_labels(celltype2.map(CELLTYPE2_LAYER_MAP).fillna('ext'))
cluster_palette = {
'L2/3': LAYER_COLORS['L2/3'],
'L4': LAYER_COLORS['L4'],
'L5/6': LAYER_COLORS['L5/6'],
'ext': LAYER_COLORS.get('ext', '#BFBFBF'),
}
legend_items = [('L2/3', 'L2/3'), ('L4', 'L4'), ('L5/6', 'L5/6'), ('ext', 'ext')]
marker_pred_csv = BENCHMARK_RESULTS_DIR / 'clustering' / 'marker_region' / 'cell_predictions.csv'
cell_pred_csv = BENCHMARK_RESULTS_DIR / 'clustering' / 'celltype2' / 'cell_predictions.csv'
cluster_values = {}
for method in METHODS:
marker_series = one_to_one_series_from_prediction_csv(marker_pred_csv, method, obs_index=obs_index, default='ext')
cell_series = one_to_one_series_from_prediction_csv(cell_pred_csv, method, obs_index=obs_index, default=np.nan)
combined = marker_series.copy()
combined.loc[cell_series.dropna().index] = cell_series.dropna()
cluster_values[method] = clean_layer_labels(combined)
fig3_summary = pd.DataFrame({
'representative_id_region': rep_slices,
'age_years': [age_by_slice[s] for s in rep_slices],
'n_cells': [int((ids == s).sum()) for s in rep_slices],
})
display(fig3_summary)
print('Best L4/RORB programs, 1-based:', {m: k + 1 for m, k in bl_k_l4.items()})
| representative_id_region | age_years | n_cells | |
|---|---|---|---|
| 0 | 6052_rep1 | 28.0 | 5550 |
| 1 | 4643_rep0 | 42.0 | 10002 |
| 2 | 5657_rep1 | 82.0 | 3267 |
| 3 | 4320_rep0 | 87.0 | 4213 |
Best L4/RORB programs, 1-based: {'MEFISTO': 2, 'Popari': 4, 'STAMP': 2, 'SpatialPCA': 4}
4.1 stGP and layer marker in representative slices
[17]:
# Representative label panels: annotated layer labels and stGP-predicted layers.
plot_layer_cluster_tiles(
adata,
celltype_values,
rep_slices=rep_slices,
palette=cluster_palette,
legend_items=legend_items,
out_dir=FIGURES_DIR,
stem='celltype2',
title='celltype',
);
[18]:
plot_layer_cluster_tiles(
adata,
cluster_values['stGP'],
rep_slices=rep_slices,
palette=cluster_palette,
legend_items=legend_items,
out_dir=FIGURES_DIR,
stem='stGP',
title='stGP',
);
[19]:
# The L4 marker RORB and the matching stGP spatial field are shown on the same sections.
plot_representative_continuous_tiles(
adata,
rorb_log,
rep_slices=rep_slices,
cmap='YlOrBr',
symmetric=False,
vmin=0.0,
vmax=4.0,
colorbar_label='log1p expression',
title='RORB',
out_dir=FIGURES_DIR,
stem='RORB_expression',
);
[20]:
plot_representative_embedding_tiles(
adata,
np.asarray(adata.obsm['X_stgp_spatial'])[:, 2],
rep_slices=rep_slices,
out_dir=FIGURES_DIR,
stem='stGP_b',
title='stGP',
signed=True,
);
4.2 Baseline Cluster/Embedding Composite
For each baseline method, the left mini-panel shows layer-like clusters and the right mini-panel shows the embedding component most correlated with RORB.
[21]:
method_order = ['SpatialPCA', 'Popari', 'MEFISTO', 'STAMP']
plot_human_cluster_embedding_composite(
adata,
rep_slices=rep_slices,
cluster_values=cluster_values,
embedding_values=baseline_embedding_values,
methods=method_order,
palette=cluster_palette,
legend_items=legend_items,
out_dir=FIGURES_DIR,
);
4.3 Benchmark Summary Panels
These panels use the id_region-keyed benchmark source tables generated above: marker-region ARI/NMI for cluster recovery, and layer marker correlations for embedding quality.
[22]:
marker_metrics_csv = BENCHMARK_RESULTS_DIR / 'clustering' / 'marker_region' / 'slice_metrics.csv'
correlation_csv = BENCHMARK_RESULTS_DIR / 'summary' / 'source_data' / 'marker_embedding_correlations.csv'
plot_marker_region_ari_nmi_combined(
marker_metrics_csv,
methods=METHODS,
out_dir=FIGURES_DIR,
);
[23]:
plot_embedding_marker_correlation_boxplots(
correlation_csv,
methods=METHODS,
layer_markers=LAYER_MARKERS,
out_dir=FIGURES_DIR,
);
[ ]: