{ "cells": [ { "cell_type": "markdown", "id": "title", "metadata": {}, "source": [ "# Injured PT stGP analysis: right kidney\n", "\n", "Fit stGP on the right biological replicate (`ident` suffix `R`) of injured proximal tubule cells. Run this notebook with the `stGP` conda environment." ] }, { "cell_type": "code", "execution_count": 1, "id": "setup", "metadata": { "execution": { "iopub.execute_input": "2026-05-31T09:36:49.837267Z", "iopub.status.busy": "2026-05-31T09:36:49.837152Z", "iopub.status.idle": "2026-05-31T09:37:11.801437Z", "shell.execute_reply": "2026-05-31T09:37:11.799667Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AnnData object with n_obs × n_vars = 76699 × 299\n", " obs: 'x_centroid', 'y_centroid', 'n_genes', 'n_counts', 'ident', 'region', 'celltype_plot', 'time', 'CN', 'injury_time_days', 'side', 'age'\n", " uns: 'CN_colors', 'celltype_plot_colors', 'ident_colors', 'neighbors', 'pca', 'umap'\n", " obsm: 'X_pca', 'X_pca_harmony', 'X_umap', 'spatial'\n", " varm: 'PCs'\n", " obsp: 'connectivities', 'distances'\n", " ident time injury_time_days side\n", " Hour4R Hour4 0.166667 R\n", "Hour12R Hour12 0.500000 R\n", " Day2R Day2 2.000000 R\n", " Day14R Day14 14.000000 R\n", " Week6R Week6 42.000000 R\n" ] } ], "source": [ "import os, 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", "\n", "warnings.filterwarnings(\"ignore\", category=FutureWarning)\n", "sys.path.insert(0, \"..\")\n", "\n", "DATA_PROC = Path(\"data/processed\")\n", "RESULTS_DIR = Path(\"Results/stgp\")\n", "FIGURES_DIR = Path(\"Figures/Inj_PT_R\")\n", "FIGURES_DIR.mkdir(parents=True, exist_ok=True)\n", "\n", "CELLTYPE = \"Inj_PT\"\n", "SIDE = \"R\"\n", "ANALYSIS_NAME = f\"{CELLTYPE}_{SIDE}\"\n", "\n", "adata_inj_pt = sc.read_h5ad(DATA_PROC / f\"{CELLTYPE}.h5ad\")\n", "adata_inj_pt.obs[\"side\"] = adata_inj_pt.obs[\"ident\"].astype(str).str[-1]\n", "adata_inj_pt = adata_inj_pt[adata_inj_pt.obs[\"injury_time_days\"] > 0].copy()\n", "adata_inj_pt = adata_inj_pt[adata_inj_pt.obs[\"side\"] == SIDE].copy()\n", "adata_inj_pt.obs[\"age\"] = adata_inj_pt.obs[\"injury_time_days\"].copy()\n", "\n", "print(adata_inj_pt)\n", "print(adata_inj_pt.obs[[\"ident\", \"time\", \"injury_time_days\", \"side\"]].drop_duplicates().sort_values(\"injury_time_days\").to_string(index=False))" ] }, { "cell_type": "code", "execution_count": 2, "id": "global-style", "metadata": { "execution": { "iopub.execute_input": "2026-05-31T09:37:11.806225Z", "iopub.status.busy": "2026-05-31T09:37:11.805408Z", "iopub.status.idle": "2026-05-31T09:37:11.812211Z", "shell.execute_reply": "2026-05-31T09:37:11.811512Z" } }, "outputs": [], "source": [ "import matplotlib as mpl\n", "\n", "mpl.rcParams.update({\n", " \"font.family\": \"sans-serif\",\n", " \"font.sans-serif\": [\"Arial\", \"Helvetica\", \"DejaVu Sans\"],\n", " \"font.size\": 11,\n", " \"pdf.fonttype\": 42,\n", " \"ps.fonttype\": 42,\n", " \"svg.fonttype\": \"none\",\n", " \"figure.dpi\": 150,\n", " \"savefig.dpi\": 400,\n", " \"figure.facecolor\": \"white\",\n", " \"axes.facecolor\": \"white\",\n", " \"axes.linewidth\": 1.2,\n", " \"axes.spines.top\": False,\n", " \"axes.spines.right\": False,\n", " \"axes.labelsize\": 11,\n", " \"axes.titlesize\": 12,\n", " \"xtick.direction\": \"out\",\n", " \"ytick.direction\": \"out\",\n", " \"legend.frameon\": False,\n", " \"legend.fontsize\": 9,\n", " \"lines.linewidth\": 1.5,\n", "})" ] }, { "cell_type": "code", "execution_count": 3, "id": "prepare-stgp", "metadata": { "execution": { "iopub.execute_input": "2026-05-31T09:37:11.815079Z", "iopub.status.busy": "2026-05-31T09:37:11.814944Z", "iopub.status.idle": "2026-05-31T09:37:15.671043Z", "shell.execute_reply": "2026-05-31T09:37:15.670026Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/byual/.conda/envs/stGP/lib/python3.11/functools.py:909: UserWarning: zero-centering a sparse array/matrix densifies it.\n", " return dispatch(args[0].__class__)(*args, **kw)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Prepared Inj_PT_R: 76699 cells across 5 slices\n", " slice age n_cells\n", " Hour4R 0.166667 18287\n", "Hour12R 0.500000 38001\n", " Day2R 2.000000 19408\n", " Day14R 14.000000 820\n", " Week6R 42.000000 183\n" ] } ], "source": [ "import time\n", "from stgp.estimation import fit_pfactor_auto\n", "from stgp.kernels import (\n", " bandwidth_select_spatial,\n", " bandwidth_select_temporal,\n", " build_K_age,\n", " build_K_spa_list_from_stacked,\n", ")\n", "from stgp.preprocessing import standardize_coords_list\n", "\n", "OUT_DIR = RESULTS_DIR / ANALYSIS_NAME\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_inj_pt.obs[\"injury_time_days\"], errors=\"coerce\").to_numpy(float)\n", "groups = adata_inj_pt.obs[\"ident\"].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_prep = adata_inj_pt.copy()\n", "sc.pp.scale(adata_prep)\n", "Y_list = [adata_prep.X[ix] for ix in idx_per_group]\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)\n", "ages = ages[sort_ord]\n", "slices = uniq[sort_ord]\n", "idx_sorted = [idx_per_group[i] for i in sort_ord]\n", "nlist = nlist[sort_ord]\n", "Y_list = [Y_list[i] for i in sort_ord]\n", "\n", "coords_list = standardize_coords_list([adata_inj_pt.obsm[\"spatial\"][ix] for ix in idx_per_group])\n", "coords_list = [coords_list[i] for i in sort_ord]\n", "\n", "print(f\"Prepared {ANALYSIS_NAME}: {adata_inj_pt.n_obs} cells across {len(slices)} slices\")\n", "print(pd.DataFrame({\"slice\": slices, \"age\": ages, \"n_cells\": nlist}).to_string(index=False))" ] }, { "cell_type": "code", "execution_count": 4, "id": "build-kernels", "metadata": { "execution": { "iopub.execute_input": "2026-05-31T09:37:15.674269Z", "iopub.status.busy": "2026-05-31T09:37:15.674112Z", "iopub.status.idle": "2026-05-31T09:37:35.772810Z", "shell.execute_reply": "2026-05-31T09:37:35.770652Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "gamma_spa = 0.4106 | gamma_age = 0.3904" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "gamma_spa = bandwidth_select_spatial(coords_list, frac=0.01, rho=0.5)\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=\"ar1\", 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": 5, "id": "fit-stgp", "metadata": { "execution": { "iopub.execute_input": "2026-05-31T09:37:35.777030Z", "iopub.status.busy": "2026-05-31T09:37:35.776785Z", "iopub.status.idle": "2026-05-31T12:04:06.868900Z", "shell.execute_reply": "2026-05-31T12:04:06.866940Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[sweep=001] dW_rel=1.656e-01 dTheta_rel=1.727e-01 time=5.957e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=002] dW_rel=1.568e-01 dTheta_rel=2.711e-02 time=4.321e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=003] dW_rel=8.584e-02 dTheta_rel=2.277e-02 time=4.079e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=004] dW_rel=6.787e-02 dTheta_rel=1.343e-02 time=3.357e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=005] dW_rel=7.023e-02 dTheta_rel=1.177e-02 time=3.423e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=006] dW_rel=1.876e-02 dTheta_rel=8.113e-03 time=2.782e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=007] dW_rel=6.750e-02 dTheta_rel=9.243e-03 time=2.930e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=008] dW_rel=1.176e-01 dTheta_rel=2.817e-02 time=4.225e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=009] dW_rel=1.425e-01 dTheta_rel=2.328e-02 time=4.059e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=010] dW_rel=2.780e-02 dTheta_rel=1.921e-02 time=4.274e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=011] dW_rel=7.550e-02 dTheta_rel=2.859e-02 time=3.203e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=012] dW_rel=1.575e-02 dTheta_rel=3.185e-03 time=1.747e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=013] dW_rel=7.294e-03 dTheta_rel=4.125e-03 time=1.259e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=014] dW_rel=4.750e-03 dTheta_rel=3.447e-03 time=1.283e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=015] dW_rel=3.394e-03 dTheta_rel=2.579e-03 time=9.984e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=016] dW_rel=2.509e-03 dTheta_rel=1.905e-03 time=1.208e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=017] dW_rel=1.885e-03 dTheta_rel=1.430e-03 time=1.382e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=018] dW_rel=1.429e-03 dTheta_rel=1.088e-03 time=1.590e+02\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=019] dW_rel=1.088e-03 dTheta_rel=8.401e-04 time=7.327e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=020] dW_rel=8.303e-04 dTheta_rel=6.476e-04 time=7.566e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=021] dW_rel=6.372e-04 dTheta_rel=5.003e-04 time=2.745e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=022] dW_rel=4.894e-04 dTheta_rel=3.863e-04 time=2.772e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=023] dW_rel=3.768e-04 dTheta_rel=2.983e-04 time=4.946e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=024] dW_rel=2.905e-04 dTheta_rel=2.308e-04 time=5.718e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=025] dW_rel=2.237e-04 dTheta_rel=1.813e-04 time=3.391e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=026] dW_rel=1.735e-04 dTheta_rel=1.388e-04 time=1.904e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=027] dW_rel=1.335e-04 dTheta_rel=1.068e-04 time=1.577e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=028] dW_rel=1.033e-04 dTheta_rel=8.551e-05 time=1.583e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[auto_rank prune] p=8 -> 6 (dropped factors [6, 7])\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=001] dW_rel=2.848e-01 dTheta_rel=3.871e-02 time=5.849e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=002] dW_rel=1.247e-01 dTheta_rel=3.304e-02 time=6.729e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=003] dW_rel=4.483e-02 dTheta_rel=2.465e-02 time=4.213e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=004] dW_rel=1.371e-02 dTheta_rel=9.353e-03 time=3.109e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=005] dW_rel=7.102e-03 dTheta_rel=4.718e-03 time=2.319e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=006] dW_rel=4.290e-03 dTheta_rel=2.881e-03 time=2.372e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=007] dW_rel=2.964e-03 dTheta_rel=1.976e-03 time=1.926e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=008] dW_rel=2.226e-03 dTheta_rel=1.447e-03 time=2.113e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=009] dW_rel=1.741e-03 dTheta_rel=1.096e-03 time=1.925e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=010] dW_rel=1.383e-03 dTheta_rel=8.476e-04 time=1.903e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=011] dW_rel=1.107e-03 dTheta_rel=6.660e-04 time=1.628e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=012] dW_rel=8.859e-04 dTheta_rel=5.235e-04 time=1.572e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=013] dW_rel=7.123e-04 dTheta_rel=4.176e-04 time=1.349e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=014] dW_rel=5.735e-04 dTheta_rel=3.310e-04 time=1.370e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=015] dW_rel=4.614e-04 dTheta_rel=2.665e-04 time=1.359e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=016] dW_rel=3.713e-04 dTheta_rel=2.142e-04 time=1.257e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=017] dW_rel=2.994e-04 dTheta_rel=1.692e-04 time=1.174e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=018] dW_rel=2.414e-04 dTheta_rel=1.377e-04 time=1.187e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=019] dW_rel=1.947e-04 dTheta_rel=1.112e-04 time=1.086e+01\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=020] dW_rel=1.569e-04 dTheta_rel=8.966e-05 time=8.831e+00\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=021] dW_rel=1.266e-04 dTheta_rel=7.218e-05 time=8.632e+00\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[sweep=022] dW_rel=1.022e-04 dTheta_rel=5.805e-05 time=8.557e+00\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Runtime: 8790.9s | programs selected: 6\n", "Saved: Results/stgp/Inj_PT_R/stgp_result.pkl\n" ] } ], "source": [ "t0 = time.perf_counter()\n", "res = fit_pfactor_auto(\n", " Y_list=Y_list,\n", " Nlist=nlist,\n", " K_age=K_age,\n", " Kspa_list=K_spa_list,\n", " p_max=10,\n", " k=15,\n", " inner_rank1_tol=1e-4,\n", " rel_improve_total_tol=0.002,\n", " backfit_tol=1e-4,\n", " prune_energy_frac=0.005,\n", " random_state=0,\n", " verbose=1,\n", ")\n", "print(f\"Runtime: {time.perf_counter() - t0:.1f}s | programs selected: {res['W'].shape[0]}\")\n", "\n", "res[\"gamma_age\"] = gamma_age\n", "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": 6, "id": "attach-scores", "metadata": { "execution": { "iopub.execute_input": "2026-05-31T12:04:06.873573Z", "iopub.status.busy": "2026-05-31T12:04:06.873220Z", "iopub.status.idle": "2026-05-31T12:05:17.304718Z", "shell.execute_reply": "2026-05-31T12:05:17.302777Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saved: Results/stgp/Inj_PT_R/adata_with_scores.h5ad\n", "Saved: Results/stgp/Inj_PT_R/W.csv\n" ] } ], "source": [ "ADATA_PATH = OUT_DIR / \"adata_with_scores.h5ad\"\n", "\n", "adata = adata_inj_pt.copy()\n", "all_idx = np.concatenate(idx_sorted)\n", "H_arr = np.empty_like(res[\"H\"])\n", "b_arr = np.empty_like(res[\"b\"])\n", "H_arr[all_idx] = res[\"H\"]\n", "b_arr[all_idx] = res[\"b\"]\n", "\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.tolist(),\n", " ages=ages.tolist(),\n", " gamma_age=float(res[\"gamma_age\"]),\n", " 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", "p_sel = res[\"W\"].shape[0]\n", "W_df = pd.DataFrame(\n", " res[\"W\"],\n", " index=[f\"stGP{j + 1}\" for j in range(p_sel)],\n", " columns=adata.var_names.astype(str),\n", ")\n", "W_df.to_csv(OUT_DIR / \"W.csv\")\n", "print(f\"Saved: {OUT_DIR / 'W.csv'}\")" ] } ], "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 }