OSCR

StPedf: Cell trajectory inference of spatial transcriptomics via spatial proximity embedding and spatial density-adaptive fusion.

Code ↔ Paper

7 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 7 matches
  1. [1] § 5. Methods › 5.5. Cell density-aware adaptive weighting mechanism and low-confidence region protection ↔ TI/single_time/transition.py, lines 130–233 · score 0.84 · low confidence region, density adaptive, Low density, spatial weights, hole, connectivity
  2. [2] § 2. Results › 2.1. Overview ↔ TI/single_time/transition.py, lines 130–233 · score 0.80 · fused cost, adaptive transition matrix, cost matrix, high density, low density, spatial distance
  3. [3] § 2. Results › 2.1. Overview ↔ TI/single_time/lap.py, lines 381–439 · score 0.61 · terminal cells, vector field, transition matrix, adjacency, optimal, embedding
  4. [4] § 5. Methods › 5.4. Multi-task training objectives ↔ example.ipynb, lines 204–219 · score 0.60 · Kendall correlation, Spearman correlation, ground truth
  5. [5] § 5. Methods › 5.8. Selection of starting cells ↔ TI/single_time/velocity.py, lines 21–99 · score 0.55 · transition probability matrix, optimal transport, sequence, cells
  6. [6] § 5. Methods › 5.10. Cell velocity calculation and trajectory construction methods ↔ TI/single_time/velocity.py, lines 650–717 · score 0.53 · displacement vector, cell velocity, weighted
  7. [7] § 5. Methods › 5.10. Cell velocity calculation and trajectory construction methods ↔ TI/single_time/velocity.py, lines 212–310 · score 0.52 · numerical stability, transition probabilities, velocity, cell

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 834 lines · 25 KB · no license · 3 matches

  1. import ot
  2. import sys
  3. import scanpy as sc
  4. import pandas as pd
  5. import numpy as np
  6. from anndata import AnnData
  7. from sklearn.metrics.pairwise import euclidean_distances
  8. from sklearn.neighbors import NearestNeighbors
  9. from scipy.sparse import csr_matrix
  10. from scipy.stats import norm
  11. from numpy.random import RandomState
  12. import plotly.graph_objs as go
  13. import plotly.offline as py
  14. from ipywidgets import VBox
  15. import seaborn as sns
  16. from typing import Optional, Union, Literal, Tuple, List
  18. from .utils import nearest_neighbors, kmeans_centers
  21. def get_ot_matrix(
  22. adata: AnnData,
  23. data_type: str,
  24. alpha1: int = 1,
  25. alpha2: int = 1,
  26. random_state: Union[None, int, RandomState] = 0,
  27. ) -> np.ndarray:
  28. """
  29. Calculate transfer probabilities between cells.
  31. Using optimal transport theory based on gene expression and/or spatial location information.
  33. Parameters
  34. ----------
  35. adata
  36. An :class:`~anndata.AnnData` object.
  37. data_type
  38. The type of sequencing data.
  40. - ``'spatial'``: for the spatial transcriptome data.
  41. - ``'single-cell'``: for the single-cell sequencing data.
  43. alpha1
  44. The proportion of spatial location information.
  45. (Default: 1)
  46. alpha2
  47. The proportion of gene expression information.
  48. (Default: 1)
  49. random_state
  50. Different initial states for the pca.
  51. (Default: 0)
  53. Returns
  54. -------
  55. :class:`~numpy.ndarray`
  56. Cell transition probability matrix.
  57. """
  58. if "X_pca" not in adata.obsm:
  59. print("X_pca is not in adata.obsm, automatically do PCA first.")
  60. sc.tl.pca(adata, svd_solver="arpack", random_state=random_state)
  61. newdata = adata.obsm["X_pca"]
  62. newdata2 = newdata.copy()
  64. if data_type == "spatial":
  65. newcoor = adata.obsm["X_spatial"]
  66. newcoor2 = newcoor.copy()
  68. # calculate physical distance
  69. ed_coor = euclidean_distances(newcoor, newcoor2, squared=True)
  70. m1 = ed_coor / sum(sum(ed_coor))
  71. # calculate gene expression PCA space distance
  72. ed_gene = euclidean_distances(newdata, newdata2, squared=True)
  73. m2 = ed_gene / sum(sum(ed_gene))
  75. M = alpha1 * m1 + alpha2 * m2
  76. M /= M.max()
  77. row, col = np.diag_indices_from(M)
  78. M[row, col] = M.max() * 1000000
  80. elif data_type == "single-cell":
  81. ed = euclidean_distances(newdata, newdata2, squared=True)
  82. M = ed / sum(sum(ed))
  83. M /= M.max()
  84. row, col = np.diag_indices_from(M)
  85. M[row, col] = M.max() * 1000000
  87. else:
  88. sys.exit(
  89. "Please give the right data type, choose from 'spatial' or 'single-cell'."
  90. )
  92. a, b = (
  93. np.ones((adata.n_obs,)) / adata.n_obs,
  94. np.ones((adata.n_obs,)) / adata.n_obs,
  95. )
  96. lambd = 1e-1
  97. Gs = np.array(ot.sinkhorn(a, b, M, lambd))
  99. return Gs
  102. def set_start_cells(
  103. adata: AnnData,
  104. select_way: Literal["coordinates", "cell_type"],
  105. cell_type: Optional[str] = None,
  106. start_point: Optional[Tuple[int, int]] = None,
  107. basis: str = "spatial",
  108. split: bool = False,
  109. n_clusters: int = 2,
  110. n_neigh: int = 5,
  111. ) -> list:
  112. """
  113. Use coordinates or cell type to manually select starting cells.
  115. Parameters
  116. ----------
  117. adata
  118. An :class:`~anndata.AnnData` object.
  119. select_way
  120. Ways to select starting cells.
  122. (1) ``'cell_type'``: select by cell type.
  123. (2) ``'coordinates'``: select by coordinates.
  125. cell_type
  126. Restrict the cell type of starting cells.
  127. (Deafult: None)
  128. start_point
  129. The coordinates of the start point in 'coordinates' mode.
  130. basis
  131. The basis in `adata.obsm` to store position information.
  132. split
  133. Whether to split the specific type of cells into several small clusters according to cell density.
  134. n_clsuters
  135. The number of cluster centers after splitting.
  136. n_neigh
  137. The number of neighbors next to the start point/cluster center selected as the starting cell.
  139. Returns
  140. -------
  141. list
  142. The index number of selected starting cells.
  143. """
  144. if select_way == "coordinates":
  145. if start_point is None:
  146. raise ValueError(
  147. f"`start_point` must be specified in the 'coordinates' mode."
  148. )
  150. start_cells = nearest_neighbors(start_point, adata.obsm["X_" + basis], n_neigh)[
  151. 0
  152. ]
  154. if cell_type is not None:
  155. type_cells = np.where(adata.obs["cluster"] == cell_type)[0]
  156. start_cells = set(start_cells).intersection(set(type_cells))
  158. elif select_way == "cell_type":
  159. if cell_type is None:
  160. raise ValueError("in 'cell_type' mode, `cell_type` cannot be None.")
  162. start_cells = np.where(adata.obs["cluster"] == cell_type)[0]
  164. if split == True:
  165. mask = adata.obs["cluster"] == cell_type
  166. cell_coords = adata.obsm["X_" + basis][mask]
  167. cluster_centers = kmeans_centers(cell_coords, n_clusters=n_clusters)
  169. select_cluster_coords = adata.obsm["X_" + basis].copy()
  170. select_cluster_coords[np.logical_not(mask)] = 1e10
  171. start_cells = nearest_neighbors(
  172. cluster_centers, select_cluster_coords, n_neigh
  173. ).flatten()
  174. else:
  175. raise ValueError(f"`select_way` must choose from 'coordinates' or 'cell_type'.")
  177. return list(start_cells)
  180. def get_ptime(adata: AnnData, start_cells: list):
  181. """
  182. Get the cell pseudotime based on transition probabilities from initial cells.
  184. Parameters
  185. ----------
  186. adata
  187. An :class:`~anndata.AnnData` object.
  188. start_cells
  189. List of index numbers of starting cells.
  191. Returns
  192. -------
  193. :class:`~numpy.ndarray`
  194. Ptime correspongding to cells.
  195. """
  196. select_trans = adata.obsp["trans"][start_cells]
  197. cell_tran = np.sum(select_trans, axis=0)
  198. adata.obs["tran"] = cell_tran
  199. cell_tran_sort = list(np.argsort(cell_tran))
  200. cell_tran_sort = cell_tran_sort[::-1]
  202. ptime = pd.Series(dtype="float32", index=adata.obs.index)
  203. for i in range(adata.n_obs):
  204. ptime[cell_tran_sort[i]] = i / (adata.n_obs - 1)
  206. return ptime.values
  208. from scipy.sparse import issparse, csr_matrix
  209. import numpy as np
  210. from sklearn.neighbors import NearestNeighbors
  212. def get_neigh_trans(
  213. adata: AnnData, basis: str, n_neigh_pos: int = 10, n_neigh_embed: int = 0
  214. ):
  215. """
  216. Get transport neighbors using both spatial positions and SEDR embeddings
  218. Parameters
  219. ----------
  220. adata
  221. An :class:`~anndata.AnnData` object.
  222. basis
  223. The basis used in visualizing the cell position.
  224. n_neigh_pos
  225. Number of neighbors based on spatial positions.
  226. n_neigh_embed
  227. Number of neighbors based on SEDR embeddings.
  229. Returns
  230. -------
  231. :class:`~scipy.sparse._csr.csr_matrix`
  232. Sparse matrix of transition probabilities for selected neighbors.
  233. """
  234. if n_neigh_pos == 0 and n_neigh_embed == 0:
  235. raise ValueError(
  236. "Number of position and embedding neighbors cannot both be zero."
  237. )
  239. # Ensure numerical stability of transition matrix
  240. trans_mat = adata.obsp["trans"]
  241. if np.any(trans_mat < 0):
  242. trans_mat = np.maximum(trans_mat, 0)
  243. row_sums = trans_mat.sum(axis=1)
  244. if np.any(row_sums == 0):
  245. trans_mat += 1e-10
  246. row_sums = trans_mat.sum(axis=1)
  247. trans_mat = trans_mat / row_sums[:, np.newaxis]
  248. adata.obsp["trans"] = trans_mat
  250. # Spatial position neighbors
  251. if n_neigh_pos:
  252. nn_pos = NearestNeighbors(n_neighbors=n_neigh_pos, n_jobs=-1)
  253. nn_pos.fit(adata.obsm["X_" + basis])
  254. _, neigh_pos = nn_pos.kneighbors(adata.obsm["X_" + basis])
  255. neigh_pos = neigh_pos[:, 1:] # exclude self
  257. # SEDR embedding neighbors
  258. if n_neigh_embed:
  259. if "StPedf" not in adata.obsm:
  260. raise ValueError("StPedf embedding not found in adata.obsm")
  262. nn_embed = NearestNeighbors(n_neighbors=n_neigh_embed, n_jobs=-1)
  263. nn_embed.fit(adata.obsm["StPedf"])
  264. _, neigh_embed = nn_embed.kneighbors(adata.obsm["StPedf"])
  265. neigh_embed = neigh_embed[:, 1:] # exclude self
  267. # Build neighbor list
  268. neigh_list = []
  269. for i in range(adata.n_obs):
  270. neighbors = []
  271. if n_neigh_pos:
  272. neighbors.extend(neigh_pos[i])
  273. if n_neigh_embed:
  274. neighbors.extend(neigh_embed[i])
  276. # Deduplicate and exclude self
  277. neighbors = np.unique(neighbors)
  278. neighbors = neighbors[neighbors != i]
  279. neigh_list.append(neighbors)
  281. # Build sparse transition probability matrix
  282. indptr = [0]
  283. indices = []
  284. csr_data = []
  286. for i in range(adata.n_obs):
  287. neighbors = neigh_list[i]
  288. if len(neighbors) == 0:
  289. continue
  291. # Get transition probabilities and normalize
  292. probs = adata.obsp["trans"][i, neighbors]
  293. if issparse(probs):
  294. probs = probs.toarray().flatten()
  296. prob_sum = probs.sum()
  297. if prob_sum > 0:
  298. probs = probs / prob_sum
  299. else:
  300. probs = np.ones_like(probs) / len(probs)
  302. indices.extend(neighbors)
  303. csr_data.extend(probs)
  304. indptr.append(len(indices))
  306. trans_neigh_csr = csr_matrix(
  307. (csr_data, indices, indptr), shape=(adata.n_obs, adata.n_obs)
  308. )
  310. return trans_neigh_csr
  313. '''def get_neigh_trans(
  314. adata: AnnData, basis: str, n_neigh_pos: int = 10, n_neigh_gene: int = 0
  315. ):
  316. """
  317. Get the transport neighbors from two ways, position and/or gene expression
  319. Parameters
  320. ----------
  321. adata
  322. An :class:`~anndata.AnnData` object.
  323. basis
  324. The basis used in visualizing the cell position.
  325. n_neigh_pos
  326. Number of neighbors based on cell positions such as spatial or umap coordinates.
  327. (Default: 10)
  328. n_neigh_gene
  329. Number of neighbors based on gene expression (PCA).
  330. (Default: 0)
  332. Returns
  333. -------
  334. :class:`~scipy.sparse._csr.csr_matrix`
  335. A sparse matrix composed of transition probabilities of selected neighbor cells.
  336. """
  337. if n_neigh_pos == 0 and n_neigh_gene == 0:
  338. raise ValueError(
  339. "the number of position neighbors and gene neighbors cannot be zero at the same time."
  340. )
  342. if n_neigh_pos:
  343. nn = NearestNeighbors(n_neighbors=n_neigh_pos, n_jobs=-1)
  344. nn.fit(adata.obsm["X_" + basis])
  345. dist_pos, neigh_pos = nn.kneighbors(adata.obsm["X_" + basis])
  346. dist_pos = dist_pos[:, 1:]
  347. neigh_pos = neigh_pos[:, 1:]
  349. neigh_pos_list = []
  350. for i in range(adata.n_obs):
  351. idx = neigh_pos[i] # embedding上的邻居
  352. idx2 = neigh_pos[idx] # embedding上邻居的邻居
  353. idx2 = np.setdiff1d(idx2, i)
  355. neigh_pos_list.append(np.unique(np.concatenate([idx, idx2])))
  356. # neigh_pos_list.append(idx)
  358. if n_neigh_gene:
  359. if "X_pca" not in adata.obsm:
  360. print("X_pca is not in adata.obsm, automatically do PCA first.")
  361. sc.tl.pca(adata)
  362. sc.pp.neighbors(
  363. adata, use_rep="X_pca", key_added="X_pca", n_neighbors=n_neigh_gene
  364. )
  366. neigh_gene = adata.obsp["distances"].indices.reshape(
  367. -1, adata.uns["neighbors"]["params"]["n_neighbors"] - 1
  368. )
  370. indptr = [0]
  371. indices = []
  372. csr_data = []
  373. count = 0
  374. for i in range(adata.n_obs):
  375. if n_neigh_pos == 0:
  376. n_all = neigh_gene[i]
  377. elif n_neigh_gene == 0:
  378. n_all = neigh_pos_list[i]
  379. else:
  380. n_all = np.unique(np.concatenate([neigh_pos_list[i], neigh_gene[i]]))
  381. count += len(n_all)
  382. indptr.append(count)
  383. indices.extend(n_all)
  384. csr_data.extend(
  385. adata.obsp["trans"][i][n_all]
  386. / (adata.obsp["trans"][i][n_all].sum()) # normalize
  387. )
  389. trans_neigh_csr = csr_matrix(
  390. (csr_data, indices, indptr), shape=(adata.n_obs, adata.n_obs)
  391. )
  393. return trans_neigh_csr'''
  396. '''def get_velocity(
  397. adata: AnnData,
  398. basis: str,
  399. n_neigh_pos: int = 10,
  400. n_neigh_gene: int = 0,
  401. grid_num=50,
  402. smooth=0.5,
  403. density=1.0,
  404. ) -> tuple:
  405. """
  406. Get the velocity of each cell.
  408. The speed can be determined in terms of the cell location and/or gene expression.
  410. Parameters
  411. ----------
  412. adata
  413. An :class:`~anndata.AnnData` object.
  414. basis
  415. The label of cell coordinates, for example, `umap` or `spatial`.
  416. n_neigh_pos
  417. Number of neighbors based on cell positions such as spatial or umap coordinates.
  418. (Default: 10)
  419. n_neigh_gene
  420. Number of neighbors based on gene expression.
  421. (Default: 0)
  423. Returns
  424. -------
  425. tuple
  426. The grid coordinates and cell velocities on each grid to draw the streamplot figure.
  427. """
  428. adata.obsp["trans_neigh_csr"] = get_neigh_trans(
  429. adata, basis, n_neigh_pos, n_neigh_gene
  430. )
  432. position = adata.obsm["X_" + basis]
  433. V = np.zeros(position.shape) # 速度为2维
  435. for cell in range(adata.n_obs): # 循环每个细胞
  436. cell_u = 0.0 # 初始化细胞速度
  437. cell_v = 0.0
  438. x1 = position[cell][0] # 初始化细胞坐标
  439. y1 = position[cell][1]
  440. for neigh in adata.obsp["trans_neigh_csr"][cell].indices: # 针对每个邻居
  441. p = adata.obsp["trans_neigh_csr"][cell, neigh]
  442. if (
  443. adata.obs["ptime"][neigh] < adata.obs["ptime"][cell]
  444. ): # 若邻居的ptime小于当前的，则概率反向
  445. p = -p
  447. x2 = position[neigh][0]
  448. y2 = position[neigh][1]
  450. # 正交向量确定速度方向，乘上概率确定速度大小
  451. sub_u = p * (x2 - x1) / (np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2))
  452. sub_v = p * (y2 - y1) / (np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2))
  453. cell_u += sub_u
  454. cell_v += sub_v
  455. V[cell][0] = cell_u / adata.obsp["trans_neigh_csr"][cell].indptr[1]
  456. V[cell][1] = cell_v / adata.obsp["trans_neigh_csr"][cell].indptr[1]
  457. adata.obsm["velocity_" + basis] = V
  458. print(f"The velocity of cells store in 'velocity_{basis}'.")
  460. P_grid, V_grid = get_velocity_grid(
  461. adata,
  462. P=position,
  463. V=adata.obsm["velocity_" + basis],
  464. grid_num=grid_num,
  465. smooth=smooth,
  466. density=density,
  467. )
  468. return P_grid, V_grid'''
  469. def get_velocity(
  470. adata: AnnData,
  471. basis: str,
  472. n_neigh_pos: int = 10,
  473. n_neigh_embed: int = 5,
  474. grid_num=50,
  475. smooth=0.5,
  476. density=1.0,
  477. ) -> tuple:
  478. """
  479. Calculate cell velocity using both spatial positions and SEDR embeddings
  481. Parameters
  482. ----------
  483. adata
  484. An :class:`~anndata.AnnData` object.
  485. basis
  486. Coordinate system label (e.g., 'spatial').
  487. n_neigh_pos
  488. Number of spatial position neighbors.
  489. n_neigh_embed
  490. Number of SEDR embedding neighbors.
  491. grid_num
  492. Resolution for grid-based visualization.
  493. smooth
  494. Smoothing factor for visualization.
  495. density
  496. Density factor for visualization.
  498. Returns
  499. -------
  500. tuple
  501. Grid positions and velocities for visualization.
  502. """
  503. # Get neighbor matrix (based on spatial positions and SEDR embeddings)
  504. adata.obsp["trans_neigh_csr"] = get_neigh_trans(
  505. adata, basis, n_neigh_pos, n_neigh_embed
  506. )
  508. position = adata.obsm["X_" + basis]
  509. V = np.zeros(position.shape)
  511. # Add small value to prevent division by zero
  512. epsilon = 1e-10
  514. for cell in range(adata.n_obs):
  515. cell_u = 0.0
  516. cell_v = 0.0
  517. x1, y1 = position[cell]
  519. neighbors = adata.obsp["trans_neigh_csr"][cell].indices
  520. n_neighbors = len(neighbors)
  522. if n_neighbors == 0:
  523. V[cell] = [0.0, 0.0]
  524. continue
  526. for neigh in neighbors:
  527. p = adata.obsp["trans_neigh_csr"][cell, neigh]
  529. # Pseudotime direction correction
  530. if adata.obs["ptime"][neigh] < adata.obs["ptime"][cell]:
  531. p = -p
  533. x2, y2 = position[neigh]
  535. # Calculate distance and direction
  536. dx = x2 - x1
  537. dy = y2 - y1
  538. dist = np.sqrt(dx**2 + dy**2) + epsilon
  540. # Calculate velocity components
  541. sub_u = p * dx / dist
  542. sub_v = p * dy / dist
  544. # Numerical stability check
  545. if not np.isfinite(sub_u):
  546. sub_u = 0.0
  547. if not np.isfinite(sub_v):
  548. sub_v = 0.0
  550. cell_u += sub_u
  551. cell_v += sub_v
  553. # Average velocity
  554. V[cell][0] = cell_u / n_neighbors
  555. V[cell][1] = cell_v / n_neighbors
  557. # Final check
  558. if not np.isfinite(V[cell][0]):
  559. V[cell][0] = 0.0
  560. if not np.isfinite(V[cell][1]):
  561. V[cell][1] = 0.0
  563. # Add small noise to prevent all zeros
  564. noise = np.random.normal(0, 1e-6, V.shape)
  565. V += noise
  567. adata.obsm["velocity_" + basis] = V
  568. print(f"Cell velocities stored in 'velocity_{basis}'")
  570. # Grid processing for visualization
  571. P_grid, V_grid = get_velocity_grid(
  572. adata,
  573. P=position,
  574. V=V,
  575. grid_num=grid_num,
  576. smooth=smooth,
  577. density=density,
  578. )
  579. return P_grid, V_grid
  581. def get_2_biggest_clusters(adata):
  582. """Automatic selection of starting cells to determine the direction of cell trajectories
  584. Args:
  585. adata (anndata):
  587. Returns:
  588. tuple: Contains 2 cluster with maximum sum of transition probabilities
  589. """
  590. clusters = np.unique(adata.obs["cluster"])
  591. cluster_trans = pd.DataFrame(index=clusters, columns=clusters)
  592. for start_cluster in clusters:
  593. for end_cluster in clusters:
  594. if start_cluster == end_cluster:
  595. cluster_trans.loc[start_cluster, end_cluster] = 0
  596. continue
  597. starts = adata.obs["cluster"] == start_cluster
  598. ends = adata.obs["cluster"] == end_cluster
  599. cluster_trans.loc[start_cluster, end_cluster] = (
  600. np.sum(adata.obsp["trans"][starts][:, ends])
  601. / np.sum(starts)
  602. / np.sum(ends)
  603. )
  604. highest_2_clusters = cluster_trans.stack().astype(float).idxmax()
  605. return highest_2_clusters
  608. def auto_get_start_cluster(adata, clusters: Optional[list] = None):
  609. """
  610. Select the start cluster with the largest sum of transfer probability
  612. Parameters
  613. ----------
  614. adata
  615. Anndata
  616. clusters, list
  617. Give clusters to find, by default None, each cluster will be traversed and calculated
  619. Returns
  620. -------
  621. str
  622. One cluster with maximum sum of transition probabilities
  623. """
  624. if clusters == None:
  625. clusters = np.unique(adata.obs["cluster"])
  626. cluster_prob_sum = {}
  627. for cluster in clusters:
  628. start_cells = set_start_cells(adata, select_way="cell_type", cell_type=cluster)
  629. adata.obs["ptime"] = get_ptime(adata, start_cells)
  630. cell_time_sort = adata.obs["ptime"].values.argsort()
  632. prob_sum = 0
  633. for i in range(len(cell_time_sort) - 1):
  634. pre = cell_time_sort[i]
  635. next = cell_time_sort[i + 1]
  636. prob = adata.obsp["trans"][pre, next]
  637. prob_sum += prob
  638. cluster_prob_sum[cluster] = prob_sum
  639. highest_cluster = max(cluster_prob_sum, key=cluster_prob_sum.get)
  641. print(
  642. "The auto selecting cluster is: '"
  643. + highest_cluster
  644. + "'. If there is a large discrepancy with the known biological knowledge, please manually select the starting cluster."
  645. )
  647. return highest_cluster
  650. def get_velocity_grid(
  651. adata,
  652. P: np.ndarray,
  653. V: np.ndarray,
  654. grid_num: int = 50,
  655. smooth: float = 0.5,
  656. density: float = 1.0,
  657. ) -> tuple:
  658. """
  659. Convert cell velocity to grid velocity for streamline display
  661. The visualization of vector field borrows idea from scTour: https://github.com/LiQian-XC/sctour/blob/main/sctour.
  663. Parameters
  664. ----------
  665. P
  666. The position of cells.
  667. V
  668. The velocity of cells.
  669. smooth
  670. The factor for scale in Gaussian pdf.
  671. (Default: 0.5)
  672. density
  673. grid density
  674. (Default: 1.0)
  675. Returns
  676. ----------
  677. tuple
  678. The embedding and unitary displacement vectors in grid level.
  679. """
  680. grids = []
  681. for dim in range(P.shape[1]):
  682. m, M = np.min(P[:, dim]), np.max(P[:, dim])
  683. m = m - 0.01 * np.abs(M - m)
  684. M = M + 0.01 * np.abs(M - m)
  685. gr = np.linspace(m, M, int(grid_num * density))
  686. grids.append(gr)
  688. meshes = np.meshgrid(*grids)
  689. P_grid = np.vstack([i.flat for i in meshes]).T
  691. n_neighbors = int(P.shape[0] / grid_num)
  692. nn = NearestNeighbors(n_neighbors=n_neighbors, n_jobs=-1)
  693. nn.fit(P)
  694. dists, neighs = nn.kneighbors(P_grid)
  696. scale = np.mean([grid[1] - grid[0] for grid in grids]) * smooth
  697. weight = norm.pdf(x=dists, scale=scale)
  698. p_mass = weight.sum(1)
  700. V_grid = (V[neighs] * weight[:, :, None]).sum(1)
  701. V_grid /= np.maximum(1, p_mass)[:, None]
  703. P_grid = np.stack(grids)
  704. ns = P_grid.shape[1]
  705. V_grid = V_grid.T.reshape(2, ns, ns)
  707. mass = np.sqrt((V_grid * V_grid).sum(0))
  708. min_mass = 1e-5
  709. min_mass = np.clip(min_mass, None, np.percentile(mass, 99) * 0.01)
  710. cutoff = mass < min_mass
  712. V_grid[0][cutoff] = np.nan
  714. adata.uns["P_grid"] = P_grid
  715. adata.uns["V_grid"] = V_grid
  717. return P_grid, V_grid
  720. class Lasso:
  721. """
  722. Lasso an region of interest (ROI) based on spatial cluster.
  724. Parameters
  725. ----------
  726. adata
  727. An :class:`~anndata.AnnData` object.
  728. """
  730. __sub_index = []
  731. sub_cells = []
  733. def __init__(self, adata):
  734. self.adata = adata
  736. def vi_plot(
  737. self,
  738. basis: str = "spatial",
  739. cell_type: Optional[str] = None,
  740. ):
  741. """
  742. Plot figures.
  744. Parameters
  745. ----------
  746. basis
  747. The basis in `adata.obsm` to store position information.
  748. (Deafult: 'spatial')
  749. cell_type
  750. Restrict the cell type of starting cells.
  751. (Deafult: None)
  753. Returns
  754. -------
  755. The container of cell scatter plot and table.
  756. """
  757. cell_types = self.adata.obs["cluster"].unique()
  758. colors = sns.color_palette(n_colors=len(cell_types)).as_hex()
  759. cluster_color = dict(zip(cell_types, colors))
  760. self.adata.uns["cluster_color"] = cluster_color
  762. df = pd.DataFrame()
  763. df["group_ID"] = self.adata.obs_names
  764. df["labels"] = self.adata.obs["cluster"].values
  765. df["spatial_0"] = self.adata.obsm["X_" + basis][:, 0]
  766. df["spatial_1"] = self.adata.obsm["X_" + basis][:, 1]
  767. df["color"] = df.labels.map(self.adata.uns["cluster_color"])
  769. py.init_notebook_mode()
  771. f = go.FigureWidget(
  772. [
  773. go.Scatter(
  774. x=df["spatial_0"],
  775. y=df["spatial_1"],
  776. mode="markers",
  777. marker_color=df["color"],
  778. )
  779. ]
  780. )
  781. scatter = f.data[0]
  782. f.layout.plot_bgcolor = "rgb(255,255,255)"
  783. f.layout.autosize = False
  785. axis_dict = dict(
  786. showticklabels=True,
  787. autorange=True,
  788. )
  789. f.layout.yaxis = axis_dict
  790. f.layout.xaxis = axis_dict
  791. f.layout.width = 600
  792. f.layout.height = 600
  794. # Create a table FigureWidget that updates on selection from points in the scatter plot of f
  795. t = go.FigureWidget(
  796. [
  797. go.Table(
  798. header=dict(
  799. values=["group_ID", "labels", "spatial_0", "spatial_1"],
  800. fill=dict(color="#C2D4FF"),
  801. align=["left"] * 5,
  802. ),
  803. cells=dict(
  804. values=[
  805. df[col]
  806. for col in ["group_ID", "labels", "spatial_0", "spatial_1"]
  807. ],
  808. fill=dict(color="#F5F8FF"),
  809. align=["left"] * 5,
  810. ),
  811. )
  812. ]
  813. )
  815. def selection_fn(trace, points, selector):
  817. t.data[0].cells.values = [
  818. df.loc[points.point_inds][col]
  819. for col in ["group_ID", "labels", "spatial_0", "spatial_1"]
  820. ]
  822. Lasso.__sub_index = t.data[0].cells.values[0]
  823. Lasso.sub_cells = np.where(self.adata.obs.index.isin(Lasso.__sub_index))[0]
  825. if cell_type is not None:
  826. type_cells = np.where(self.adata.obs["cluster"] == cell_type)[0]
  827. Lasso.sub_cells = sorted(
  828. set(Lasso.sub_cells).intersection(set(type_cells))
  829. )
  831. scatter.on_selection(selection_fn)
  833. # Put everything together
  834. return VBox((f, t))

