Excitatory Neuron Gene Programs in the Aging Human Brain

This tutorial walks through the full stGP pipeline on the human 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. stGP decomposes gene expression in each slice into p latent programs—each with a non-negative gene loading vector (W), a cell-level activity score (H), a spatially smooth field (b), and an age-varying amplitude (α)—enabling joint discovery of spatial and age-related transcriptional programs in the aging brain.

1. Setup

[1]:

import sys, warnings, pickle
from pathlib import Path

import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt
from IPython.display import display

from plots import set_nature_style

set_nature_style()
warnings.filterwarnings("ignore", category=FutureWarning)
sys.path.insert(0, "..")

# Paths
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 for this tutorial
adata_ext = sc.read_h5ad(DATA_PROC / f"{CELLTYPE}.h5ad")

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 program

  • H (cell scores, N × p): overall activity of each program in each cell

  • b (spatial field, N × p): spatially smooth residual component

  • α (age effect, p × S): how program amplitude varies across slices/ages

  • θ (GP hyperparameters, p × 2): spatial amplitude and noise fraction per 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]:

if PKL_PATH.exists():
    with open(PKL_PATH, "rb") as f:
        res = pickle.load(f)
    print(f"Loaded: {PKL_PATH}")
else:
    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]}")

    # ── Save results ─────────────────────────────────────────────────────────
    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}")

Loaded: 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

[8]:

from plots import plot_W_program_heatmap

fig = plot_W_program_heatmap(
    W_df,
    out=FIGURES_DIR / "W_heatmap_vertical.png",
    dpi=400,
    orientation="vertical",
)

[9]:

from plots import plot_W_program_heatmap

fig = plot_W_program_heatmap(
    W_df,
    out=FIGURES_DIR / "W_heatmap.png",
    dpi=400, orientation = "horizontal"
)

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.

[10]:

# 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)

[11]:

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,
    celltype="Excitatory Neuron", age_unit="years",
    ncols=4, fg_dot_size=5.0, dpi=300,
)
for j, fig in enumerate(figs):
    fig.savefig(FIGURES_DIR / f"spatial_stGP{j+1}.png", dpi=300, bbox_inches="tight")
    plt.close(fig)

[12]:

figs[0]

[12]:

<Figure size 3120x2310 with 13 Axes>

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.

[13]:

stgp_info   = adata.uns['stgp']
ages_slices = np.array(stgp_info['ages'])
alpha       = np.array(stgp_info['alpha'])        # (p, n_slices)
alpha_lower = np.array(stgp_info['alpha_lower'])
alpha_upper = np.array(stgp_info['alpha_upper'])

COLOR = '#2C7FB8'
order = np.argsort(ages_slices)
t = ages_slices[order]

for j in range(p_sel):
    a  = alpha[j][order]
    lo = alpha_lower[j][order]
    hi = alpha_upper[j][order]

    fig, ax = plt.subplots(figsize=(4.5, 4), constrained_layout=True)

    ax.fill_between(t, lo, hi, alpha=0.18, color=COLOR)
    ax.plot(t, lo, lw=0.8, ls='--', color=COLOR, alpha=0.55)
    ax.plot(t, hi, lw=0.8, ls='--', color=COLOR, alpha=0.55)
    ax.plot(t, a, lw=1.8, color=COLOR)
    ax.scatter(t, a, s=32, color=COLOR, zorder=3, label='Posterior mean')
    ax.axhline(0, color='0.6', lw=0.7, ls=':')

    ax.set_xlabel('Age (yr)')
    ax.set_ylabel('Age effect α')
    ax.legend(fontsize=9)
    fig.savefig(FIGURES_DIR / f"alpha_trajectory_stGP{j+1}.png", dpi=400, bbox_inches="tight")
    plt.close(fig)

3.4 GP Parameters (θ)

[14]:

theta = np.asarray(stgp_info["theta"], dtype=float)
theta

[14]:

array([[ 38.81529414,  15.40473634],
       [973.40129965,  12.72626063],
       [223.0523848 ,  11.20525559],
       [ 84.11608287, 147.19819587]])

[15]:

theta      = np.array(stgp_info['theta'])   # (p, 2): [amplitude, noise_frac]
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).

[16]:

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 celltype2 layer labels.

All outputs are written under Figure/ext/benchmark.

[17]:

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"

if "adata" not in globals() or "X_stgp_spatial" not in adata.obsm:
    adata = sc.read_h5ad(ADATA_PATH)

if "adata_prep" not in globals():
    adata_temp = adata_ext.copy()
    sc.pp.normalize_total(adata_temp, target_sum=1e3)
    adata_prep = sc.pp.log1p(adata_temp, copy=True)

if "slices" not in globals():
    if "stgp" in adata.uns and "groups" in adata.uns["stgp"]:
        slices = np.array(adata.uns["stgp"]["groups"])
    else:
        slice_ids = adata.obs["id_region"].astype(str)
        slices = np.array(
            sorted(
                pd.unique(slice_ids),
                key=lambda sid: float(adata.obs.loc[slice_ids == sid, "age"].iloc[0]),
            )
        )

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 
[18]:

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",
)

print(f"Benchmark outputs saved under: {benchmark_outputs['benchmark_dir']}")
display(output_counts.to_frame())

