{ "cells": [ { "cell_type": "markdown", "id": "md-title", "metadata": {}, "source": [ "# Excitatory Neuron Gene Programs in the Aging Human Brain\n", "\n", "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.\n", "\n", "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.\n" ] }, { "cell_type": "markdown", "id": "md-sec1", "metadata": {}, "source": [ "## 1. Setup" ] }, { "cell_type": "code", "execution_count": 1, "id": "setup", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:17:35.936231Z", "iopub.status.busy": "2026-07-04T22:17:35.935979Z", "iopub.status.idle": "2026-07-04T22:17:46.904751Z", "shell.execute_reply": "2026-07-04T22:17:46.903954Z" } }, "outputs": [], "source": [ "import sys, warnings, pickle\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import scanpy as sc\n", "import matplotlib.pyplot as plt\n", "from IPython.display import display\n", "\n", "from plots import set_nature_style\n", "\n", "set_nature_style()\n", "warnings.filterwarnings(\"ignore\", category=FutureWarning)\n", "sys.path.insert(0, \"..\")\n", "\n", "# Paths\n", "RAW_DATA_DIR = Path(\"/import/home2/share/byual/HumanBrainMERFISH+sc_Nature2025_Jeffries\")\n", "DATA_QC = Path(\"data/qc/human_merfish_qc.h5ad\")\n", "DATA_PROC = Path(\"data/processed\")\n", "RESULTS_DIR = Path(\"Results/stgp\")\n", "FIGURES_DIR = Path(\"Figure/ext\")\n", "FIGURES_DIR.mkdir(parents=True, exist_ok=True)\n", "\n", "CELLTYPE = \"ext\" # target cell type for this tutorial\n", "adata_ext = sc.read_h5ad(DATA_PROC / f\"{CELLTYPE}.h5ad\")" ] }, { "cell_type": "markdown", "id": "md-sec3", "metadata": {}, "source": [ "---\n", "## 2. Fitting stGP\n", "\n", "stGP models gene expression in each tissue slice as a sum of **p latent programs**.\n", "Each program is characterised by:\n", "- **W** (gene loadings, `p × G`): non-negative gene weights defining the program\n", "- **H** (cell scores, `N × p`): overall activity of each program in each cell\n", "- **b** (spatial field, `N × p`): spatially smooth residual component\n", "- **α** (age effect, `p × S`): how program amplitude varies across slices/ages\n", "- **θ** (GP hyperparameters, `p × 2`): spatial amplitude and noise fraction per program\n", "\n", "The model rank **p** is selected automatically by greedy forward selection." ] }, { "cell_type": "code", "execution_count": 2, "id": "run-stgp", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:17:46.908103Z", "iopub.status.busy": "2026-07-04T22:17:46.907743Z", "iopub.status.idle": "2026-07-04T22:17:51.483174Z", "shell.execute_reply": "2026-07-04T22:17:51.481883Z" } }, "outputs": [], "source": [ "import time\n", "from stgp.estimation import fit_pfactor_auto\n", "from stgp.kernels import (\n", " bandwidth_select_spatial, bandwidth_select_temporal,\n", " build_K_age, build_K_spa_list_from_stacked\n", ")\n", "from stgp.preprocessing import standardize_coords_list, log1p_norm_centered_list\n", "\n", "OUT_DIR = RESULTS_DIR / CELLTYPE\n", "OUT_DIR.mkdir(parents=True, exist_ok=True)\n", "PKL_PATH = OUT_DIR / \"stgp_result.pkl\"\n", "\n", "age_arr = pd.to_numeric(adata_ext.obs[\"age\"], errors=\"coerce\").to_numpy(float)\n", "groups = adata_ext.obs[\"id_region\"].astype(str).to_numpy()\n", "uniq, inv = np.unique(groups, return_inverse=True)\n", "idx_per_group = [np.sort(np.where(inv == t)[0]) for t in range(len(uniq))]\n", "\n", "adata_temp = adata_ext.copy()\n", "sc.pp.normalize_total(adata_temp, target_sum=1e3)\n", "adata_prep = sc.pp.log1p(adata_temp, copy = True)\n", "\n", "Y_list = [adata_prep.X[ix].toarray() for ix in idx_per_group]\n", "Y_list, _ = log1p_norm_centered_list(Y_list, target_sum = 1000)\n", "nlist = np.array([len(ix) for ix in idx_per_group])\n", "ages = np.array([age_arr[ix[0]] for ix in idx_per_group])\n", "sort_ord = np.argsort(ages); ages = ages[sort_ord]\n", "slices = uniq.copy(); slices = slices[sort_ord]\n", "nlist = nlist[sort_ord]\n", "Y_list = [Y_list[i] for i in sort_ord]" ] }, { "cell_type": "code", "execution_count": 3, "id": "ddd48513", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:17:51.485845Z", "iopub.status.busy": "2026-07-04T22:17:51.485596Z", "iopub.status.idle": "2026-07-04T22:17:58.813467Z", "shell.execute_reply": "2026-07-04T22:17:58.811017Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " gamma_spa = 0.1473 | gamma_age = 1.0997\n" ] } ], "source": [ "# ── Build GP kernels ─────────────────────────────────────────────────────\n", "coords_list = standardize_coords_list([adata_ext.obsm[\"spatial\"][ix] for ix in idx_per_group])\n", "coords_list = [coords_list[i] for i in sort_ord]\n", "gamma_spa = bandwidth_select_spatial(coords_list, frac=0.01, rho=0.6)\n", "gamma_age = bandwidth_select_temporal(ages, rho=np.exp(-1.5))\n", "print(f\" gamma_spa = {gamma_spa:.4f} | gamma_age = {gamma_age:.4f}\")\n", "\n", "K_age = build_K_age(ages, gamma_age, kernel=\"rbf\", standardize=True)\n", "K_spa_list = build_K_spa_list_from_stacked(\n", " np.vstack(coords_list), nlist, gamma_spa, standardize=False, jitter=1e-6\n", ")" ] }, { "cell_type": "code", "execution_count": 4, "id": "ed1ea14e", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:17:58.818049Z", "iopub.status.busy": "2026-07-04T22:17:58.817768Z", "iopub.status.idle": "2026-07-04T22:17:58.844774Z", "shell.execute_reply": "2026-07-04T22:17:58.843787Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Loaded: Results/stgp/ext/stgp_result.pkl\n" ] } ], "source": [ "if PKL_PATH.exists():\n", " with open(PKL_PATH, \"rb\") as f:\n", " res = pickle.load(f)\n", " print(f\"Loaded: {PKL_PATH}\")\n", "else:\n", " t0 = time.perf_counter()\n", " res = fit_pfactor_auto(\n", " Y_list=Y_list, Nlist=nlist, K_age=K_age, Kspa_list=K_spa_list,\n", " p_max=10, k=15,\n", " inner_rank1_tol=1e-4, rel_improve_total_tol=0.002, backfit_tol=1e-4, prune_energy_frac = 0.005,\n", " random_state=0, verbose=1,\n", " )\n", " print(f\"Runtime: {time.perf_counter() - t0:.1f}s | programs selected: {res['W'].shape[0]}\")\n", "\n", " # ── Save results ─────────────────────────────────────────────────────────\n", " res[\"gamma_age\"] = gamma_age; res[\"gamma_spa\"] = gamma_spa\n", " with open(PKL_PATH, \"wb\") as f:\n", " pickle.dump(res, f)\n", " print(f\"Saved: {PKL_PATH}\")" ] }, { "cell_type": "code", "execution_count": 5, "id": "attach-scores", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:17:58.846563Z", "iopub.status.busy": "2026-07-04T22:17:58.846421Z", "iopub.status.idle": "2026-07-04T22:18:04.316843Z", "shell.execute_reply": "2026-07-04T22:18:04.315247Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saved: Results/stgp/ext/adata_with_scores.h5ad\n" ] } ], "source": [ "ADATA_PATH = OUT_DIR / \"adata_with_scores.h5ad\"\n", "\n", "age_arr = pd.to_numeric(adata_ext.obs[\"age\"], errors=\"coerce\").to_numpy(float)\n", "groups = adata_ext.obs[\"id_region\"].astype(str).to_numpy()\n", "uniq, inv = np.unique(groups, return_inverse=True)\n", "idx_per_group = [np.sort(np.where(inv == t)[0]) for t in range(len(uniq))]\n", "\n", "# Apply the same age-ascending sort used during model fitting\n", "_ages_raw = np.array([age_arr[ix[0]] for ix in idx_per_group])\n", "sort_ord = np.argsort(_ages_raw)\n", "idx_sorted = [idx_per_group[i] for i in sort_ord] # cell indices in age order\n", "slices_sorted = uniq[sort_ord] # id_region in age order\n", "ages_sorted = _ages_raw[sort_ord] # ages ascending\n", "\n", "adata = adata_ext.copy()\n", "all_idx = np.concatenate(idx_sorted) # res[\"H\"] rows follow this order\n", "H_arr = np.empty_like(res[\"H\"]); H_arr[all_idx] = res[\"H\"]\n", "b_arr = np.empty_like(res[\"b\"]); b_arr[all_idx] = res[\"b\"]\n", "adata.obsm[\"X_stgp\"] = H_arr.astype(np.float32)\n", "adata.obsm[\"X_stgp_spatial\"] = b_arr.astype(np.float32)\n", "adata.uns[\"stgp\"] = dict(\n", " groups=slices_sorted.tolist(), ages=ages_sorted.tolist(),\n", " gamma_age=float(res[\"gamma_age\"]), gamma_spa=float(res[\"gamma_spa\"]),\n", " p_selected=res[\"W\"].shape[0],\n", " alpha=np.asarray(res[\"alpha\"]).tolist(),\n", " alpha_lower=np.asarray(res[\"alpha_lower\"]).tolist(),\n", " alpha_upper=np.asarray(res[\"alpha_upper\"]).tolist(),\n", " theta=np.asarray(res[\"theta\"]).tolist(),\n", " sigma2e=float(res.get(\"sigma2e\", np.nan)),\n", ")\n", "adata.write_h5ad(str(ADATA_PATH), compression=\"gzip\")\n", "print(f\"Saved: {ADATA_PATH}\")\n", "\n", "# Also write W.csv for enrichment\n", "p_sel = res[\"W\"].shape[0]\n", "W_df = pd.DataFrame(res[\"W\"],\n", " index=[f\"stGP{j+1}\" for j in range(p_sel)],\n", " columns=adata.var_names.astype(str))\n", "W_df.to_csv(OUT_DIR / \"W.csv\")" ] }, { "cell_type": "code", "execution_count": 6, "id": "load-results", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:04.320387Z", "iopub.status.busy": "2026-07-04T22:18:04.320063Z", "iopub.status.idle": "2026-07-04T22:18:04.619185Z", "shell.execute_reply": "2026-07-04T22:18:04.618359Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Loaded: 61162 cells | 4 programs | 12 slices\n" ] } ], "source": [ "# ── Reload fitted outputs ──────────────────────────────────────────────────\n", "ADATA_PATH = OUT_DIR / \"adata_with_scores.h5ad\"\n", "adata = sc.read_h5ad(str(ADATA_PATH))\n", "W_df = pd.read_csv(OUT_DIR / \"W.csv\", index_col=0)\n", "stgp_info = adata.uns[\"stgp\"]\n", "p_sel = stgp_info[\"p_selected\"]\n", "slices = np.array(stgp_info[\"groups\"])\n", "W_df.index = [f\"stGP{i+1}\" for i in range(len(W_df))]\n", "print(f\"Loaded: {adata.n_obs} cells | {p_sel} programs | {len(slices)} slices\")" ] }, { "cell_type": "markdown", "id": "md-sec4-new", "metadata": {}, "source": [ "---\n", "## 3. Model Outputs\n", "\n", "### 3.1 Gene Loadings (W matrix)\n", "\n", "Each row of **W** is a non-negative weight vector over the measured genes.\n", "Positive-weight genes define the molecular identity of each program; the\n", "magnitude reflects the gene's contribution to that program.\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "gene-loadings", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:04.621196Z", "iopub.status.busy": "2026-07-04T22:18:04.621065Z", "iopub.status.idle": "2026-07-04T22:18:04.628734Z", "shell.execute_reply": "2026-07-04T22:18:04.628167Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Top 10 genes per program:\n", " stGP1: CBLN2, CUX2, LAMP5, GRIK4, COL19A1, NEUROD1, ONECUT2, SYN3, EPHB1, C1QL3\n", " stGP2: CBLN2, FEZF2, HS3ST4, GRIK4, COL19A1, MOG, PDZRN4, NEUROD6, TBR1, SORCS3\n", " stGP3: RORB, NEUROD6, PLCH1, CNTN5, ZMAT4, CUX2, SORCS2, PVALB, FEZF2, FAM241B\n", " stGP4: AP1G2, NOXA1, HSF4, PDIA2, CORO6, RGS11, NEIL1, NPM2, TMEM145, CHRD\n" ] } ], "source": [ "W_df = pd.read_csv(OUT_DIR / \"W.csv\", index_col=0)\n", "W_df.index = [f\"stGP{i+1}\" for i in range(len(W_df))]\n", "\n", "print(\"Top 10 genes per program:\")\n", "for prog, row in W_df.iterrows():\n", " top = row[row > 0].sort_values(ascending=False).head(10)\n", " print(f\" {prog}: {', '.join(top.index.tolist())}\")" ] }, { "cell_type": "code", "execution_count": 8, "id": "693e5105", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:04.630391Z", "iopub.status.busy": "2026-07-04T22:18:04.630256Z", "iopub.status.idle": "2026-07-04T22:18:05.643615Z", "shell.execute_reply": "2026-07-04T22:18:05.641941Z" } }, "outputs": [], "source": [ "from plots import plot_W_program_heatmap\n", "\n", "fig = plot_W_program_heatmap(\n", " W_df,\n", " out=FIGURES_DIR / \"W_heatmap_vertical.png\",\n", " dpi=400,\n", " orientation=\"vertical\",\n", ")" ] }, { "cell_type": "code", "execution_count": 9, "id": "W-heatmap", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:05.645967Z", "iopub.status.busy": "2026-07-04T22:18:05.645699Z", "iopub.status.idle": "2026-07-04T22:18:06.041782Z", "shell.execute_reply": "2026-07-04T22:18:06.040591Z" } }, "outputs": [], "source": [ "from plots import plot_W_program_heatmap\n", "\n", "fig = plot_W_program_heatmap(\n", " W_df,\n", " out=FIGURES_DIR / \"W_heatmap.png\",\n", " dpi=400, orientation = \"horizontal\"\n", ")" ] }, { "cell_type": "markdown", "id": "md-sec4-2", "metadata": {}, "source": [ "### 3.2 Spatial Gene-Program Maps\n", "\n", "The **spatial field b** captures the within-slice smooth variation of each program.\n", "We tile all tissue sections ordered by donor age so any age-related spatial patterns\n", "become visible." ] }, { "cell_type": "code", "execution_count": 10, "id": "add206d0", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:06.044010Z", "iopub.status.busy": "2026-07-04T22:18:06.043871Z", "iopub.status.idle": "2026-07-04T22:18:06.083739Z", "shell.execute_reply": "2026-07-04T22:18:06.082845Z" } }, "outputs": [], "source": [ "# CUX2: canonical L2/3 excitatory-neuron marker – used below as a reference layer\n", "sc.pl.embedding(adata[adata.obs['age']==28].copy(), basis = 'spatial', color = 'CUX2', vmax = 10)" ] }, { "cell_type": "code", "execution_count": 11, "id": "spatial-maps", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:06.085848Z", "iopub.status.busy": "2026-07-04T22:18:06.085721Z", "iopub.status.idle": "2026-07-04T22:18:10.678887Z", "shell.execute_reply": "2026-07-04T22:18:10.677820Z" } }, "outputs": [], "source": [ "from plots import plot_stgp_spatial_programs\n", "\n", "scores_df = pd.DataFrame(\n", " adata.obsm[\"X_stgp_spatial\"],\n", " index=adata.obs_names,\n", " columns=[f\"stGP{j+1}\" for j in range(p_sel)],\n", ")\n", "\n", "figs = plot_stgp_spatial_programs(\n", " stgp_adata=adata, scores=scores_df,\n", " celltype=\"Excitatory Neuron\", age_unit=\"years\",\n", " ncols=4, fg_dot_size=5.0, dpi=300,\n", ")\n", "for j, fig in enumerate(figs):\n", " fig.savefig(FIGURES_DIR / f\"spatial_stGP{j+1}.png\", dpi=300, bbox_inches=\"tight\")\n", " plt.close(fig)" ] }, { "cell_type": "code", "execution_count": 12, "id": "a1581d3a", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:10.682062Z", "iopub.status.busy": "2026-07-04T22:18:10.681891Z", "iopub.status.idle": "2026-07-04T22:18:10.686398Z", "shell.execute_reply": "2026-07-04T22:18:10.685795Z" } }, "outputs": [ { "data": { "text/plain": [ "
" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "figs[0]" ] }, { "cell_type": "markdown", "id": "md-sec4-3", "metadata": {}, "source": [ "### 3.3 Age Trajectories (α)\n", "\n", "**α(t)** is the posterior mean age effect of each program — it quantifies how the\n", "program amplitude changes across the human lifespan (15–87 yr).\n", "The shaded band shows the 95% posterior credible interval." ] }, { "cell_type": "code", "execution_count": 13, "id": "age-trajectories", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:10.688847Z", "iopub.status.busy": "2026-07-04T22:18:10.688576Z", "iopub.status.idle": "2026-07-04T22:18:11.609321Z", "shell.execute_reply": "2026-07-04T22:18:11.608541Z" } }, "outputs": [], "source": [ "stgp_info = adata.uns['stgp']\n", "ages_slices = np.array(stgp_info['ages'])\n", "alpha = np.array(stgp_info['alpha']) # (p, n_slices)\n", "alpha_lower = np.array(stgp_info['alpha_lower'])\n", "alpha_upper = np.array(stgp_info['alpha_upper'])\n", "\n", "COLOR = '#2C7FB8'\n", "order = np.argsort(ages_slices)\n", "t = ages_slices[order]\n", "\n", "for j in range(p_sel):\n", " a = alpha[j][order]\n", " lo = alpha_lower[j][order]\n", " hi = alpha_upper[j][order]\n", "\n", " fig, ax = plt.subplots(figsize=(4.5, 4), constrained_layout=True)\n", "\n", " ax.fill_between(t, lo, hi, alpha=0.18, color=COLOR)\n", " ax.plot(t, lo, lw=0.8, ls='--', color=COLOR, alpha=0.55)\n", " ax.plot(t, hi, lw=0.8, ls='--', color=COLOR, alpha=0.55)\n", " ax.plot(t, a, lw=1.8, color=COLOR)\n", " ax.scatter(t, a, s=32, color=COLOR, zorder=3, label='Posterior mean')\n", " ax.axhline(0, color='0.6', lw=0.7, ls=':')\n", "\n", " ax.set_xlabel('Age (yr)')\n", " ax.set_ylabel('Age effect α')\n", " ax.legend(fontsize=9)\n", " fig.savefig(FIGURES_DIR / f\"alpha_trajectory_stGP{j+1}.png\", dpi=400, bbox_inches=\"tight\")\n", " plt.close(fig)" ] }, { "cell_type": "markdown", "id": "md-sec4-4", "metadata": {}, "source": [ "### 3.4 GP Parameters (θ)" ] }, { "cell_type": "code", "execution_count": 14, "id": "2e722771", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:11.612021Z", "iopub.status.busy": "2026-07-04T22:18:11.611852Z", "iopub.status.idle": "2026-07-04T22:18:11.615133Z", "shell.execute_reply": "2026-07-04T22:18:11.614609Z" } }, "outputs": [ { "data": { "text/plain": [ "array([[ 38.81529414, 15.40473634],\n", " [973.40129965, 12.72626063],\n", " [223.0523848 , 11.20525559],\n", " [ 84.11608287, 147.19819587]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "theta = np.asarray(stgp_info[\"theta\"], dtype=float)\n", "theta" ] }, { "cell_type": "code", "execution_count": 15, "id": "theta-bar", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:11.616819Z", "iopub.status.busy": "2026-07-04T22:18:11.616688Z", "iopub.status.idle": "2026-07-04T22:18:12.034473Z", "shell.execute_reply": "2026-07-04T22:18:12.033905Z" } }, "outputs": [], "source": [ "theta = np.array(stgp_info['theta']) # (p, 2): [amplitude, noise_frac]\n", "prog_names = [f'stGP{j+1}' for j in range(p_sel)]\n", "prog_colors = plt.cm.tab10.colors[:p_sel]\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), constrained_layout=True)\n", "\n", "for j, (col, name) in enumerate(zip(prog_colors, prog_names)):\n", " axes[0].bar(j, theta[j, 0], color=col, edgecolor='white', linewidth=0.6)\n", " axes[1].bar(j, theta[j, 1], color=col, edgecolor='white', linewidth=0.6)\n", "\n", "for ax in axes:\n", " ax.set_xticks(range(p_sel))\n", " ax.set_xticklabels(prog_names, rotation=30, ha='right')\n", " ax.set_xlabel('Program')\n", "\n", "axes[0].set_ylabel(r'$\\sigma_{\\mathrm{age}}^2$')\n", "axes[0].set_title(r'Temporal variance component ($\\sigma_{\\mathrm{age}}^2$)')\n", "axes[1].set_ylabel(r'$\\tau_{\\mathrm{spa}}^2$')\n", "axes[1].set_title(r'Spatial variance component ($\\tau_{\\mathrm{spa}}^2$)')\n", "\n", "fig.savefig(FIGURES_DIR / \"variance_components.png\", dpi=400, bbox_inches=\"tight\")" ] }, { "cell_type": "markdown", "id": "md-sec4-6", "metadata": {}, "source": [ "### 3.5 Spatial Visualisation of a Single Slice\n", "\n", "We inspect one tissue slice in detail, showing the **spatial b field** of each program\n", "(the smooth within-slice component)." ] }, { "cell_type": "code", "execution_count": 16, "id": "single-slice", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:12.036704Z", "iopub.status.busy": "2026-07-04T22:18:12.036594Z", "iopub.status.idle": "2026-07-04T22:18:13.008303Z", "shell.execute_reply": "2026-07-04T22:18:13.007040Z" } }, "outputs": [], "source": [ "slice_ages = [(adata.obs.loc[adata.obs['id_region'] == sid, 'age'].iloc[0], sid)\n", " for sid in adata.obs['id_region'].unique()]\n", "slice_ages.sort()\n", "_, example_slice = slice_ages[len(slice_ages) // 2]\n", "\n", "sub = adata[adata.obs['id_region'].astype(str) == example_slice].copy()\n", "age_val = sub.obs['age'].iloc[0]\n", "\n", "fig, axes = plt.subplots(1, p_sel, figsize=(4.5 * p_sel, 4.5), constrained_layout=True)\n", "b = sub.obsm['X_stgp_spatial']\n", "xy = np.asarray(sub.obsm['spatial'])\n", "for j, ax in enumerate(np.atleast_1d(axes)):\n", " v99 = np.nanpercentile(np.abs(b[:, j]), 99)\n", " sc_ref = ax.scatter(xy[:, 0], xy[:, 1], c=b[:, j],\n", " cmap='RdBu_r', vmin=-v99, vmax=v99,\n", " s=10, linewidths=0, rasterized=True)\n", " ax.set_aspect('equal'); ax.axis('off')\n", " ax.set_title(f'stGP{j+1}')\n", " plt.colorbar(sc_ref, ax=ax, shrink=0.7, pad=0.01)\n", "fig.savefig(FIGURES_DIR / f\"spatial_b_{example_slice}.png\", dpi=400, bbox_inches=\"tight\")" ] }, { "cell_type": "markdown", "id": "md-benchmark", "metadata": {}, "source": [ "---\n", "## 4. Benchmarking Analysis\n", "\n", "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.\n", "\n", "This section compares stGP against STAMP, MEFISTO, Popari, and SpatialPCA using three complementary views:\n", "\n", "- marker-gene/program correlations for CUX2, RORB, and HS3ST4;\n", "- spatial embedding panels and representative slice panels;\n", "- clustering recovery against marker-derived layer labels and high-resolution `celltype2` layer labels.\n", "\n", "All outputs are written under `Figure/ext/benchmark`." ] }, { "cell_type": "code", "execution_count": 17, "id": "benchmark-load-inputs", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:13.011207Z", "iopub.status.busy": "2026-07-04T22:18:13.011031Z", "iopub.status.idle": "2026-07-04T22:18:15.915667Z", "shell.execute_reply": "2026-07-04T22:18:15.915108Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
datasetcellsslicesembedding
0stGP6116212X_stgp_spatial
1STAMP6116212X_stamp
2MEFISTO6116212X_mefisto
3Popari6116212X
4SpatialPCA6116212X_spatialpca
\n", "
" ], "text/plain": [ " dataset cells slices embedding\n", "0 stGP 61162 12 X_stgp_spatial\n", "1 STAMP 61162 12 X_stamp\n", "2 MEFISTO 61162 12 X_mefisto\n", "3 Popari 61162 12 X\n", "4 SpatialPCA 61162 12 X_spatialpca" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from benchmarking_ext import BASELINE_OBSM_KEYS, LAYER_MARKERS, LAYER_SAFE, METHODS, run_ext_benchmarking\n", "\n", "OUT_DIR = RESULTS_DIR / CELLTYPE\n", "ADATA_PATH = OUT_DIR / \"adata_with_scores.h5ad\"\n", "\n", "if \"adata\" not in globals() or \"X_stgp_spatial\" not in adata.obsm:\n", " adata = sc.read_h5ad(ADATA_PATH)\n", "\n", "if \"adata_prep\" not in globals():\n", " adata_temp = adata_ext.copy()\n", " sc.pp.normalize_total(adata_temp, target_sum=1e3)\n", " adata_prep = sc.pp.log1p(adata_temp, copy=True)\n", "\n", "if \"slices\" not in globals():\n", " if \"stgp\" in adata.uns and \"groups\" in adata.uns[\"stgp\"]:\n", " slices = np.array(adata.uns[\"stgp\"][\"groups\"])\n", " else:\n", " slice_ids = adata.obs[\"id_region\"].astype(str)\n", " slices = np.array(\n", " sorted(\n", " pd.unique(slice_ids),\n", " key=lambda sid: float(adata.obs.loc[slice_ids == sid, \"age\"].iloc[0]),\n", " )\n", " )\n", "\n", "BL_DIR = Path(\"Results/baselines\")\n", "baseline_adatas = {\n", " \"STAMP\": sc.read_h5ad(BL_DIR / \"stamp_k=3/ext/adata_with_scores.h5ad\"),\n", " \"MEFISTO\": sc.read_h5ad(BL_DIR / \"mefisto/ext/adata_with_scores.h5ad\"),\n", " \"Popari\": sc.read_h5ad(BL_DIR / \"popari/ext/res_popari.h5ad\"),\n", " \"SpatialPCA\": sc.read_h5ad(BL_DIR / \"spatialpca/ext/adata_with_scores.h5ad\"),\n", "}\n", "\n", "for method, obsm_key in BASELINE_OBSM_KEYS.items():\n", " if obsm_key not in baseline_adatas[method].obsm:\n", " raise KeyError(f\"{method} is missing .obsm[{obsm_key!r}]\")\n", " if baseline_adatas[method].n_obs != adata.n_obs:\n", " raise ValueError(f\"{method} has {baseline_adatas[method].n_obs} cells; expected {adata.n_obs}\")\n", "\n", "input_summary = pd.DataFrame(\n", " [\n", " {\"dataset\": \"stGP\", \"cells\": adata.n_obs, \"slices\": len(slices), \"embedding\": \"X_stgp_spatial\"},\n", " *[\n", " {\n", " \"dataset\": method,\n", " \"cells\": baseline_adatas[method].n_obs,\n", " \"slices\": baseline_adatas[method].obs[\"id_region\"].nunique(),\n", " \"embedding\": BASELINE_OBSM_KEYS[method],\n", " }\n", " for method in BASELINE_OBSM_KEYS\n", " ],\n", " ]\n", ")\n", "\n", "display(input_summary)" ] }, { "cell_type": "code", "execution_count": 18, "id": "benchmark-run", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:18:15.917197Z", "iopub.status.busy": "2026-07-04T22:18:15.917075Z", "iopub.status.idle": "2026-07-04T22:24:12.800597Z", "shell.execute_reply": "2026-07-04T22:24:12.799766Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Benchmark outputs saved under: Figure/ext/benchmark\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
n_outputs
correlation_figures4
spatial_figures24
cluster_figures36
metric_figures10
summary_figures13
\n", "
" ], "text/plain": [ " n_outputs\n", "correlation_figures 4\n", "spatial_figures 24\n", "cluster_figures 36\n", "metric_figures 10\n", "summary_figures 13" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "BENCHMARK_DIR = FIGURES_DIR / \"benchmark\"\n", "\n", "benchmark_outputs = run_ext_benchmarking(\n", " adata=adata,\n", " adata_prep=adata_prep,\n", " baseline_adatas=baseline_adatas,\n", " benchmark_dir=BENCHMARK_DIR,\n", " slices=slices,\n", " methods=METHODS,\n", " dpi=400,\n", ")\n", "\n", "output_counts = pd.Series(\n", " {\n", " \"correlation_figures\": len(benchmark_outputs[\"correlation_figures\"]),\n", " \"spatial_figures\": len(benchmark_outputs[\"spatial_figures\"]),\n", " \"cluster_figures\": len(benchmark_outputs[\"cluster_figures\"]),\n", " \"metric_figures\": len(benchmark_outputs[\"metric_figures\"]),\n", " \"summary_figures\": len(benchmark_outputs[\"summary_figures\"]),\n", " },\n", " name=\"n_outputs\",\n", ")\n", "\n", "print(f\"Benchmark outputs saved under: {benchmark_outputs['benchmark_dir']}\")\n", "display(output_counts.to_frame())" ] }, { "cell_type": "code", "execution_count": 19, "id": "benchmark-review-source-data", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:24:12.802927Z", "iopub.status.busy": "2026-07-04T22:24:12.802794Z", "iopub.status.idle": "2026-07-04T22:24:12.822334Z", "shell.execute_reply": "2026-07-04T22:24:12.821823Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Unnamed: 0layermarker_genemethodmeanmedianstd
00L2/3CUX2MEFISTO0.7493110.8158600.207515
11L2/3CUX2Popari0.3629300.2485150.239036
22L2/3CUX2STAMP0.6103710.7375650.281830
33L2/3CUX2SpatialPCA0.5459550.4639470.222554
44L2/3CUX2stGP0.7968510.8385560.155620
55L4RORBMEFISTO0.6043040.6128080.096769
66L4RORBPopari0.2832650.2420780.265295
77L4RORBSTAMP0.3716380.3426780.179853
88L4RORBSpatialPCA0.2579410.3014750.158959
99L4RORBstGP0.6746250.6995950.097605
1010L5/6HS3ST4MEFISTO0.4471950.4681780.102479
1111L5/6HS3ST4Popari0.1258160.0564750.137783
1212L5/6HS3ST4STAMP0.2270080.3343050.238527
1313L5/6HS3ST4SpatialPCA0.3252660.3231950.140685
1414L5/6HS3ST4stGP0.5153560.5515360.101769
\n", "
" ], "text/plain": [ " Unnamed: 0 layer marker_gene method mean median std\n", "0 0 L2/3 CUX2 MEFISTO 0.749311 0.815860 0.207515\n", "1 1 L2/3 CUX2 Popari 0.362930 0.248515 0.239036\n", "2 2 L2/3 CUX2 STAMP 0.610371 0.737565 0.281830\n", "3 3 L2/3 CUX2 SpatialPCA 0.545955 0.463947 0.222554\n", "4 4 L2/3 CUX2 stGP 0.796851 0.838556 0.155620\n", "5 5 L4 RORB MEFISTO 0.604304 0.612808 0.096769\n", "6 6 L4 RORB Popari 0.283265 0.242078 0.265295\n", "7 7 L4 RORB STAMP 0.371638 0.342678 0.179853\n", "8 8 L4 RORB SpatialPCA 0.257941 0.301475 0.158959\n", "9 9 L4 RORB stGP 0.674625 0.699595 0.097605\n", "10 10 L5/6 HS3ST4 MEFISTO 0.447195 0.468178 0.102479\n", "11 11 L5/6 HS3ST4 Popari 0.125816 0.056475 0.137783\n", "12 12 L5/6 HS3ST4 STAMP 0.227008 0.334305 0.238527\n", "13 13 L5/6 HS3ST4 SpatialPCA 0.325266 0.323195 0.140685\n", "14 14 L5/6 HS3ST4 stGP 0.515356 0.551536 0.101769" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
methodn_slicesraw_ariraw_nmiraw_acc
4stGP120.3227100.3118160.687081
0MEFISTO120.2728900.2768350.661968
2STAMP120.2511360.2257010.626740
3SpatialPCA120.1488980.1694660.543434
1Popari120.1078900.1090650.502050
\n", "
" ], "text/plain": [ " method n_slices raw_ari raw_nmi raw_acc\n", "4 stGP 12 0.322710 0.311816 0.687081\n", "0 MEFISTO 12 0.272890 0.276835 0.661968\n", "2 STAMP 12 0.251136 0.225701 0.626740\n", "3 SpatialPCA 12 0.148898 0.169466 0.543434\n", "1 Popari 12 0.107890 0.109065 0.502050" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
methodn_slicesraw_ariraw_nmiraw_acc
4stGP90.5880800.5501210.829895
0MEFISTO90.5656880.5315910.820396
2STAMP90.4070570.3897120.729004
3SpatialPCA90.2339730.2628110.606660
1Popari90.1798290.1646450.540688
\n", "
" ], "text/plain": [ " method n_slices raw_ari raw_nmi raw_acc\n", "4 stGP 9 0.588080 0.550121 0.829895\n", "0 MEFISTO 9 0.565688 0.531591 0.820396\n", "2 STAMP 9 0.407057 0.389712 0.729004\n", "3 SpatialPCA 9 0.233973 0.262811 0.606660\n", "1 Popari 9 0.179829 0.164645 0.540688" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "source_dir = Path(benchmark_outputs[\"source_dir\"]) / \"summary\" / \"source_data\"\n", "\n", "corr_summary = pd.read_csv(source_dir / \"marker_embedding_correlation_summary.csv\")\n", "display(corr_summary)\n", "\n", "for ground_truth, metrics_df in benchmark_outputs[\"cluster_metrics\"].items():\n", " display(\n", " metrics_df\n", " .groupby(\"method\", as_index=False)\n", " .agg(\n", " n_slices=(\"id_region\", \"nunique\"),\n", " raw_ari=(\"raw_ari\", \"mean\"),\n", " raw_nmi=(\"raw_nmi\", \"mean\"),\n", " raw_acc=(\"raw_acc\", \"mean\"),\n", " )\n", " .sort_values(\"raw_ari\", ascending=False)\n", " )" ] }, { "cell_type": "markdown", "id": "md-sec4-7", "metadata": {}, "source": [ "---\n", "## 5. Pathway Enrichment Analysis\n", "\n", "We use **gseapy** to run over-representation analysis (ORA) on the\n", "positive-weight genes of each stGP program, testing against GO Biological\n", "Process and GO Cellular Component gene sets. For excitatory neurons, we\n", "expect programs enriched in synaptic transmission, neuronal differentiation,\n", "axon guidance, and age-related processes such as protein folding stress and\n", "mitochondrial dysfunction.\n", "\n", "> **Prerequisites:** Download MSigDB gene-set files from https://www.gsea-msigdb.org/\n", "> and place them under `data/genesets/`:\n", "> - `c5.go.bp.v2026.1.Hs.symbols.gmt`\n", "> - `c5.go.cc.v2026.1.Hs.symbols.gmt`\n" ] }, { "cell_type": "code", "execution_count": 20, "id": "enrichment-run", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:24:12.823981Z", "iopub.status.busy": "2026-07-04T22:24:12.823878Z", "iopub.status.idle": "2026-07-04T22:24:14.962671Z", "shell.execute_reply": "2026-07-04T22:24:14.961699Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "stGP1: 15 active genes\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "stGP2: 15 active genes\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "stGP3: 15 active genes\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "stGP4: 15 active genes\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "All enrichment results saved to Figure/ext/enrichment_results.csv\n" ] } ], "source": [ "from plots import plot_enrichment_dotplot, truncate_colormap as _trunc_cmap\n", "import gseapy as gp\n", "\n", "gene_sets = {\n", " \"GO Biological Process\": \"data/genesets/c5.go.bp.v2026.1.Hs.symbols.gmt\",\n", " \"GO Cellular Component\": \"data/genesets/c5.go.cc.v2026.1.Hs.symbols.gmt\",\n", "}\n", "cmaps = {\n", " \"GO Biological Process\": _trunc_cmap(\"Reds\", 0.22, 0.78),\n", " \"GO Cellular Component\": _trunc_cmap(\"Purples\", 0.25, 0.76),\n", "}\n", "background_genes = list(W_df.columns)\n", "\n", "all_res = []\n", "for program in W_df.index:\n", " gene_list = W_df.loc[program].pipe(lambda s: s[s > 0].sort_values(ascending=False).index.tolist())\n", " print(f\"{program}: {len(gene_list)} active genes\")\n", "\n", " fig, axes = plt.subplots(len(gene_sets), 1,\n", " figsize=(7.5, 2.8 * len(gene_sets)), constrained_layout=True)\n", " for ax, (set_name, gmt_file) in zip(np.atleast_1d(axes), gene_sets.items()):\n", " enr = gp.enrich(gene_list=gene_list, gene_sets=gmt_file,\n", " background=background_genes, verbose=False)\n", " res = enr.res2d.copy()\n", " res[\"program\"] = program\n", " all_res.append(res)\n", " plot_enrichment_dotplot(res, ax, set_name, cmap=cmaps[set_name])\n", "\n", " fig.savefig(FIGURES_DIR / f\"{program}_enrichment.png\", dpi=300, bbox_inches=\"tight\")\n", " plt.show()\n", "\n", "enr_df = pd.concat(all_res, ignore_index=True)\n", "enr_df.to_csv(FIGURES_DIR / \"enrichment_results.csv\", index=False)\n", "print(f\"\\nAll enrichment results saved to {FIGURES_DIR}/enrichment_results.csv\")" ] }, { "cell_type": "code", "execution_count": 21, "id": "ac84927d", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T22:24:14.964590Z", "iopub.status.busy": "2026-07-04T22:24:14.964472Z", "iopub.status.idle": "2026-07-04T22:24:14.970076Z", "shell.execute_reply": "2026-07-04T22:24:14.969472Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Gene_setTermOverlapP-valueAdjusted P-valueOdds RatioCombined ScoreGenesprogram
41c5.go.cc.v2026.1.Hs.symbols.gmtGOCC_GLUTAMATERGIC_SYNAPSE6/170.0000640.00687015.439359149.042913C1QL3;SYN3;GRIK4;APOE;CBLN2;EPHB1stGP1
91c5.go.cc.v2026.1.Hs.symbols.gmtGOCC_SYNAPSE8/450.0005830.0208057.05688952.551218C1QL3;SV2C;LAMP5;SYN3;GRIK4;APOE;CBLN2;EPHB1stGP1
92c5.go.cc.v2026.1.Hs.symbols.gmtGOCC_SYNAPTIC_CLEFT3/40.0004620.02080550.306667386.405222C1QL3;CBLN2;APOEstGP1
\n", "
" ], "text/plain": [ " Gene_set Term Overlap \\\n", "41 c5.go.cc.v2026.1.Hs.symbols.gmt GOCC_GLUTAMATERGIC_SYNAPSE 6/17 \n", "91 c5.go.cc.v2026.1.Hs.symbols.gmt GOCC_SYNAPSE 8/45 \n", "92 c5.go.cc.v2026.1.Hs.symbols.gmt GOCC_SYNAPTIC_CLEFT 3/4 \n", "\n", " P-value Adjusted P-value Odds Ratio Combined Score \\\n", "41 0.000064 0.006870 15.439359 149.042913 \n", "91 0.000583 0.020805 7.056889 52.551218 \n", "92 0.000462 0.020805 50.306667 386.405222 \n", "\n", " Genes program \n", "41 C1QL3;SYN3;GRIK4;APOE;CBLN2;EPHB1 stGP1 \n", "91 C1QL3;SV2C;LAMP5;SYN3;GRIK4;APOE;CBLN2;EPHB1 stGP1 \n", "92 C1QL3;CBLN2;APOE stGP1 " ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "all_res[1][all_res[1]['Adjusted P-value']<0.05]" ] }, { "cell_type": "code", "execution_count": null, "id": "e6ba5e2b", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "stGP", "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.11.15" } }, "nbformat": 4, "nbformat_minor": 5 }