velocity.py at commit 57c3873, no license · at the source

Overview

Authors: Yuan Zhang1, Ziyan Sun1, Zhixin Shi1, Mengdi Nan1, Yuhan Fu1, Qing Ren1, Jie Gao1
ORCID iDs: Jie Gao
  1. School of Science, Jiangnan University, Wuxi, Jiangsu, China
Institutions: Jiangnan University (China)
Journal: PLoS computational biology, volume 22, issue 6, article e1014346
Dates: received 29 December 2025; accepted 19 May 2026; published online 5 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014346 · PMID 42247447 · PMCID PMC13240877 · OpenAlex W7163716196
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism)
Methods: Smoothing, state filtering, decompositions, Machine learning, Connectivity
MeSH: Spatial Transcriptomics*, Transcriptome*, Algorithms, Animals, Cell Differentiation, Computational Biology, Computer Simulation, Humans, Neural Networks, Computer (* major topic)
Topic: Single-cell and spatial transcriptomics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Innovative Research Group Project of the National Natural Science Foundation of China (11831015); National Natural Science Foundation of China (11831015, 12271216, 92370131)
Citations: not cited yet (Europe PMC); 40 references in the paper

Abstract

Spatial transcriptomics is transforming our multidimensional understanding of cellular spatial organization and its functional mechanisms in processes such as development and disease by systematically resolving the spatial heterogeneity of gene expression within tissues. To delve deeper into the dynamic processes underlying spatial expression patterns, spatial trajectory inference integrates genetic and spatial information to reconstruct the spatial developmental trajectories of cells within tissues. This approach reveals the patterns of differentiation and dynamic changes as cellular states evolve continuously along spatial axes. However, existing methods often struggle to uniformly model the complex, nonlinear interactions between high-dimensional gene expression and spatial coordinates. Here, we introduce StPedf, whose core lies in employing a neural network with a masking mechanism to capture complex nonlinear interactions between high-dimensional genes and spatial positions. It further leverages spatial proximity information as a guiding cue, dynamically and adaptively adjusting the embedding of gene and spatial information and the weighting of spatial proximity information based on spatial density. This enables trajectory inference guided by spatial information. This enables optimal transport to derive intercellular transition matrices, reconstruct cellular differentiation trajectories, and construct pseudo-spatiotemporal maps. StPedf demonstrates superior performance over existing methods on five structurally distinct simulated datasets. Using StPedf, we successfully mapped distinct lineages in the spatial trajectories of telencephalon regeneration in the Ambystoma mexicanum, multiple malignant lineages expanding within primary tumors, and developmental spatial trajectories and pseudo-spatiotemporal maps in human dorsolateral prefrontal cortex (DLPFC). StPedf significantly enhances the accuracy and interpretability of spatial trajectory inference, providing critical technical support for revealing the dynamic patterns of cellular fate transitions within tissue microenvironments.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 7 matches between paragraphs and lines of code.

YuanZ0316/StPedf

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 57c3873d95126c76f8bd28d5f141f7416c749464, 8 April 2026
Languages: Python (16), Jupyter (3)
Size: 42 files, 19 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, environment (environment.yml), 3 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (16 files), scikit-learn (14 files), Scanpy (12 files), SciPy (10 files), Matplotlib (9 files), pandas (9 files), anndata (8 files), PyTorch (8 files), seaborn (5 files), NetworkX (3 files), statsmodels (2 files), Plotly (1 file), rpy2 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
20 files
At the source: github.com/YuanZ0316/StPedf

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 19 scripts, each with its path and the digest of its content;
  • 7 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability

Reproducible codes: The StPedf method is implemented as a Python package, publicly available at GitHub: https://github.com/YuanZ0316/StPedf Data availability statement: Simulated datasets: The simulated datasets can be downloaded from GitHub: https://github.com/YuanZ0316/StPedf/tree/main/Simulated_dataset Public datasets: The original public data used in this paper can be accessed through the following links: (1) Ambystoma mexicanum brain regeneration dataset (CNP0002068): https://db.cngb.org/data_resources/project/CNP0002068 (2) Intrahepatic cholangiocarcinoma (ICC) dataset (CNP0002199): https://db.cngb.org/data_resources/project/CNP0002199 (3) Human DLPFC 10x Visium dataset: http://research.libd.org/spatialLIBD/.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Funding: added Innovative Research Group Project of the National Natural Science Foundation of China: 11831015; National Natural Science Foundation of China: 11831015, 12271216, 92370131

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 9 MeSH terms, 38 references.

Cite

This paper

Zhang, Y., Sun, Z., Shi, Z., Nan, M., Fu, Y., Ren, Q., & Gao, J. (2026). StPedf: Cell trajectory inference of spatial transcriptomics via spatial proximity embedding and spatial density-adaptive fusion. PLoS computational biology, 22(6), e1014346. https://doi.org/10.1371/journal.pcbi.1014346

BibTeX

@article{zhang2026stpedf,
author = {Zhang, Yuan and Sun, Ziyan and Shi, Zhixin and Nan, Mengdi and Fu, Yuhan and Ren, Qing and Gao, Jie},
title = {{StPedf: Cell trajectory inference of spatial transcriptomics via spatial proximity embedding and spatial density-adaptive fusion}},
journal = {PLoS computational biology},
year = {2026},
month = jun,
volume = {22},
number = {6},
pages = {e1014346},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014346},
url = {https://doi.org/10.1371/journal.pcbi.1014346},
pmid = {42247447},
pmcid = {PMC13240877}
}

RIS

TY - JOUR
AU - Zhang, Yuan
AU - Sun, Ziyan
AU - Shi, Zhixin
AU - Nan, Mengdi
AU - Fu, Yuhan
AU - Ren, Qing
AU - Gao, Jie
TI - StPedf: Cell trajectory inference of spatial transcriptomics via spatial proximity embedding and spatial density-adaptive fusion
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/06/05
VL - 22
IS - 6
SP - e1014346
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014346
UR - https://doi.org/10.1371/journal.pcbi.1014346
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014346",
"type": "article-journal",
"title": "StPedf: Cell trajectory inference of spatial transcriptomics via spatial proximity embedding and spatial density-adaptive fusion",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Zhang",
"given": "Yuan"
},
{
"family": "Sun",
"given": "Ziyan"
},
{
"family": "Shi",
"given": "Zhixin"
},
{
"family": "Nan",
"given": "Mengdi"
},
{
"family": "Fu",
"given": "Yuhan"
},
{
"family": "Ren",
"given": "Qing"
},
{
"family": "Gao",
"given": "Jie"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "6",
"page": "e1014346",
"DOI": "10.1371/journal.pcbi.1014346",
"PMID": "42247447",
"PMCID": "PMC13240877",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014346",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
5
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/advs.77003 [code]
SemanticST: A Scalable Multi-Contextual Graph Learning Framework for Uncovering Spatial Niches and Robust Multi-Sample Integration in Spatial Transcriptomics.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: rpy2, anndata, Scanpy, 7 other tools, genetics / omics, 7 references
[2] doi:10.1093/bib/bbag298 [code]
Empowering multifaceted analysis of spatial transcriptomics data with RGAST.
Journal: Briefings in bioinformatics
In common: rpy2, anndata, Scanpy, 8 other tools, genetics / omics, 5 references
[3] doi:10.21203/rs.3.rs-9676637/v1 [code]
A Comprehensive Benchmarking of Spatial Deconvolution and Domain Detection Methods across Diverse Tissues and Spatial Transcriptomic Technologies
Journal: Research Square (preprint)
In common: rpy2, anndata, Scanpy, 7 other tools, genetics / omics, 5 references
[4] doi:10.1038/s41592-026-03194-8 [code]
Beyond benchmarking: an expert-guided consensus approach to spatially aware clustering.
Journal: Nature methods
In common: rpy2, anndata, Scanpy, 7 other tools, genetics / omics, 5 references
[5] doi:10.1093/bioinformatics/btag540 [code]
Deciphering spatial heterogeneity by multimodal spatial transcriptomics modelling with SpatialModal.
Journal: Bioinformatics (Oxford, England)
In common: rpy2, anndata, Scanpy, 7 other tools, genetics / omics, 5 references
[6] doi:10.1016/j.isci.2026.116055 [code]
Mapping the transcriptional diversity of calcium signaling in the mouse and human brain.
Journal: iScience
In common: rpy2, anndata, Scanpy, 10 other tools, genetics / omics, 1 reference
[7] doi:10.1093/bioinformatics/btag430 [code]
SPIDER: spatially integrated denoising via embedding regularization with single cell supervision.
Journal: Bioinformatics (Oxford, England)
In common: rpy2, anndata, Scanpy, 6 other tools, genetics / omics, 5 references
[8] doi:10.1186/s13073-026-01704-z [code]
Gene expression profiling enables refined parcellation of cortical layers in the heterogeneous human cerebral cortex.
Journal: Genome medicine
In common: anndata, Scanpy, PyTorch, 6 other tools, db.cngb.org/data_resources/project, genetics / omics, 2 references
[9] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: rpy2, anndata, Scanpy, 8 other tools, genetics / omics, 2 references
[10] doi:10.1016/j.stemcr.2026.103015 [code]
Brain injury reactivates a developmental program driving genesis and integration of transient LGE-class interneurons.
Journal: Stem cell reports
In common: anndata, Scanpy, NetworkX, 8 other tools, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.