/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(

Benchmark outputs saved under: Figure/ext/benchmark
n_outputs
correlation_figures 4
spatial_figures 24
cluster_figures 36
metric_figures 10
summary_figures 13 
[19]:

source_dir = Path(benchmark_outputs["source_dir"]) / "summary" / "source_data"

corr_summary = pd.read_csv(source_dir / "marker_embedding_correlation_summary.csv")
display(corr_summary)

for ground_truth, metrics_df in benchmark_outputs["cluster_metrics"].items():
    display(
        metrics_df
        .groupby("method", as_index=False)
        .agg(
            n_slices=("id_region", "nunique"),
            raw_ari=("raw_ari", "mean"),
            raw_nmi=("raw_nmi", "mean"),
            raw_acc=("raw_acc", "mean"),
        )
        .sort_values("raw_ari", ascending=False)
    )
Unnamed: 0 layer marker_gene method mean median std
0 0 L2/3 CUX2 MEFISTO 0.749311 0.815860 0.207515
1 1 L2/3 CUX2 Popari 0.362930 0.248515 0.239036
2 2 L2/3 CUX2 STAMP 0.610371 0.737565 0.281830
3 3 L2/3 CUX2 SpatialPCA 0.545955 0.463947 0.222554
4 4 L2/3 CUX2 stGP 0.796851 0.838556 0.155620
5 5 L4 RORB MEFISTO 0.604304 0.612808 0.096769
6 6 L4 RORB Popari 0.283265 0.242078 0.265295
7 7 L4 RORB STAMP 0.371638 0.342678 0.179853
8 8 L4 RORB SpatialPCA 0.257941 0.301475 0.158959
9 9 L4 RORB stGP 0.674625 0.699595 0.097605
10 10 L5/6 HS3ST4 MEFISTO 0.447195 0.468178 0.102479
11 11 L5/6 HS3ST4 Popari 0.125816 0.056475 0.137783
12 12 L5/6 HS3ST4 STAMP 0.227008 0.334305 0.238527
13 13 L5/6 HS3ST4 SpatialPCA 0.325266 0.323195 0.140685
14 14 L5/6 HS3ST4 stGP 0.515356 0.551536 0.101769
method n_slices raw_ari raw_nmi raw_acc
4 stGP 12 0.322710 0.311816 0.687081
0 MEFISTO 12 0.272890 0.276835 0.661968
2 STAMP 12 0.251136 0.225701 0.626740
3 SpatialPCA 12 0.148898 0.169466 0.543434
1 Popari 12 0.107890 0.109065 0.502050
method n_slices raw_ari raw_nmi raw_acc
4 stGP 9 0.588080 0.550121 0.829895
0 MEFISTO 9 0.565688 0.531591 0.820396
2 STAMP 9 0.407057 0.389712 0.729004
3 SpatialPCA 9 0.233973 0.262811 0.606660
1 Popari 9 0.179829 0.164645 0.540688

5. Pathway Enrichment Analysis

We use gseapy to run over-representation analysis (ORA) on the positive-weight genes of each stGP program, testing against GO Biological Process and GO Cellular Component gene sets. For excitatory neurons, we expect programs enriched in synaptic transmission, neuronal differentiation, axon guidance, and age-related processes such as protein folding stress and mitochondrial dysfunction.

Prerequisites: Download MSigDB gene-set files from https://www.gsea-msigdb.org/ and place them under data/genesets/:

  • c5.go.bp.v2026.1.Hs.symbols.gmt

  • c5.go.cc.v2026.1.Hs.symbols.gmt

[20]:

from plots import plot_enrichment_dotplot, truncate_colormap as _trunc_cmap
import gseapy as gp

gene_sets = {
    "GO Biological Process": "data/genesets/c5.go.bp.v2026.1.Hs.symbols.gmt",
    "GO Cellular Component": "data/genesets/c5.go.cc.v2026.1.Hs.symbols.gmt",
}
cmaps = {
    "GO Biological Process": _trunc_cmap("Reds",    0.22, 0.78),
    "GO Cellular Component": _trunc_cmap("Purples", 0.25, 0.76),
}
background_genes = list(W_df.columns)

all_res = []
for program in W_df.index:
    gene_list = W_df.loc[program].pipe(lambda s: s[s > 0].sort_values(ascending=False).index.tolist())
    print(f"{program}: {len(gene_list)} active genes")

    fig, axes = plt.subplots(len(gene_sets), 1,
                              figsize=(7.5, 2.8 * len(gene_sets)), constrained_layout=True)
    for ax, (set_name, gmt_file) in zip(np.atleast_1d(axes), gene_sets.items()):
        enr = gp.enrich(gene_list=gene_list, gene_sets=gmt_file,
                        background=background_genes, verbose=False)
        res = enr.res2d.copy()
        res["program"] = program
        all_res.append(res)
        plot_enrichment_dotplot(res, ax, set_name, cmap=cmaps[set_name])

    fig.savefig(FIGURES_DIR / f"{program}_enrichment.png", dpi=300, bbox_inches="tight")
    plt.show()

enr_df = pd.concat(all_res, ignore_index=True)
enr_df.to_csv(FIGURES_DIR / "enrichment_results.csv", index=False)
print(f"\nAll enrichment results saved to {FIGURES_DIR}/enrichment_results.csv")

stGP1: 15 active genes
stGP2: 15 active genes
stGP3: 15 active genes
stGP4: 15 active genes

All enrichment results saved to Figure/ext/enrichment_results.csv

[21]:

all_res[1][all_res[1]['Adjusted P-value']<0.05]

[21]:
Gene_set Term Overlap P-value Adjusted P-value Odds Ratio Combined Score Genes program
41 c5.go.cc.v2026.1.Hs.symbols.gmt GOCC_GLUTAMATERGIC_SYNAPSE 6/17 0.000064 0.006870 15.439359 149.042913 C1QL3;SYN3;GRIK4;APOE;CBLN2;EPHB1 stGP1
91 c5.go.cc.v2026.1.Hs.symbols.gmt GOCC_SYNAPSE 8/45 0.000583 0.020805 7.056889 52.551218 C1QL3;SV2C;LAMP5;SYN3;GRIK4;APOE;CBLN2;EPHB1 stGP1
92 c5.go.cc.v2026.1.Hs.symbols.gmt GOCC_SYNAPTIC_CLEFT 3/4 0.000462 0.020805 50.306667 386.405222 C1QL3;CBLN2;APOE stGP1 
[ ]: