{ "cells": [ { "cell_type": "markdown", "id": "title", "metadata": {}, "source": [ "# Downstream clustering and trajectory inference for stGP Microglia\n", "\n", "This notebook uses the saved stGP outputs for the mouse brain MERFISH microglia analysis. It focuses on one `mouse_id` at a time, clusters cells in the stGP spatial embedding with a symmetric KNN graph and spectral clustering." ] }, { "cell_type": "markdown", "id": "setup-heading", "metadata": {}, "source": [ "## 1. Setup and load stGP outputs" ] }, { "cell_type": "code", "execution_count": 1, "id": "setup-load", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:06:12.394737Z", "iopub.status.busy": "2026-07-04T23:06:12.394604Z", "iopub.status.idle": "2026-07-04T23:06:20.523499Z", "shell.execute_reply": "2026-07-04T23:06:20.522977Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/byual/anaconda3/envs/NicheScope/lib/python3.9/site-packages/pyslingshot/slingshot.py:11: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n", " from tqdm.autonotebook import tqdm\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Reading NicheScope compatibility copy: /import/home4/byual/stGP-0529/RealData_MouseBrainMERFISH/Results/stgp/Microglia/downstream_cluster/adata_with_scores_niche_scope_compat.h5ad\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "AnnData object with n_obs × n_vars = 3626 × 220\n", " obs: 'volume', 'center_x', 'center_y', 'min_x', 'min_y', 'max_x', 'max_y', 'transcript_count', 'num_detected_genes', 'barcodeCount', 'mouse_id', 'slide_id', 'cohort', 'age', 'batch', 'celltype', 'region', 'subregion', 'bbox_area', 'n_genes_by_counts', 'total_counts'\n", " var: 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'\n", " uns: 'neighbors', 'pca', 'preprocess_info', 'stgp', 'umap'\n", " obsm: 'X_pca', 'X_stgp', 'X_stgp_spatial', 'X_umap', 'spatial'\n", " varm: 'PCs'\n", " obsp: 'connectivities', 'distances'\n", "Analysis selection: mouse_97; n_obs=3626, n_vars=220\n", "mouse_id slide_id age batch\n", " 97 B17 32.6 B\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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
mouse_idageslide_idn_cells
37799732.6B173626
35589328.5B53113
28977518.8B42870
33358624.6B22826
17634230.9A92796
6641412.9A102789
31078019.8B42738
2229574.3B172711
1213.8A82543
15413826.7A72543
11103021.4A62538
399610134.5B32518
2448616.6B52498
23375.4A92497
19824633.2A82434
13243323.5A102398
451119.8A72390
8701915.5A62258
2006533.4B32195
26857015.8B22136
\n", "
" ], "text/plain": [ " mouse_id age slide_id n_cells\n", "3779 97 32.6 B17 3626\n", "3558 93 28.5 B5 3113\n", "2897 75 18.8 B4 2870\n", "3335 86 24.6 B2 2826\n", "1763 42 30.9 A9 2796\n", "664 14 12.9 A10 2789\n", "3107 80 19.8 B4 2738\n", "2229 57 4.3 B17 2711\n", "12 1 3.8 A8 2543\n", "1541 38 26.7 A7 2543\n", "1110 30 21.4 A6 2538\n", "3996 101 34.5 B3 2518\n", "2448 61 6.6 B5 2498\n", "233 7 5.4 A9 2497\n", "1982 46 33.2 A8 2434\n", "1324 33 23.5 A10 2398\n", "451 11 9.8 A7 2390\n", "870 19 15.5 A6 2258\n", "2006 53 3.4 B3 2195\n", "2685 70 15.8 B2 2136" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import shutil, warnings\n", "from pathlib import Path\n", "\n", "import h5py\n", "import numpy as np\n", "import pandas as pd\n", "import scanpy as sc\n", "import scipy.sparse as sp\n", "import matplotlib.pyplot as plt\n", "from sklearn.neighbors import NearestNeighbors\n", "from sklearn.cluster import SpectralClustering\n", "from pyslingshot import Slingshot\n", "\n", "warnings.filterwarnings('ignore', category=FutureWarning)\n", "plt.rcParams['figure.dpi'] = 120\n", "\n", "ROOT = Path.cwd()\n", "CELLTYPE = 'Microglia'\n", "STGP_DIR = ROOT / 'Results' / 'stgp' / CELLTYPE\n", "ADATA_PATH = STGP_DIR / 'adata_with_scores.h5ad'\n", "OUT_DIR = STGP_DIR / 'downstream_cluster'\n", "OUT_DIR.mkdir(parents=True, exist_ok=True)\n", "\n", "# Spatial coordinates are section-local, so run downstream trajectory on one mouse/section.\n", "# Set MOUSE_ID = None to use all microglia, but spatial overlays across mice are not recommended.\n", "MOUSE_ID = '97'\n", "TRUTH_KEY = 'subregion'\n", "ROOT_TRUTH_LABEL = 'CC/ACO'\n", "START_CLUSTER = None # None selects the stGP domain most enriched for ROOT_TRUTH_LABEL.\n", "N_EPOCHS = 10\n", "RANDOM_STATE = 1234\n", "\n", "def _prepare_nichescope_h5ad(path: Path, cache_dir: Path) -> Path:\n", " \"\"\"Return a NicheScope-readable copy when the stGP h5ad contains null encodings.\"\"\"\n", " fixed = cache_dir / f'{path.stem}_niche_scope_compat.h5ad'\n", " if fixed.exists() and fixed.stat().st_mtime >= path.stat().st_mtime:\n", " return fixed\n", "\n", " shutil.copy2(path, fixed)\n", " with h5py.File(fixed, 'a') as h5:\n", " # NicheScope's anndata cannot read the null-encoded Scanpy log1p base.\n", " if 'uns/log1p/base' in h5:\n", " del h5['uns/log1p/base']\n", " return fixed\n", "\n", "def read_h5ad_compat(path: Path, cache_dir: Path):\n", " \"\"\"Read saved stGP AnnData from the NicheScope environment.\"\"\"\n", " fixed = _prepare_nichescope_h5ad(path, cache_dir)\n", " print(f'Reading NicheScope compatibility copy: {fixed}')\n", " return sc.read_h5ad(fixed)\n", "\n", "if not ADATA_PATH.exists():\n", " raise FileNotFoundError(f'Missing saved stGP AnnData: {ADATA_PATH}')\n", "\n", "adata_full = read_h5ad_compat(ADATA_PATH, OUT_DIR)\n", "if 'X_stgp_spatial' not in adata_full.obsm:\n", " raise KeyError('Expected adata.obsm[\"X_stgp_spatial\"] from the stGP Microglia output.')\n", "if TRUTH_KEY not in adata_full.obs:\n", " raise KeyError(f'Expected adata.obs[{TRUTH_KEY!r}] for clustering comparison.')\n", "\n", "if MOUSE_ID is None:\n", " adata = adata_full.copy()\n", " label_suffix = 'all_microglia'\n", "else:\n", " mask = adata_full.obs['mouse_id'].astype(str).to_numpy() == str(MOUSE_ID)\n", " if not mask.any():\n", " available = sorted(adata_full.obs['mouse_id'].astype(str).unique())\n", " raise ValueError(f'MOUSE_ID={MOUSE_ID!r} not found. Available mouse_id values: {available}')\n", " adata = adata_full[mask].copy()\n", " label_suffix = f'mouse_{MOUSE_ID}'\n", "\n", "print(adata)\n", "print(f'Analysis selection: {label_suffix}; n_obs={adata.n_obs}, n_vars={adata.n_vars}')\n", "print(adata.obs[['mouse_id', 'slide_id', 'age', 'batch']].drop_duplicates().to_string(index=False))\n", "display(adata_full.obs.groupby(['mouse_id', 'age', 'slide_id']).size().reset_index(name='n_cells').sort_values('n_cells', ascending=False).head(20))" ] }, { "cell_type": "markdown", "id": "cluster-heading", "metadata": {}, "source": [ "## 2. KNN spectral clustering in the stGP spatial embedding" ] }, { "cell_type": "code", "execution_count": 2, "id": "knn-spectral-cluster", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:06:20.525192Z", "iopub.status.busy": "2026-07-04T23:06:20.524824Z", "iopub.status.idle": "2026-07-04T23:06:22.109979Z", "shell.execute_reply": "2026-07-04T23:06:22.109131Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "truth=subregion, n_clusters=4, k_nn=60\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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
subregionCC/ACOCTX_L1/MENCTX_L2/3CTX_L4/5/6STR_CP/ACBSTR_LS/NDBVEN
stGP
1192951471935193311
25133657184242790
3319567196102
42802011910
\n", "
" ], "text/plain": [ "subregion CC/ACO CTX_L1/MEN CTX_L2/3 CTX_L4/5/6 STR_CP/ACB STR_LS/NDB \\\n", "stGP \n", "1 192 95 147 193 519 331 \n", "2 5 133 657 184 242 79 \n", "3 319 5 6 7 196 10 \n", "4 280 2 0 1 19 1 \n", "\n", "subregion VEN \n", "stGP \n", "1 1 \n", "2 0 \n", "3 2 \n", "4 0 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Microglia has a single celltype label, so subregion is used as the comparison label.\n", "b = np.asarray(adata.obsm['X_stgp_spatial'])\n", "if b.shape[0] != adata.n_obs:\n", " b = b.T\n", "if b.shape[0] != adata.n_obs:\n", " raise ValueError(f'X_stgp_spatial shape {b.shape} is not aligned with n_obs={adata.n_obs}')\n", "\n", "truth = adata.obs[TRUTH_KEY].astype('category').cat.remove_unused_categories()\n", "n_clusters = truth.cat.categories.size\n", "n_clusters = 4\n", "if n_clusters < 2:\n", " raise ValueError(f'Need at least 2 observed {TRUTH_KEY} categories for spectral clustering.')\n", "k_nn = min(max(2, int(np.round(np.sqrt(b.shape[0])))), b.shape[0] - 1)\n", "print(f'truth={TRUTH_KEY}, n_clusters={n_clusters}, k_nn={k_nn}')\n", "\n", "nn = NearestNeighbors(n_neighbors=k_nn + 1, metric='euclidean').fit(b).kneighbors(return_distance=False)[:, 1:]\n", "rows = np.repeat(np.arange(nn.shape[0]), k_nn)\n", "cols = nn.ravel()\n", "knn_graph = sp.csr_matrix((np.ones(rows.size), (rows, cols)), shape=(b.shape[0], b.shape[0]))\n", "knn_graph = knn_graph.maximum(knn_graph.T)\n", "\n", "clusterlabel = SpectralClustering(\n", " n_clusters=n_clusters,\n", " affinity='precomputed',\n", " assign_labels='kmeans',\n", " random_state=RANDOM_STATE,\n", ").fit_predict(knn_graph) + 1\n", "\n", "adata.obs['stGP_domain'] = pd.Categorical(clusterlabel.astype(str))\n", "ctab = pd.crosstab(adata.obs['stGP_domain'], truth, rownames=['stGP'], colnames=[TRUTH_KEY])\n", "display(ctab)\n", "ctab.to_csv(OUT_DIR / f'{label_suffix}_stgp_domain_{TRUTH_KEY}_crosstab.csv')" ] }, { "cell_type": "code", "execution_count": 3, "id": "61eb00f6", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:06:22.112537Z", "iopub.status.busy": "2026-07-04T23:06:22.112382Z", "iopub.status.idle": "2026-07-04T23:06:22.115666Z", "shell.execute_reply": "2026-07-04T23:06:22.115294Z" } }, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 52417 × 220\n", " obs: 'volume', 'center_x', 'center_y', 'min_x', 'min_y', 'max_x', 'max_y', 'transcript_count', 'num_detected_genes', 'barcodeCount', 'mouse_id', 'slide_id', 'cohort', 'age', 'batch', 'celltype', 'region', 'subregion', 'bbox_area', 'n_genes_by_counts', 'total_counts'\n", " var: 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'\n", " uns: 'neighbors', 'pca', 'preprocess_info', 'stgp', 'umap'\n", " obsm: 'X_pca', 'X_stgp', 'X_stgp_spatial', 'X_umap', 'spatial'\n", " varm: 'PCs'\n", " obsp: 'connectivities', 'distances'" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata_full" ] }, { "cell_type": "markdown", "id": "trajectory-heading", "metadata": {}, "source": [ "## 3. Slingshot trajectory inference" ] }, { "cell_type": "code", "execution_count": 4, "id": "run-slingshot", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:06:22.117674Z", "iopub.status.busy": "2026-07-04T23:06:22.117561Z", "iopub.status.idle": "2026-07-04T23:06:49.478450Z", "shell.execute_reply": "2026-07-04T23:06:49.477074Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cluster categories: ['1', '2', '3', '4']\n", "Start cluster: 4 (start_node=3)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\r\n", " 0%| | 0/10 [00:00\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
clusterlabelsubregionslingPseudotime_1
3615378900220100285-972CTX_L2/318.680730
3615378900220100290-972CTX_L2/321.746658
3615378900220200026-972CTX_L2/321.590094
3615378900230100262-972CTX_L1/MEN22.114490
3615378900230100309-972CTX_L2/322.111445
\n", "" ], "text/plain": [ " clusterlabel subregion slingPseudotime_1\n", "3615378900220100285-97 2 CTX_L2/3 18.680730\n", "3615378900220100290-97 2 CTX_L2/3 21.746658\n", "3615378900220200026-97 2 CTX_L2/3 21.590094\n", "3615378900230100262-97 2 CTX_L1/MEN 22.114490\n", "3615378900230100309-97 2 CTX_L2/3 22.111445" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs['slingPseudotime_1'] = sl.unified_pseudotime\n", "\n", "if sl.curves is not None and sl.cell_weights is not None:\n", " for l_idx, curve in enumerate(sl.curves):\n", " pt = curve.pseudotimes_interp.copy()\n", " weight = sl.cell_weights[:, l_idx].copy()\n", " pt[weight <= 0] = np.nan\n", " adata.obs[f'slingPseudotime_{l_idx + 1}'] = pt\n", " adata.obs[f'slingCurveWeight_{l_idx + 1}'] = weight\n", "\n", "adata.obs[['clusterlabel', TRUTH_KEY, 'slingPseudotime_1']].head()" ] }, { "cell_type": "markdown", "id": "outputs-heading", "metadata": {}, "source": [ "## 4. Visualise and save outputs" ] }, { "cell_type": "code", "execution_count": 6, "id": "plot-spatial", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:06:49.494394Z", "iopub.status.busy": "2026-07-04T23:06:49.494265Z", "iopub.status.idle": "2026-07-04T23:06:49.877960Z", "shell.execute_reply": "2026-07-04T23:06:49.877210Z" } }, "outputs": [], "source": [ "xy = np.asarray(adata.obsm['spatial'])\n", "pt = adata.obs['slingPseudotime_1'].to_numpy()\n", "\n", "x = xy[:, 0]\n", "y = -xy[:, 1]\n", "\n", "fig = plt.figure(figsize=(6, 6), constrained_layout=True)\n", "gs = fig.add_gridspec(\n", " 2, 2,\n", " width_ratios=[1, 0.08],\n", " height_ratios=[1, 1]\n", ")\n", "\n", "ax0 = fig.add_subplot(gs[0, 0])\n", "ax1 = fig.add_subplot(gs[1, 0], sharex=ax0, sharey=ax0)\n", "leg_ax = fig.add_subplot(gs[0, 1])\n", "cax = fig.add_subplot(gs[1, 1])\n", "\n", "for label in adata.obs['clusterlabel'].cat.categories:\n", " mask = adata.obs['clusterlabel'].astype(str).to_numpy() == str(label)\n", " ax0.scatter(x[mask], y[mask], s=5, linewidths=0, rasterized=True, label=str(label))\n", "\n", "scat = ax1.scatter(x, y, c=pt, s=5, cmap='viridis', linewidths=0, rasterized=True)\n", "\n", "for ax in [ax0, ax1]:\n", " ax.set_aspect('equal')\n", " ax.axis('off')\n", " ax.set_xlim(x.min(), x.max())\n", " ax.set_ylim(y.max(), y.min()) # 上下翻转两个共享空间子图\n", "\n", "leg_ax.axis('off')\n", "handles, labels = ax0.get_legend_handles_labels()\n", "leg_ax.legend(handles, labels, title='clusterlabel', markerscale=2, loc='upper left')\n", "\n", "plt.colorbar(scat, cax=cax, label='Slingshot pseudotime')\n", "\n", "fig.savefig(OUT_DIR / f'{label_suffix}_slingshot_spatial.png', dpi=300, bbox_inches='tight')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 7, "id": "plot-arrow-overlay", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:06:49.880800Z", "iopub.status.busy": "2026-07-04T23:06:49.880660Z", "iopub.status.idle": "2026-07-04T23:06:50.156160Z", "shell.execute_reply": "2026-07-04T23:06:50.155385Z" } }, "outputs": [], "source": [ "# Arrow overlay: local spatial pseudotime gradient, from lower to higher pseudotime.\n", "GRIDNUM = 10\n", "MIN_CELLS_PER_GRID = 10\n", "ARROW_LENGTH_FRAC = 0.75\n", "\n", "xy = np.asarray(adata.obsm['spatial'])\n", "pt = adata.obs['slingPseudotime_1'].to_numpy(dtype=float)\n", "valid = np.isfinite(pt) & np.all(np.isfinite(xy), axis=1)\n", "xy_valid = xy[valid]\n", "pt_valid = pt[valid]\n", "\n", "x_edges = np.linspace(xy_valid[:, 0].min(), xy_valid[:, 0].max(), GRIDNUM + 1)\n", "y_edges = np.linspace(xy_valid[:, 1].min(), xy_valid[:, 1].max(), GRIDNUM + 1)\n", "x_bin = np.clip(np.digitize(xy_valid[:, 0], x_edges) - 1, 0, GRIDNUM - 1)\n", "y_bin = np.clip(np.digitize(xy_valid[:, 1], y_edges) - 1, 0, GRIDNUM - 1)\n", "\n", "mean_pt = np.full((GRIDNUM, GRIDNUM), np.nan)\n", "mean_xy = np.full((GRIDNUM, GRIDNUM, 2), np.nan)\n", "for i in range(GRIDNUM):\n", " for j in range(GRIDNUM):\n", " m = (x_bin == i) & (y_bin == j)\n", " if int(m.sum()) >= MIN_CELLS_PER_GRID:\n", " mean_pt[i, j] = float(np.nanmean(pt_valid[m]))\n", " mean_xy[i, j] = np.nanmean(xy_valid[m], axis=0)\n", "\n", "arrow_start, arrow_vec = [], []\n", "for i in range(GRIDNUM):\n", " for j in range(GRIDNUM):\n", " if not np.isfinite(mean_pt[i, j]):\n", " continue\n", " grad = np.zeros(2, dtype=float)\n", " center = mean_xy[i, j]\n", " for di in (-1, 0, 1):\n", " for dj in (-1, 0, 1):\n", " if di == 0 and dj == 0:\n", " continue\n", " ni, nj = i + di, j + dj\n", " if ni < 0 or ni >= GRIDNUM or nj < 0 or nj >= GRIDNUM or not np.isfinite(mean_pt[ni, nj]):\n", " continue\n", " direction = mean_xy[ni, nj] - center\n", " dist = np.linalg.norm(direction)\n", " if dist > 0:\n", " grad += (mean_pt[ni, nj] - mean_pt[i, j]) * direction / dist\n", " norm = np.linalg.norm(grad)\n", " if norm > 0:\n", " arrow_start.append(center)\n", " arrow_vec.append(grad / norm)\n", "\n", "arrow_start = np.asarray(arrow_start)\n", "arrow_vec = np.asarray(arrow_vec)\n", "cell_size = min(np.diff(x_edges).mean(), np.diff(y_edges).mean())\n", "arrow_vec = arrow_vec * cell_size * ARROW_LENGTH_FRAC\n", "\n", "fig, ax = plt.subplots(figsize=(7, 7), constrained_layout=True)\n", "for label in adata.obs['clusterlabel'].cat.categories:\n", " mask = adata.obs['clusterlabel'].astype(str).to_numpy() == str(label)\n", " ax.scatter(xy[mask, 0], xy[mask, 1], s=5, linewidths=0, rasterized=True, label=str(label))\n", "if len(arrow_start) > 0:\n", " ax.quiver(arrow_start[:, 0], arrow_start[:, 1], arrow_vec[:, 0], arrow_vec[:, 1], angles='xy', scale_units='xy', scale=1, color='black', width=0.006, headwidth=4.5, headlength=6.0, headaxislength=5.0, label='pseudotime gradient')\n", "ax.set_aspect('equal')\n", "ax.invert_yaxis()\n", "ax.axis('off')\n", "ax.set_title('Spatial pseudotime direction overlay')\n", "ax.legend(title='clusterlabel', markerscale=2, bbox_to_anchor=(1.02, 1), loc='upper left')\n", "fig.savefig(OUT_DIR / f'{label_suffix}_slingshot_arrow_overlay.png', dpi=300, bbox_inches='tight')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 8, "id": "plot-umap", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:06:50.158584Z", "iopub.status.busy": "2026-07-04T23:06:50.158448Z", "iopub.status.idle": "2026-07-04T23:07:05.090692Z", "shell.execute_reply": "2026-07-04T23:07:05.089647Z" } }, "outputs": [], "source": [ "sc.pp.neighbors(adata, use_rep='X_DRM', n_neighbors=k_nn)\n", "sc.tl.umap(adata, random_state=RANDOM_STATE)\n", "sc.pl.umap(adata, color=['clusterlabel', 'slingPseudotime_1', TRUTH_KEY], size=10, wspace=0.4, show=False)\n", "plt.savefig(OUT_DIR / f'{label_suffix}_slingshot_umap.png', dpi=300, bbox_inches='tight')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "id": "save-results", "metadata": { "execution": { "iopub.execute_input": "2026-07-04T23:07:05.093128Z", "iopub.status.busy": "2026-07-04T23:07:05.092980Z", "iopub.status.idle": "2026-07-04T23:07:05.620412Z", "shell.execute_reply": "2026-07-04T23:07:05.619663Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saved trajectory AnnData: /import/home4/byual/stGP-0529/RealData_MouseBrainMERFISH/Results/stgp/Microglia/downstream_cluster/mouse_97_stgp_cluster_slingshot.h5ad\n", "Saved trajectory obs: /import/home4/byual/stGP-0529/RealData_MouseBrainMERFISH/Results/stgp/Microglia/downstream_cluster/mouse_97_stgp_cluster_slingshot_obs.csv\n" ] } ], "source": [ "OUT_H5AD = OUT_DIR / f'{label_suffix}_stgp_cluster_slingshot.h5ad'\n", "OUT_OBS = OUT_DIR / f'{label_suffix}_stgp_cluster_slingshot_obs.csv'\n", "\n", "adata.write_h5ad(OUT_H5AD, compression='gzip')\n", "adata.obs.to_csv(OUT_OBS)\n", "\n", "print(f'Saved trajectory AnnData: {OUT_H5AD.resolve()}')\n", "print(f'Saved trajectory obs: {OUT_OBS.resolve()}')" ] } ], "metadata": { "kernelspec": { "display_name": "NicheScope", "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.9.23" } }, "nbformat": 4, "nbformat_minor": 5 }