OSCR

ReliST: A model-agnostic risk layer for spatial transcriptomics deconvolution.

Code ↔ Paper

24 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 24 matches
  1. [1] § STAR★Methods › Method details › Pseudo-spatial known-composition benchmark ↔ scripts/run_revision_known_composition_benchmark.py, lines 970–1033 · score 0.99 · training reference h5ad, diffuse mixture, low depth, binomial thinning, uniform mixture, pseudo bulk
  2. [2] § STAR★Methods › Method details › Pseudo-spatial known-composition benchmark ↔ scripts/run_revision_known_composition_benchmark.py, lines 35–116 · score 0.99 · BrNum, donor disjoint split, cellType_k, known composition benchmark, smooth pseudo spatial, random seed
  3. [3] § Results › Known-composition pseudo-spots directly link ReliST risk to true deconvolution error ↔ scripts/run_revision_known_composition_benchmark.py, lines 35–116 · score 0.94 · BrNum, cellType_k, known composition benchmark, smooth pseudo spatial, minimum cell, DLPFC snRNA
  4. [4] § STAR★Methods › Method details › Reference perturbation stress tests ↔ scripts/run_revision_reference_perturbation.py, lines 268–347 · score 0.94 · reference composition column, coarsened truth, gene dropout, reference contamination, reference expression profile, Reference perturbation
  5. [5] § STAR★Methods › Method details › Revision analysis reproducibility ↔ scripts/build_revision_figures.py, lines 83–134 · score 0.92 · revision component ablation, revision threshold sensitivity, multimodel eval, revision reference perturbation, revision known composition, build revision
  6. [6] § STAR★Methods › Method details › Component ablation ↔ scripts/run_revision_component_ablation.py, lines 176–218 · score 0.90 · Component ablation, ablation score, single component, full risk, keep fractions, phi_local
  7. [7] § STAR★Methods › Method details › Risk score calculation ↔ scripts/run_revision_known_composition_benchmark.py, lines 910–968 · score 0.89 · local_uncertainty_risk_score, reference_risk_score, revision known composition, phi_local, phi_uncertainty, phi_reference
  8. [8] § STAR★Methods › Method details › Risk score calculation ↔ scripts/run_revision_known_composition_multimodel_eval.py, lines 486–601 · score 0.89 · local_uncertainty_risk_score, reference_risk_score, revision known composition, phi_local, phi_uncertainty, phi_reference
  9. [9] § STAR★Methods › Method details › Component ablation ↔ scripts/run_revision_reference_perturbation_component_ablation.py, lines 111–205 · score 0.83 · Component ablation, phi_local, phi_uncertainty, phi_reference, full risk, Delta
  10. [10] § STAR★Methods › Method details › Threshold sensitivity and use-case guidance ↔ scripts/run_revision_threshold_sensitivity.py, lines 203–246 · score 0.80 · Threshold sensitivity, coverage risk curves, illustrative review budgets, risk fractions, policies, error reduction
  11. [11] § Results › Known-composition pseudo-spots directly link ReliST risk to true deconvolution error ↔ scripts/run_revision_threshold_sensitivity.py, lines 203–246 · score 0.80 · threshold sensitivity, coverage risk curve, illustrative review budgets, risk_score, keep fractions, cutoffs
  12. [12] § Results › DLPFC supports risk scores as proxy-anchored diagnostic signals ↔ src/st_risk/reporting/anchored_validation.py, lines 399–455 · score 0.74 · Retained proxy ratios, signature_residual, marker_discordance, coverage risk curve, layer_guess, reliability curve
  13. [13] § STAR★Methods › Method details › Confidence and uncertainty baseline comparison ↔ scripts/build_revision_figures.py, lines 1–80 · score 0.73 · Cross model disagreement, abundance entropy risk, inverse top, ReliST, ambiguity, margin
  14. [14] § STAR★Methods › Method details › Revision analysis reproducibility ↔ scripts/run_revision_known_composition_multimodel_eval.py, lines 55–96 · score 0.73 · multimodel eval, revision known composition, seed, benchmark
  15. [15] § Results › Common-feature controls and goal-based decision use define support tiers ↔ src/st_risk/reporting/decision_support.py, lines 266–361 · score 0.70 · decision support bundle, low disagreement, high disagreement, risk maps, trusted, filters
  16. [16] § Results › Known-composition pseudo-spots directly link ReliST risk to true deconvolution error ↔ scripts/run_revision_reference_perturbation.py, lines 268–347 · score 0.66 · gene dropout, reference contamination, reference perturbation, coarsening, excitatory, inhibitory
  17. [17] § Results › DLPFC supports risk scores as proxy-anchored diagnostic signals ↔ src/st_risk/reporting/anchored_validation.py, lines 1–31 · score 0.66 · Maynard prototype distance, anchored validation, signature residual, discordance, proxies, layer
  18. [18] § Results › DLPFC supports risk scores as proxy-anchored diagnostic signals ↔ scripts/run_revision_uncertainty_baseline_eval.py, lines 155–179 · score 0.64 · signature_residual, marker_discordance, layer_guess, Maynard, AUCs, anchored
  19. [19] § STAR★Methods › Quantification and statistical analysis › Validation metrics, resampling, and artifact-aware interpretation ↔ src/st_risk/reporting/anchored_validation.py, lines 1–31 · score 0.63 · Maynard prototype distance, signature residual, discordance, anchored, validation, layer
  20. [20] § STAR★Methods › Method details › Base models and canonical output contract ↔ scripts/build_revision_figures.py, lines 389–430 · score 0.57 · risk axis score, contract fairness, AUC, proxy
  21. [21] § STAR★Methods › Method details › Reproducible risk scoring ↔ scripts/run_risk_scoring.py, lines 163–236 · score 0.54 · reference_subsampling_instability, phi_reference, discordance, subsets, score, risk
  22. [22] § STAR★Methods › Method details › Reproducible risk scoring ↔ src/st_risk/eval/reference_eval.py, lines 180–207 · score 0.54 · reference_subsampling_instability, phi_reference, discordance, subsets, score, risk
  23. [23] § Results › Known-composition pseudo-spots directly link ReliST risk to true deconvolution error ↔ scripts/run_revision_component_ablation.py, lines 176–218 · score 0.52 · component ablation, known composition benchmark, phi_reference, reference perturbation, matched, pseudo
  24. [24] § STAR★Methods › Method details › Reproducible risk scoring ↔ scripts/run_revision_known_composition_benchmark.py, lines 848–908 · score 0.51 · inverse distance weighting, matrix, neighborhoods, abundance, spot, scoring

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 · 1,037 lines · 45 KB · MIT · 5 matches

  1. from __future__ import annotations
  3. import argparse
  4. import json
  5. from pathlib import Path
  7. import anndata as ad
  8. import numpy as np
  9. import pandas as pd
  10. from scipy import sparse
  11. from scipy.optimize import nnls
  12. from scipy.stats import pearsonr, spearmanr
  13. from sklearn.decomposition import TruncatedSVD
  14. from sklearn.metrics import average_precision_score, roc_auc_score
  16. from st_risk.eval.reference_eval import (
  17. compute_reference_marker_scores,
  18. reference_marker_discordance_proxy,
  19. reference_signature_residual_proxy,
  20. reference_subsampling_instability,
  21. select_signature_markers,
  22. subset_markers,
  23. )
  24. from st_risk.models.base import BaseSpatialModelOutput
  25. from st_risk.paths import ensure_results_layout, project_root, results_file, set_selected_run
  26. from st_risk.risk.features import ambiguity_score, build_feature_table
  27. from st_risk.risk.neighbors import inverse_distance_weights, knn_indices
  28. from st_risk.risk.score import grouped_zscore, sigmoid
  29. from st_risk.risk.stability import gene_subsample_stability, ridge_project_celltype_proportions, row_normalize
  32. DEFAULT_RUN_ID = "2026-06-20-dlpfc-known-composition-v2-donor-disjoint"
  35. def parse_args() -> argparse.Namespace:
  36. default_source_root = project_root() / "results" / "dlpfc_rctd"
  37. selected_run = (default_source_root / "selected_run.txt").read_text(encoding="utf-8").strip()
  38. default_source_run = default_source_root / "runs" / selected_run
  40. parser = argparse.ArgumentParser(
  41. description="Build a DLPFC pseudo-spot known-composition benchmark for revision analyses."
  42. )
  43. parser.add_argument(
  44. "--reference-h5ad",
  45. type=Path,
  46. default=project_root() / "data" / "Human DLPFC" / "ready" / "dlpfc_snrna_ref.h5ad",
  47. help="snRNA reference h5ad used to draw pseudo-spots.",
  48. )
  49. parser.add_argument(
  50. "--reference-label-column",
  51. default="cellType_k",
  52. help="Reference cell-type column used as the known composition label.",
  53. )
  54. parser.add_argument(
  55. "--split-column",
  56. default="BrNum",
  57. help="Reference obs column used for donor-disjoint train/simulation split.",
  58. )
  59. parser.add_argument(
  60. "--split-mode",
  61. default="donor_disjoint",
  62. choices=("donor_disjoint", "cell_random"),
  63. help="How to split reference cells into signature-building and held-out simulation pools.",
  64. )
  65. parser.add_argument(
  66. "--min-split-cells-per-type",
  67. type=int,
  68. default=20,
  69. help="Minimum cells per retained type on each side of a donor-disjoint split.",
  70. )
  71. parser.add_argument(
  72. "--source-run-dir",
  73. type=Path,
  74. default=default_source_run,
  75. help="Existing DLPFC run providing the selected gene list and optional cell-type order.",
  76. )
  77. parser.add_argument(
  78. "--output-root",
  79. type=Path,
  80. default=project_root() / "results" / "revision_known_composition_benchmark",
  81. help="Result root for this revision benchmark.",
  82. )
  83. parser.add_argument("--run-id", default=DEFAULT_RUN_ID, help="Run id under output-root/runs/.")
  84. parser.add_argument("--n-spots", type=int, default=1200, help="Number of pseudo-spots to generate.")
  85. parser.add_argument("--n-regions", type=int, default=6, help="Number of smooth pseudo-spatial regions.")
  86. parser.add_argument("--cells-per-spot-min", type=int, default=4, help="Minimum cells per pseudo-spot.")
  87. parser.add_argument("--cells-per-spot-max", type=int, default=10, help="Maximum cells per pseudo-spot.")
  88. parser.add_argument("--min-cells-per-type", type=int, default=80, help="Minimum reference cells per retained type.")
  89. parser.add_argument("--train-fraction", type=float, default=0.5, help="Fraction of cells used to build signatures.")
  90. parser.add_argument("--layer", default="counts", help="Reference h5ad layer containing raw counts.")
  91. parser.add_argument("--marker-top-k", type=int, default=25, help="Top signature markers per cell type.")
  92. parser.add_argument("--min-positive-markers", type=int, default=10, help="Minimum positive markers per cell type.")
  93. parser.add_argument(
  94. "--marker-subset-mode",
  95. default="odd",
  96. choices=("all", "odd", "even", "top_half", "bottom_half"),
  97. help="Marker subset used for reference_subsampling_instability.",
  98. )
  99. parser.add_argument("--reference-repeats", type=int, default=8, help="Marker subsampling repeats.")
  100. parser.add_argument("--reference-fraction", type=float, default=0.5, help="Marker subsampling fraction.")
  101. parser.add_argument(
  102. "--projection-method",
  103. default="nnls",
  104. choices=("nnls", "ridge"),
  105. help="Nonnegative projection method used as the lightweight deconvolution baseline.",
  106. )
  107. parser.add_argument(
  108. "--stability-repeats",
  109. type=int,
  110. default=0,
  111. help="Optional gene-subsampling stability repeats. Default 0 keeps this revision benchmark reference-centered.",
  112. )
  113. parser.add_argument("--stability-gene-fraction", type=float, default=0.8, help="Gene fraction if stability is enabled.")
  114. parser.add_argument("--ridge-lambda", type=float, default=1e-3, help="Ridge penalty for signature projection.")
  115. parser.add_argument("--random-state", type=int, default=20260620, help="Random seed.")
  116. return parser.parse_args()
  119. def _read_gene_list(source_run_dir: Path) -> list[str]:
  120. used_genes_path = source_run_dir / "tables" / "base_model_used_genes.csv"
  121. if not used_genes_path.exists():
  122. used_genes_path = source_run_dir / "tables" / "cell2location_used_genes.csv"
  123. if used_genes_path.exists():
  124. return pd.read_csv(used_genes_path)["gene"].astype(str).tolist()
  126. signatures_path = source_run_dir / "tables" / "reference_signatures_means.csv"
  127. if not signatures_path.exists():
  128. raise FileNotFoundError(
  129. f"Could not find base_model_used_genes.csv or reference_signatures_means.csv under {source_run_dir}"
  130. )
  131. return pd.read_csv(signatures_path, index_col=0).index.astype(str).tolist()
  134. def _read_source_celltype_order(source_run_dir: Path) -> list[str] | None:
  135. signatures_path = source_run_dir / "tables" / "reference_signatures_means.csv"
  136. if not signatures_path.exists():
  137. return None
  138. signatures = pd.read_csv(signatures_path, index_col=0, nrows=1)
  139. return signatures.columns.astype(str).tolist()
  142. def _normalize_log_cp10k(matrix: np.ndarray | sparse.spmatrix) -> np.ndarray:
  143. if sparse.issparse(matrix):
  144. matrix = matrix.toarray()
  145. values = np.asarray(matrix, dtype=np.float32)
  146. library = values.sum(axis=1, keepdims=True)
  147. safe_library = np.where(np.isclose(library, 0.0), 1.0, library)
  148. return np.log1p((values / safe_library) * 1e4).astype(np.float32)
  151. def _as_dense_vector(matrix: np.ndarray | sparse.spmatrix) -> np.ndarray:
  152. if sparse.issparse(matrix):
  153. matrix = matrix.toarray()
  154. return np.asarray(matrix, dtype=np.float32).reshape(-1)
  157. def _select_reference_matrix(
  158. reference_h5ad: Path,
  159. *,
  160. layer: str,
  161. label_col: str,
  162. split_col: str | None,
  163. requested_genes: list[str],
  164. source_celltype_order: list[str] | None,
  165. min_cells_per_type: int,
  166. ) -> tuple[pd.DataFrame, sparse.spmatrix | np.ndarray, list[str], list[str], pd.DataFrame]:
  167. adata = ad.read_h5ad(reference_h5ad, backed="r")
  168. if label_col not in adata.obs.columns:
  169. raise KeyError(f"{label_col} is not present in {reference_h5ad}")
  170. if split_col and split_col not in adata.obs.columns:
  171. raise KeyError(f"{split_col} is not present in {reference_h5ad}")
  172. if layer not in adata.layers:
  173. raise KeyError(f"{layer} is not present in layers of {reference_h5ad}")
  175. var_lookup = {str(gene).lower(): str(gene) for gene in adata.var_names.astype(str)}
  176. genes = [var_lookup[str(gene).lower()] for gene in requested_genes if str(gene).lower() in var_lookup]
  177. if len(genes) < 50:
  178. raise ValueError(f"Only {len(genes)} requested genes were found in {reference_h5ad}; need at least 50.")
  180. labels_all = adata.obs[label_col].astype(str)
  181. counts_by_type = labels_all.value_counts()
  182. source_order = pd.Index(source_celltype_order or [], dtype=str)
  183. if source_celltype_order is None:
  184. celltypes = counts_by_type.loc[counts_by_type >= min_cells_per_type].index.astype(str).tolist()
  185. else:
  186. celltypes = [
  187. celltype
  188. for celltype in source_celltype_order
  189. if celltype in counts_by_type.index and int(counts_by_type[celltype]) >= min_cells_per_type
  190. ]
  191. if len(celltypes) < 3:
  192. raise ValueError("At least three retained cell types are required for this benchmark.")
  194. inclusion_table = (
  195. counts_by_type.rename_axis("celltype")
  196. .reset_index(name="n_reference_cells")
  197. .assign(
  198. in_source_celltype_order=lambda frame: frame["celltype"].isin(source_order).astype(bool),
  199. passes_min_cells=lambda frame: frame["n_reference_cells"].ge(min_cells_per_type),
  200. retained_before_split=lambda frame: frame["celltype"].isin(celltypes).astype(bool),
  201. )
  202. .sort_values(["retained_before_split", "n_reference_cells"], ascending=[False, False])
  203. .reset_index(drop=True)
  204. )
  206. keep_mask = labels_all.isin(celltypes).to_numpy()
  207. obs_columns = [label_col]
  208. if split_col:
  209. obs_columns.append(split_col)
  210. obs = adata.obs.loc[keep_mask, obs_columns].copy()
  211. obs[label_col] = obs[label_col].astype(str)
  212. if split_col:
  213. obs[split_col] = obs[split_col].astype(str)
  214. matrix = adata[keep_mask, genes].layers[layer]
  215. if sparse.issparse(matrix):
  216. matrix = matrix.tocsr()
  217. else:
  218. matrix = np.asarray(matrix, dtype=np.float32)
  219. if hasattr(adata, "file") and adata.file is not None:
  220. adata.file.close()
  221. return obs, matrix, genes, celltypes, inclusion_table
  224. def _split_reference_cells(
  225. obs: pd.DataFrame,
  226. *,
  227. label_col: str,
  228. celltypes: list[str],
  229. train_fraction: float,
  230. rng: np.random.Generator,
  231. ) -> tuple[dict[str, np.ndarray], dict[str, np.ndarray]]:
  232. train: dict[str, np.ndarray] = {}
  233. simulate: dict[str, np.ndarray] = {}
  234. fraction = float(np.clip(train_fraction, 0.2, 0.8))
  235. labels = obs[label_col].astype(str).to_numpy()
  236. for celltype in celltypes:
  237. positions = np.flatnonzero(labels == celltype)
  238. shuffled = rng.permutation(positions)
  239. n_train = int(round(len(shuffled) * fraction))
  240. n_train = min(max(n_train, 20), len(shuffled) - 1)
  241. train[celltype] = np.sort(shuffled[:n_train])
  242. simulate[celltype] = np.sort(shuffled[n_train:])
  243. if simulate[celltype].size == 0:
  244. simulate[celltype] = train[celltype]
  245. return train, simulate
  248. def _split_reference_cells_donor_disjoint(
  249. obs: pd.DataFrame,
  250. *,
  251. label_col: str,
  252. split_col: str,
  253. celltypes: list[str],
  254. train_fraction: float,
  255. min_cells_per_side: int,
  256. rng: np.random.Generator,
  257. ) -> tuple[dict[str, np.ndarray], dict[str, np.ndarray], list[str], dict[str, object]]:
  258. labels = obs[label_col].astype(str).to_numpy()
  259. units = obs[split_col].astype(str).to_numpy()
  260. unique_units = np.asarray(pd.Index(units).unique().astype(str), dtype=object)
  261. if unique_units.size < 2:
  262. raise ValueError(f"Need at least two unique {split_col} values for donor-disjoint split.")
  264. fraction = float(np.clip(train_fraction, 0.2, 0.8))
  265. n_train_units = int(round(unique_units.size * fraction))
  266. n_train_units = min(max(n_train_units, 1), unique_units.size - 1)
  268. best: tuple[list[str], set[str], set[str]] | None = None
  269. best_score = (-1, -1)
  270. for _ in range(500):
  271. shuffled = rng.permutation(unique_units)
  272. train_units = set(map(str, shuffled[:n_train_units]))
  273. simulate_units = set(map(str, shuffled[n_train_units:]))
  274. train_mask = np.asarray([unit in train_units for unit in units], dtype=bool)
  275. simulate_mask = np.asarray([unit in simulate_units for unit in units], dtype=bool)
  276. retained: list[str] = []
  277. total_cells = 0
  278. for celltype in celltypes:
  279. label_mask = labels == celltype
  280. n_train = int((label_mask & train_mask).sum())
  281. n_sim = int((label_mask & simulate_mask).sum())
  282. if n_train >= min_cells_per_side and n_sim >= min_cells_per_side:
  283. retained.append(celltype)
  284. total_cells += n_train + n_sim
  285. score = (len(retained), total_cells)
  286. if score > best_score:
  287. best = (retained, train_units, simulate_units)
  288. best_score = score
  289. if len(retained) == len(celltypes):
  290. break
  292. if best is None:
  293. raise RuntimeError("Could not construct a donor-disjoint split.")
  294. retained_celltypes, train_units, simulate_units = best
  295. if len(retained_celltypes) < 3:
  296. raise ValueError(
  297. "Fewer than three cell types passed the donor-disjoint split filters; "
  298. "try lowering --min-split-cells-per-type or using --split-mode cell_random."
  299. )
  301. train: dict[str, np.ndarray] = {}
  302. simulate: dict[str, np.ndarray] = {}
  303. train_mask = np.asarray([unit in train_units for unit in units], dtype=bool)
  304. simulate_mask = np.asarray([unit in simulate_units for unit in units], dtype=bool)
  305. for celltype in retained_celltypes:
  306. label_mask = labels == celltype
  307. train[celltype] = np.flatnonzero(label_mask & train_mask)
  308. simulate[celltype] = np.flatnonzero(label_mask & simulate_mask)
  310. metadata = {
  311. "split_mode": "donor_disjoint",
  312. "split_column": split_col,
  313. "train_units": sorted(train_units),
  314. "simulate_units": sorted(simulate_units),
  315. "n_train_units": int(len(train_units)),
  316. "n_simulate_units": int(len(simulate_units)),
  317. "min_split_cells_per_type": int(min_cells_per_side),
  318. "n_retained_celltypes_after_split": int(len(retained_celltypes)),
  319. "dropped_after_split": [celltype for celltype in celltypes if celltype not in retained_celltypes],
  320. }
  321. return train, simulate, retained_celltypes, metadata
  324. def _augment_inclusion_table(
  325. inclusion_table: pd.DataFrame,
  326. *,
  327. celltypes_after_split: list[str],
  328. train_cells: dict[str, np.ndarray],
  329. simulate_cells: dict[str, np.ndarray],
  330. ) -> pd.DataFrame:
  331. table = inclusion_table.copy()
  332. train_counts = {celltype: int(len(indices)) for celltype, indices in train_cells.items()}
  333. simulate_counts = {celltype: int(len(indices)) for celltype, indices in simulate_cells.items()}
  334. table["n_train_signature_cells"] = table["celltype"].map(train_counts).fillna(0).astype(int)
  335. table["n_heldout_simulation_cells"] = table["celltype"].map(simulate_counts).fillna(0).astype(int)
  336. table["retained_after_split"] = table["celltype"].isin(celltypes_after_split).astype(bool)
  337. table["exclusion_reason"] = "retained"
  338. table.loc[~table["passes_min_cells"], "exclusion_reason"] = "below_min_cells_per_type"
  339. table.loc[
  340. table["retained_before_split"] & ~table["retained_after_split"],
  341. "exclusion_reason",
  342. ] = "insufficient_train_or_heldout_cells_after_split"
  343. table.loc[~table["retained_before_split"] & table["passes_min_cells"], "exclusion_reason"] = "not_in_source_celltype_order"
  344. return table
  347. def _compute_signatures(
  348. matrix: sparse.spmatrix | np.ndarray,
  349. *,
  350. train_cells: dict[str, np.ndarray],
  351. genes: list[str],
  352. celltypes: list[str],
  353. ) -> pd.DataFrame:
  354. columns = {}
  355. for celltype in celltypes:
  356. normalized = _normalize_log_cp10k(matrix[train_cells[celltype], :])
  357. columns[celltype] = normalized.mean(axis=0)
  358. signatures = pd.DataFrame(columns, index=genes, dtype=float)
  359. return signatures
  362. def _build_region_prototypes(celltypes: list[str], *, n_regions: int, rng: np.random.Generator) -> pd.DataFrame:
  363. n_types = len(celltypes)
  364. rows = []
  365. for region_id in range(n_regions):
  366. alpha = np.full(n_types, 0.08, dtype=float)
  367. dominant_count = min(4, n_types)
  368. dominant = rng.choice(n_types, size=dominant_count, replace=False)
  369. alpha[dominant] = 3.0
  370. prototype = rng.dirichlet(alpha)
  371. rows.append(prototype)
  372. return pd.DataFrame(rows, columns=celltypes, index=[f"region_{i + 1}" for i in range(n_regions)])
  375. def _smooth_region_composition(
  376. *,
  377. y: int,
  378. height: int,
  379. prototypes: pd.DataFrame,
  380. rng: np.random.Generator,
  381. concentration: float,
  382. ) -> tuple[str, np.ndarray]:
  383. n_regions = prototypes.shape[0]
  384. scaled = ((y + 0.5) / max(height, 1)) * n_regions
  385. lower = int(np.floor(scaled))
  386. lower = min(max(lower, 0), n_regions - 1)
  387. upper = min(lower + 1, n_regions - 1)
  388. mix = scaled - lower
  389. base = (1.0 - mix) * prototypes.iloc[lower].to_numpy(dtype=float) + mix * prototypes.iloc[upper].to_numpy(dtype=float)
  390. alpha = np.clip(base * concentration, 0.02, None)
  391. composition = rng.dirichlet(alpha)
  392. return str(prototypes.index[lower]), composition
  395. def _sample_pseudo_spots(
  396. matrix: sparse.spmatrix | np.ndarray,
  397. *,
  398. simulate_cells: dict[str, np.ndarray],
  399. celltypes: list[str],
  400. genes: list[str],
  401. markers: dict[str, list[str]],
  402. n_spots: int,
  403. n_regions: int,
  404. cells_per_spot_min: int,
  405. cells_per_spot_max: int,
  406. rng: np.random.Generator,
  407. ) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, np.ndarray]:
  408. width = int(np.ceil(np.sqrt(n_spots)))
  409. height = int(np.ceil(n_spots / width))
  410. prototypes = _build_region_prototypes(celltypes, n_regions=n_regions, rng=rng)
  411. gene_to_idx = {gene: idx for idx, gene in enumerate(genes)}
  413. scenario_names = np.asarray(["clean", "low_depth", "marker_dropout", "diffuse_mixture"], dtype=object)
  414. scenario_probs = np.asarray([0.55, 0.15, 0.15, 0.15], dtype=float)
  416. spot_rows: list[dict[str, object]] = []
  417. true_rows: list[np.ndarray] = []
  418. pseudo_counts = np.zeros((n_spots, len(genes)), dtype=np.float32)
  419. for spot_id in range(n_spots):
  420. x = spot_id % width
  421. y = spot_id // width
  422. scenario = str(rng.choice(scenario_names, p=scenario_probs))
  423. concentration = 12.0 if scenario == "diffuse_mixture" else 80.0
  424. region_name, composition = _smooth_region_composition(
  425. y=y,
  426. height=height,
  427. prototypes=prototypes,
  428. rng=rng,
  429. concentration=concentration,
  430. )
  431. if scenario == "diffuse_mixture":
  432. composition = row_normalize((0.65 * composition + 0.35 / len(celltypes))[None, :])[0]
  434. n_cells = int(rng.integers(cells_per_spot_min, cells_per_spot_max + 1))
  435. type_counts = rng.multinomial(n_cells, composition)
  436. if type_counts.sum() == 0:
  437. type_counts[int(np.argmax(composition))] = n_cells
  438. true_fraction = type_counts / max(type_counts.sum(), 1)
  440. selected_cell_positions: list[int] = []
  441. for celltype, count in zip(celltypes, type_counts, strict=True):
  442. if count <= 0:
  443. continue
  444. pool = simulate_cells[celltype]
  445. selected = rng.choice(pool, size=int(count), replace=True)
  446. selected_cell_positions.extend(int(value) for value in selected.tolist())
  448. if selected_cell_positions:
  449. spot_counts = _as_dense_vector(matrix[selected_cell_positions, :].sum(axis=0))
  450. else:
  451. spot_counts = np.zeros(len(genes), dtype=np.float32)
  453. if scenario == "low_depth":
  454. keep_probability = float(rng.uniform(0.15, 0.45))
  455. spot_counts = rng.binomial(np.maximum(spot_counts, 0).astype(np.int64), keep_probability).astype(np.float32)
  456. elif scenario == "marker_dropout":
  457. dominant_type = celltypes[int(np.argmax(true_fraction))]
  458. marker_genes = [gene for gene in markers.get(dominant_type, []) if gene in gene_to_idx]
  459. if marker_genes:
  460. marker_idx = np.asarray([gene_to_idx[gene] for gene in marker_genes], dtype=int)
  461. spot_counts[marker_idx] = rng.binomial(
  462. np.maximum(spot_counts[marker_idx], 0).astype(np.int64),
  463. 0.25,
  464. ).astype(np.float32)
  466. pseudo_counts[spot_id, :] = spot_counts
  467. spot_name = f"pseudo_spot_{spot_id + 1:05d}"
  468. spot_rows.append(
  469. {
  470. "spot_id": spot_name,
  471. "sample_id": "dlpfc_known_composition",
  472. "x_spatial": float(x),
  473. "y_spatial": float(y),
  474. "pseudo_region": region_name,
  475. "scenario": scenario,
  476. "n_cells": int(n_cells),
  477. "library_size": float(spot_counts.sum()),
  478. }
  479. )
  480. true_rows.append(true_fraction)
  482. spot_table = pd.DataFrame(spot_rows).set_index("spot_id")
  483. true_abundance = pd.DataFrame(true_rows, index=spot_table.index, columns=celltypes)
  484. return spot_table, true_abundance, prototypes, pseudo_counts
  487. def _compute_local_heterogeneity(expression: np.ndarray, neighbors: np.ndarray, *, n_components: int = 12) -> np.ndarray:
  488. n_components = min(n_components, max(2, expression.shape[1] - 1), max(2, expression.shape[0] - 1))
  489. svd = TruncatedSVD(n_components=n_components, random_state=0)
  490. embedding = svd.fit_transform(expression)
  491. heterogeneity = np.zeros(embedding.shape[0], dtype=float)
  492. for i, row in enumerate(neighbors):
  493. valid = row[row >= 0]
  494. if len(valid) == 0:
  495. heterogeneity[i] = 1.0
  496. continue
  497. local = embedding[valid]
  498. center = local.mean(axis=0, keepdims=True)
  499. heterogeneity[i] = float(np.mean(np.sum((local - center) ** 2, axis=1))) + 1e-6
  500. return heterogeneity
  503. def _project_abundance(
  504. expression: np.ndarray,
  505. signatures: pd.DataFrame,
  506. *,
  507. method: str,
  508. ridge_lambda: float,
  509. ) -> np.ndarray:
  510. normalized_method = method.strip().lower()
  511. signature_values = signatures.to_numpy(dtype=np.float32)
  512. if normalized_method == "ridge":
  513. return ridge_project_celltype_proportions(expression, signature_values, ridge_lambda=ridge_lambda)
  514. if normalized_method != "nnls":
  515. raise ValueError(f"Unsupported projection method: {method}")
  517. projected = np.zeros((expression.shape[0], signature_values.shape[1]), dtype=np.float32)
  518. for idx, y in enumerate(expression):
  519. weights, _ = nnls(signature_values, np.asarray(y, dtype=np.float64), maxiter=1000)
  520. projected[idx, :] = weights.astype(np.float32)
  521. return row_normalize(projected)
  524. def _combine_any_features(
  525. table: pd.DataFrame,
  526. weights: dict[str, float],
  527. *,
  528. groups: pd.Series | None = None,
  529. ) -> pd.Series:
  530. linear = np.zeros(table.shape[0], dtype=float)
  531. total_weight = 0.0
  532. for name, weight in weights.items():
  533. if name not in table.columns or np.isclose(float(weight), 0.0):
  534. continue
  535. linear += float(weight) * grouped_zscore(table[name].to_numpy(dtype=float), groups=groups)
  536. total_weight += abs(float(weight))
  537. if np.isclose(total_weight, 0.0):
  538. raise ValueError("At least one non-zero available feature is required.")
  539. return pd.Series(sigmoid(linear / total_weight), index=table.index)
  542. def _abundance_baselines(predicted: pd.DataFrame) -> pd.DataFrame:
  543. values = predicted.to_numpy(dtype=float)
  544. row_sums = values.sum(axis=1, keepdims=True)
  545. probs = np.divide(values, row_sums, out=np.zeros_like(values), where=row_sums > 0)
  546. sorted_probs = np.sort(probs, axis=1)
  547. top1 = sorted_probs[:, -1] if probs.shape[1] else np.zeros(probs.shape[0], dtype=float)
  548. top2 = sorted_probs[:, -2] if probs.shape[1] > 1 else np.zeros(probs.shape[0], dtype=float)
  549. if probs.shape[1] <= 1:
  550. entropy = np.zeros(probs.shape[0], dtype=float)
  551. else:
  552. entropy = -np.sum(probs * np.log(probs + 1e-12), axis=1) / np.log(probs.shape[1])
  553. return pd.DataFrame(
  554. {
  555. "abundance_entropy_risk": entropy,
  556. "inverse_top1_margin": 1.0 - (top1 - top2),
  557. "inverse_max_abundance": 1.0 - top1,
  558. },
  559. index=predicted.index,
  560. )
  563. def _error_table(predicted: pd.DataFrame, truth: pd.DataFrame) -> pd.DataFrame:
  564. pred = predicted.loc[truth.index, truth.columns].to_numpy(dtype=float)
  565. true = truth.to_numpy(dtype=float)
  566. absolute = np.abs(pred - true)
  567. rmse = np.sqrt(np.mean((pred - true) ** 2, axis=1))
  568. numerator = (pred * true).sum(axis=1)
  569. denominator = np.linalg.norm(pred, axis=1) * np.linalg.norm(true, axis=1)
  570. cosine = np.divide(numerator, denominator, out=np.zeros_like(numerator), where=denominator > 0)
  571. return pd.DataFrame(
  572. {
  573. "l1_error": absolute.sum(axis=1),
  574. "total_variation_error": 0.5 * absolute.sum(axis=1),
  575. "rmse_error": rmse,
  576. "cosine_distance": 1.0 - np.clip(cosine, -1.0, 1.0),
  577. "dominant_mismatch": predicted.idxmax(axis=1).ne(truth.idxmax(axis=1)).astype(int).to_numpy(),
  578. },
  579. index=truth.index,
  580. )
  583. def _safe_corr(score: pd.Series, error: pd.Series, *, method: str) -> tuple[float, float]:
  584. valid = pd.concat([score, error], axis=1).dropna()
  585. if valid.shape[0] < 3 or valid.iloc[:, 0].nunique() <= 1 or valid.iloc[:, 1].nunique() <= 1:
  586. return np.nan, np.nan
  587. if method == "spearman":
  588. stat, pvalue = spearmanr(valid.iloc[:, 0], valid.iloc[:, 1])
  589. elif method == "pearson":
  590. stat, pvalue = pearsonr(valid.iloc[:, 0], valid.iloc[:, 1])
  591. else:
  592. raise ValueError(method)
  593. return float(stat), float(pvalue)
  596. def _safe_auc(score: pd.Series, labels: pd.Series) -> tuple[float, float]:
  597. valid = pd.concat([score, labels], axis=1).dropna()
  598. if valid.shape[0] < 3 or valid.iloc[:, 1].nunique() < 2 or valid.iloc[:, 0].nunique() <= 1:
  599. return np.nan, np.nan
  600. y_true = valid.iloc[:, 1].astype(int).to_numpy()
  601. y_score = valid.iloc[:, 0].astype(float).to_numpy()
  602. return float(roc_auc_score(y_true, y_score)), float(average_precision_score(y_true, y_score))
  605. def _score_error_summary(table: pd.DataFrame, *, score_cols: list[str], error_col: str) -> pd.DataFrame:
  606. error = table[error_col].astype(float)
  607. high20 = (error >= error.quantile(0.80)).astype(int)
  608. high10 = (error >= error.quantile(0.90)).astype(int)
  609. rows = []
  610. for score_col in score_cols:
  611. score = table[score_col].astype(float)
  612. spearman, spearman_p = _safe_corr(score, error, method="spearman")
  613. pearson, pearson_p = _safe_corr(score, error, method="pearson")
  614. auc20, ap20 = _safe_auc(score, high20)
  615. auc10, ap10 = _safe_auc(score, high10)
  616. low_mask = score <= score.quantile(0.20)
  617. high_mask = score >= score.quantile(0.80)
  618. rows.append(
  619. {
  620. "score_name": score_col,
  621. "n_spots": int(score.notna().sum()),
  622. "error_col": error_col,
  623. "spearman_error": spearman,
  624. "spearman_pvalue": spearman_p,
  625. "pearson_error": pearson,
  626. "pearson_pvalue": pearson_p,
  627. "auroc_top20_error": auc20,
  628. "average_precision_top20_error": ap20,
  629. "auroc_top10_error": auc10,
  630. "average_precision_top10_error": ap10,
  631. "bottom20_score_mean_error": float(error.loc[low_mask].mean()),
  632. "top20_score_mean_error": float(error.loc[high_mask].mean()),
  633. "top_minus_bottom20_error": float(error.loc[high_mask].mean() - error.loc[low_mask].mean()),
  634. }
  635. )
  636. return pd.DataFrame(rows).sort_values(["auroc_top20_error", "spearman_error"], ascending=[False, False])
  639. def _selective_error_curve(
  640. table: pd.DataFrame,
  641. *,
  642. score_cols: list[str],
  643. error_col: str,
  644. keep_fractions: tuple[float, ...] = (0.5, 0.6, 0.7, 0.8, 0.9, 1.0),
  645. ) -> pd.DataFrame:
  646. full_mean = float(table[error_col].mean())
  647. high_error_threshold = float(table[error_col].quantile(0.8))
  648. rows = []
  649. for score_col in score_cols:
  650. ordered = table.sort_values(score_col, ascending=True)
  651. for keep_fraction in keep_fractions:
  652. n_keep = max(1, int(round(ordered.shape[0] * keep_fraction)))
  653. kept = ordered.head(n_keep)
  654. mean_error = float(kept[error_col].mean())
  655. rows.append(
  656. {
  657. "score_name": score_col,
  658. "keep_fraction": float(keep_fraction),
  659. "abstain_fraction": float(1.0 - keep_fraction),
  660. "n_kept": int(n_keep),
  661. "mean_error": mean_error,
  662. "median_error": float(kept[error_col].median()),
  663. "error_reduction_vs_full": float(1.0 - (mean_error / full_mean)) if full_mean > 0 else np.nan,
  664. "high_error_fraction": float((kept[error_col] >= high_error_threshold).mean()),
  665. "full_mean_error": full_mean,
  666. }
  667. )
  668. return pd.DataFrame(rows)
  671. def _scenario_summary(table: pd.DataFrame, *, score_col: str = "risk_score") -> pd.DataFrame:
  672. grouped = table.groupby("scenario", sort=True)
  673. return (
  674. grouped.agg(
  675. n_spots=("scenario", "size"),
  676. mean_true_error=("total_variation_error", "mean"),
  677. median_true_error=("total_variation_error", "median"),
  678. mean_risk_score=(score_col, "mean"),
  679. mean_phi_local=("phi_local", "mean"),
  680. mean_phi_uncertainty=("phi_uncertainty", "mean"),
  681. mean_phi_reference=("phi_reference", "mean"),
  682. )
  683. .reset_index()
  684. .sort_values("mean_true_error", ascending=False)
  685. )
  688. def _write_benchmark_h5ad_inputs(
  689. run_dir: Path,
  690. *,
  691. obs: pd.DataFrame,
  692. matrix: sparse.spmatrix | np.ndarray,
  693. genes: list[str],
  694. train_cells: dict[str, np.ndarray],
  695. celltypes: list[str],
  696. spot_table: pd.DataFrame,
  697. pseudo_counts: np.ndarray,
  698. label_col: str,
  699. layer: str,
  700. ) -> tuple[Path, Path]:
  701. train_indices = np.sort(np.concatenate([train_cells[celltype] for celltype in celltypes]))
  702. reference_obs = obs.iloc[train_indices].copy()
  703. reference_var = pd.DataFrame(index=pd.Index(genes, name=None).astype(str))
  704. reference_counts = matrix[train_indices, :]
  705. reference = ad.AnnData(X=reference_counts.copy(), obs=reference_obs, var=reference_var)
  706. reference.layers[layer] = reference_counts.copy()
  708. pseudo_obs = spot_table.copy()
  709. pseudo_var = pd.DataFrame(index=pd.Index(genes, name=None).astype(str))
  710. pseudo = ad.AnnData(X=pseudo_counts.copy(), obs=pseudo_obs, var=pseudo_var)
  711. pseudo.layers[layer] = pseudo_counts.copy()
  712. pseudo.obsm["spatial"] = pseudo_obs[["x_spatial", "y_spatial"]].to_numpy(dtype=float)
  713. pseudo.obs["sample_id"] = pseudo.obs["sample_id"].astype(str)
  715. reference_path = results_file(run_dir, "artifacts", "known_composition_train_reference.h5ad")
  716. pseudo_path = results_file(run_dir, "artifacts", "known_composition_pseudo_visium.h5ad")
  717. reference.write_h5ad(reference_path)
  718. pseudo.write_h5ad(pseudo_path)
  719. return reference_path, pseudo_path
  722. def _write_report(
  723. run_dir: Path,
  724. *,
  725. summary: pd.DataFrame,
  726. selective: pd.DataFrame,
  727. scenario: pd.DataFrame,
  728. metadata: dict[str, object],
  729. ) -> None:
  730. best = summary.iloc[0]
  731. risk_row = summary.loc[summary["score_name"] == "risk_score"]
  732. risk_text = "not available"
  733. if not risk_row.empty:
  734. row = risk_row.iloc[0]
  735. risk_text = (
  736. f"Spearman={row['spearman_error']:.3f}, "
  737. f"AUROC(top20 error)={row['auroc_top20_error']:.3f}, "
  738. f"top-bottom20 error gap={row['top_minus_bottom20_error']:.3f}"
  739. )
  741. keep80 = selective.loc[(selective["score_name"] == "risk_score") & (selective["keep_fraction"] == 0.8)]
  742. keep80_text = "not available"
  743. if not keep80.empty:
  744. row = keep80.iloc[0]
  745. keep80_text = (
  746. f"mean error={row['mean_error']:.3f}, "
  747. f"error reduction={row['error_reduction_vs_full']:.3f}, "
  748. f"high-error fraction={row['high_error_fraction']:.3f}"
  749. )
  751. lines = [
  752. "# Revision Known-Composition Benchmark",
  753. "",
  754. "## Purpose",
  755. "",
  756. "本运行生成 DLPFC pseudo-spots（伪空间点），保留 known cell-type composition（已知细胞类型组成），用于直接评估 ReliST risk score（ReliST 风险分数）与 true deconvolution error（真实反卷积误差）的关系。",
  757. "",
  758. "## Main Result Snapshot",
  759. "",
  760. f"- pseudo-spots（伪空间点）数量：`{metadata['n_spots']}`",
  761. f"- retained cell types（保留细胞类型）：`{metadata['n_celltypes']}`",
  762. f"- split mode（拆分方式）：`{metadata['split']['split_mode']}`，split column（拆分列）：`{metadata['split'].get('split_column', 'none')}`。",
  763. f"- primary `risk_score（风险分数）`：{risk_text}",
  764. f"- 低风险 keep 80%（保留 80% 低风险点）后：{keep80_text}",
  765. f"- 当前最佳 score（分数）：`{best['score_name']}`，AUROC(top20 error)={best['auroc_top20_error']:.3f}",
  766. "",
  767. "## Caveats",
  768. "",
  769. "- 这是 pseudo-spatial known-composition benchmark（伪空间已知组成基准），可直接回答审稿人关于 true error（真实误差）的核心问题，但仍不是自然组织中的 spot-level ground truth（空间点级真实标签）。",
  770. "- 当前默认不启用 base-model perturbation stability（基础模型扰动稳定性）；`phi_stability`（稳定性特征）保留为零列，避免偏离当前 manuscript boundary（手稿边界）。",
  771. "- `risk_score`（风险分数）在本运行中定义为 `phi_local（局部特征）`、`phi_uncertainty（输出模糊性）` 和 `phi_reference（参考特征）` 的等权标准化组合；高分表示更不可靠。",
  772. "",
  773. "## Output Tables",
  774. "",
  775. "- `tables/known_composition_spot_table.csv`",
  776. "- `tables/known_composition_true_abundance.csv`",
  777. "- `tables/known_composition_predicted_abundance.csv`",
  778. "- `tables/known_composition_risk_error_table.csv`",
  779. "- `tables/known_composition_celltype_inclusion.csv`",
  780. "- `tables/known_composition_score_error_summary.csv`",
  781. "- `tables/known_composition_selective_error_curve.csv`",
  782. "- `tables/known_composition_scenario_summary.csv`",
  783. "- `tables/reference_signatures_means.csv`",
  784. "- `tables/reference_signature_markers.csv`",
  785. "- `artifacts/known_composition_train_reference.h5ad`",
  786. "- `artifacts/known_composition_pseudo_visium.h5ad`",
  787. ]
  788. (run_dir / "revision_known_composition_benchmark.md").write_text("\n".join(lines) + "\n", encoding="utf-8")
  791. def main() -> int:
  792. args = parse_args()
  793. rng = np.random.default_rng(args.random_state)
  794. run_dir = args.output_root / "runs" / args.run_id
  795. ensure_results_layout(run_dir)
  796. set_selected_run(args.output_root, args.run_id)
  798. requested_genes = _read_gene_list(args.source_run_dir)
  799. source_celltype_order = _read_source_celltype_order(args.source_run_dir)
  800. obs, matrix, genes, celltypes, inclusion_table = _select_reference_matrix(
  801. args.reference_h5ad,
  802. layer=args.layer,
  803. label_col=args.reference_label_column,
  804. split_col=args.split_column if args.split_mode == "donor_disjoint" else None,
  805. requested_genes=requested_genes,
  806. source_celltype_order=source_celltype_order,
  807. min_cells_per_type=args.min_cells_per_type,
  808. )
  809. if args.split_mode == "donor_disjoint":
  810. train_cells, simulate_cells, celltypes, split_metadata = _split_reference_cells_donor_disjoint(
  811. obs,
  812. label_col=args.reference_label_column,
  813. split_col=args.split_column,
  814. celltypes=celltypes,
  815. train_fraction=args.train_fraction,
  816. min_cells_per_side=args.min_split_cells_per_type,
  817. rng=rng,
  818. )
  819. else:
  820. train_cells, simulate_cells = _split_reference_cells(
  821. obs,
  822. label_col=args.reference_label_column,
  823. celltypes=celltypes,
  824. train_fraction=args.train_fraction,
  825. rng=rng,
  826. )
  827. split_metadata = {
  828. "split_mode": "cell_random",
  829. "split_column": None,
  830. "train_fraction": float(args.train_fraction),
  831. "n_retained_celltypes_after_split": int(len(celltypes)),
  832. "dropped_after_split": [],
  833. }
  834. inclusion_table = _augment_inclusion_table(
  835. inclusion_table,
  836. celltypes_after_split=celltypes,
  837. train_cells=train_cells,
  838. simulate_cells=simulate_cells,
  839. )
  840. signatures = _compute_signatures(matrix, train_cells=train_cells, genes=genes, celltypes=celltypes)
  841. markers_all, marker_table = select_signature_markers(
  842. signatures,
  843. top_k=args.marker_top_k,
  844. min_positive_markers=args.min_positive_markers,
  845. )
  846. markers = subset_markers(markers_all, mode=args.marker_subset_mode)
  848. spot_table, true_abundance, region_prototypes, pseudo_counts = _sample_pseudo_spots(
  849. matrix,
  850. simulate_cells=simulate_cells,
  851. celltypes=celltypes,
  852. genes=genes,
  853. markers=markers_all,
  854. n_spots=args.n_spots,
  855. n_regions=args.n_regions,
  856. cells_per_spot_min=args.cells_per_spot_min,
  857. cells_per_spot_max=args.cells_per_spot_max,
  858. rng=rng,
  859. )
  860. reference_h5ad, pseudo_visium_h5ad = _write_benchmark_h5ad_inputs(
  861. run_dir,
  862. obs=obs,
  863. matrix=matrix,
  864. genes=genes,
  865. train_cells=train_cells,
  866. celltypes=celltypes,
  867. spot_table=spot_table,
  868. pseudo_counts=pseudo_counts,
  869. label_col=args.reference_label_column,
  870. layer=args.layer,
  871. )
  872. expression = _normalize_log_cp10k(pseudo_counts)
  873. expression_df = pd.DataFrame(expression, index=spot_table.index, columns=genes)
  875. predicted_values = _project_abundance(
  876. expression,
  877. signatures,
  878. method=args.projection_method,
  879. ridge_lambda=args.ridge_lambda,
  880. )
  881. predicted_abundance = pd.DataFrame(predicted_values, index=spot_table.index, columns=celltypes)
  882. error = _error_table(predicted_abundance, true_abundance)
  884. coords = spot_table[["x_spatial", "y_spatial"]].to_numpy(dtype=float)
  885. neighbors = knn_indices(coords, k=8)
  886. weights = inverse_distance_weights(coords, neighbors)
  887. heterogeneity = _compute_local_heterogeneity(expression, neighbors)
  888. stability_predictions = None
  889. if args.stability_repeats > 0:
  890. stability_predictions = gene_subsample_stability(
  891. expression,
  892. signatures.to_numpy(dtype=np.float32),
  893. repeats=args.stability_repeats,
  894. gene_fraction=args.stability_gene_fraction,
  895. ridge_lambda=args.ridge_lambda,
  896. random_state=args.random_state,
  897. )
  899. ambiguity = pd.Series(ambiguity_score(predicted_abundance), index=predicted_abundance.index, name="phi_uncertainty")
  900. model_output = BaseSpatialModelOutput(abundance=predicted_abundance, uncertainty=ambiguity)
  901. features = build_feature_table(
  902. model_output,
  903. neighbors=neighbors,
  904. weights=weights,
  905. heterogeneity=heterogeneity,
  906. stability_predictions=stability_predictions,
  907. confidence_proxy_precomputed=True,
  908. )
  910. marker_scores = compute_reference_marker_scores(expression_df, markers)
  911. features["phi_reference"] = reference_subsampling_instability(
  912. predicted_abundance,
  913. expression_df,
  914. markers,
  915. repeats=args.reference_repeats,
  916. subset_fraction=args.reference_fraction,
  917. random_state=args.random_state,
  918. )
  919. reference_marker = reference_marker_discordance_proxy(predicted_abundance, marker_scores)
  920. reference_residual = reference_signature_residual_proxy(
  921. predicted_abundance,
  922. expression_df,
  923. signatures,
  924. genes=marker_table["gene"].astype(str).tolist(),
  925. )
  927. groups = spot_table["sample_id"].astype(str)
  928. risk_table = features.copy()
  929. risk_table["risk_score"] = _combine_any_features(
  930. risk_table,
  931. {"phi_local": 1.0, "phi_uncertainty": 1.0, "phi_reference": 1.0},
  932. groups=groups,
  933. )
  934. risk_table["reference_risk_score"] = _combine_any_features(
  935. risk_table,
  936. {"phi_uncertainty": 2.0, "phi_reference": 2.0},
  937. groups=groups,
  938. )
  939. risk_table["local_uncertainty_risk_score"] = _combine_any_features(
  940. risk_table,
  941. {"phi_local": 1.0, "phi_uncertainty": 1.0},
  942. groups=groups,
  943. )
  944. for column, values in _abundance_baselines(predicted_abundance).items():
  945. risk_table[column] = values
  946. risk_table["snrna_marker_discordance"] = reference_marker
  947. risk_table["snrna_signature_residual"] = reference_residual
  948. risk_table = pd.concat([spot_table, risk_table, error], axis=1)
  950. score_cols = [
  951. "risk_score",
  952. "reference_risk_score",
  953. "local_uncertainty_risk_score",
  954. "abundance_entropy_risk",
  955. "inverse_top1_margin",
  956. "inverse_max_abundance",
  957. "phi_local",
  958. "phi_uncertainty",
  959. "phi_reference",
  960. "snrna_marker_discordance",
  961. "snrna_signature_residual",
  962. ]
  963. if args.stability_repeats > 0:
  964. score_cols.append("phi_stability")
  966. score_summary = _score_error_summary(risk_table, score_cols=score_cols, error_col="total_variation_error")
  967. selective_curve = _selective_error_curve(risk_table, score_cols=score_cols, error_col="total_variation_error")
  968. scenario_summary = _scenario_summary(risk_table)
  970. spot_table.to_csv(results_file(run_dir, "tables", "known_composition_spot_table.csv"))
  971. true_abundance.to_csv(results_file(run_dir, "tables", "known_composition_true_abundance.csv"))
  972. predicted_abundance.to_csv(results_file(run_dir, "tables", "known_composition_predicted_abundance.csv"))
  973. risk_table.to_csv(results_file(run_dir, "tables", "known_composition_risk_error_table.csv"))
  974. inclusion_table.to_csv(results_file(run_dir, "tables", "known_composition_celltype_inclusion.csv"), index=False)
  975. score_summary.to_csv(results_file(run_dir, "tables", "known_composition_score_error_summary.csv"), index=False)
  976. selective_curve.to_csv(results_file(run_dir, "tables", "known_composition_selective_error_curve.csv"), index=False)
  977. scenario_summary.to_csv(results_file(run_dir, "tables", "known_composition_scenario_summary.csv"), index=False)
  978. region_prototypes.to_csv(results_file(run_dir, "tables", "known_composition_region_prototypes.csv"))
  979. signatures.to_csv(results_file(run_dir, "tables", "reference_signatures_means.csv"))
  980. marker_table.to_csv(results_file(run_dir, "tables", "reference_signature_markers.csv"), index=False)
  981. expression_df.to_csv(results_file(run_dir, "tables", "known_composition_expression_log_cp10k.csv"))
  983. metadata = {
  984. "run_id": args.run_id,
  985. "reference_h5ad": str(args.reference_h5ad),
  986. "reference_label_column": args.reference_label_column,
  987. "source_run_dir": str(args.source_run_dir),
  988. "n_spots": int(args.n_spots),
  989. "n_genes": int(len(genes)),
  990. "n_celltypes": int(len(celltypes)),
  991. "n_original_labels": int(inclusion_table.shape[0]),
  992. "n_retained_before_split": int(inclusion_table["retained_before_split"].sum()),
  993. "n_retained_after_split": int(inclusion_table["retained_after_split"].sum()),
  994. "celltypes": celltypes,
  995. "random_state": int(args.random_state),
  996. "train_fraction": float(args.train_fraction),
  997. "split": split_metadata,
  998. "train_reference_h5ad": str(reference_h5ad),
  999. "pseudo_visium_h5ad": str(pseudo_visium_h5ad),
  1000. "cells_per_spot_min": int(args.cells_per_spot_min),
  1001. "cells_per_spot_max": int(args.cells_per_spot_max),
  1002. "pseudo_spot_scenarios": {
  1003. "clean": "No extra degradation beyond pseudo-bulk sampling.",
  1004. "low_depth": "Binomial thinning of counts with keep_probability sampled uniformly from 0.15 to 0.45.",
  1005. "marker_dropout": "Dominant-type marker counts thinned to 25% keep probability.",
  1006. "diffuse_mixture": "Composition smoothed toward a uniform mixture before cell sampling.",
  1007. },
  1008. "spatial_coordinate_generation": (
  1009. "Pseudo-spots are placed on a regular grid. Smooth region prototypes vary along the y coordinate; "
  1010. "therefore local-structure features should be interpreted with the shuffled-coordinate/null controls "
  1011. "added in downstream revision scripts."
  1012. ),
  1013. "risk_score_definition": "equal-weight grouped-zscore combination of phi_local, phi_uncertainty, and phi_reference",
  1014. "projection_method": args.projection_method,
  1015. "stability_repeats": int(args.stability_repeats),
  1016. "reference_feature_mode": "reference_subsampling_instability",
  1017. "reference_marker_subset_mode": args.marker_subset_mode,
  1018. "primary_error_col": "total_variation_error",
  1019. "score_columns": score_cols,
  1020. "manuscript_boundary": (
  1021. "This benchmark is a known-composition validation of risk-error association. "
  1022. "It does not make natural tissue spot-level truth claims."
  1023. ),
  1024. }
  1025. results_file(run_dir, "metadata", "known_composition_benchmark.json").write_text(
  1026. json.dumps(metadata, indent=2, ensure_ascii=False),
  1027. encoding="utf-8",
  1028. )
  1029. _write_report(run_dir, summary=score_summary, selective=selective_curve, scenario=scenario_summary, metadata=metadata)
  1031. print(f"Wrote known-composition benchmark to {run_dir}")
  1032. print(json.dumps({"run_id": args.run_id, "n_spots": args.n_spots, "n_celltypes": len(celltypes)}, indent=2))
  1033. return 0
  1036. if __name__ == "__main__":
  1037. raise SystemExit(main())

run_revision_known_composition_benchmark.py at commit ae4af50, under MIT · at the source

Overview

Authors: Xinyu Zhang1, Li He2, Yu Peng3, Yijia Li4, Jianyuan Kang1, Lisheng Peng4, Yifei Xu1, Sen Lin2
ORCID iDs: Xinyu Zhang
  1. Institute of Gastroenterology, Shenzhen Traditional Chinese Medicine Hospital, The Fourth Clinical Medical College of Guangzhou University of Chinese Medicine, Shenzhen, Guangdong, China
  2. Department of Oncology, Shenzhen Traditional Chinese Medicine Hospital, The Fourth Clinical Medical College of Guangzhou University of Chinese Medicine, Shenzhen, Guangdong, China
  3. Science and Technology Innovation Center, Guangzhou University of Chinese Medicine, Guangzhou, Guangdong, China
  4. Shenzhen Traditional Chinese Medicine Hospital, The Fourth Clinical Medical College of Guangzhou University of Chinese Medicine, Shenzhen, Guangdong, China
Institutions: Guangzhou University of Chinese Medicine (China); Shenzhen Traditional Chinese Medicine Hospital (China)
Journal: iScience, volume 29, issue 9, article 117206
Dates: received 21 May 2026; accepted 30 July 2026; published online 14 August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.isci.2026.117206 · PMID 42643167 · PMCID PMC13503129 · OpenAlex W7203496055
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), mouse (organism)
Methods: Connectivity, Statistics, Preprocessing
Keywords: ReliST, spatial transcriptomics, deconvolution, risk-aware deconvolution, reliability, selective prediction, risk-bin diagnostics, reference mismatch, artifact-aware interpretation
Topic: Single-cell and spatial transcriptomics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Sanming Project of Medicine in Shenzhen (SZZYSM202211003)
Citations: not cited yet (Europe PMC); 55 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

tp5353/ReliST

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: ae4af503c975f3bfad45ab89658ad92fd25c3465, 20 June 2026
Languages: Python (57), R (1)
Size: 107 files, 58 scripts
Software Heritage: not archived
Found in: “Data and code availability”
Holds: README, license file, environment (pyproject.toml), tests, documentation
Not found: CITATION.cff, continuous integration
Tools: pandas (37 files), NumPy (32 files), scikit-learn (8 files), SciPy (7 files), anndata (6 files), PyTorch (4 files), Matplotlib (1 file), Pillow (1 file), Scanpy (1 file), seaborn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
60 files
At the source: github.com/tp5353/ReliST

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;
  • 58 scripts, each with its path and the digest of its content;
  • 24 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

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it points to the authors' code: tp5353/ReliST
  • it says that the data are available on request
  • it says that the code is available on request

Read it in the paper: doi.org/10.1016/j.isci.2026.117206.

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

  • Authors: added Xinyu Zhang (0000-0001-7077-3314); removed Xinyu Zhang
  • Funding: added Sanming Project of Medicine in Shenzhen: SZZYSM202211003

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 9 keywords, 47 references.

Cite

This paper

Zhang, X., He, L., Peng, Y., Li, Y., Kang, J., Peng, L., Xu, Y., & Lin, S. (2026). ReliST: A model-agnostic risk layer for spatial transcriptomics deconvolution. iScience, 29(9), 117206. https://doi.org/10.1016/j.isci.2026.117206

BibTeX

@article{zhang2026relist,
author = {Zhang, Xinyu and He, Li and Peng, Yu and Li, Yijia and Kang, Jianyuan and Peng, Lisheng and Xu, Yifei and Lin, Sen},
title = {{ReliST: A model-agnostic risk layer for spatial transcriptomics deconvolution}},
journal = {iScience},
year = {2026},
month = aug,
volume = {29},
number = {9},
pages = {117206},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/j.isci.2026.117206},
url = {https://doi.org/10.1016/j.isci.2026.117206},
pmid = {42643167},
pmcid = {PMC13503129}
}

RIS

TY - JOUR
AU - Zhang, Xinyu
AU - He, Li
AU - Peng, Yu
AU - Li, Yijia
AU - Kang, Jianyuan
AU - Peng, Lisheng
AU - Xu, Yifei
AU - Lin, Sen
TI - ReliST: A model-agnostic risk layer for spatial transcriptomics deconvolution
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/08/14
VL - 29
IS - 9
SP - 117206
SN - 2589-0042
PB - Elsevier
DO - 10.1016/j.isci.2026.117206
UR - https://doi.org/10.1016/j.isci.2026.117206
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.isci.2026.117206",
"type": "article-journal",
"title": "ReliST: A model-agnostic risk layer for spatial transcriptomics deconvolution",
"container-title": "iScience",
"author": [
{
"family": "Zhang",
"given": "Xinyu"
},
{
"family": "He",
"given": "Li"
},
{
"family": "Peng",
"given": "Yu"
},
{
"family": "Li",
"given": "Yijia"
},
{
"family": "Kang",
"given": "Jianyuan"
},
{
"family": "Peng",
"given": "Lisheng"
},
{
"family": "Xu",
"given": "Yifei"
},
{
"family": "Lin",
"given": "Sen"
}
],
"container-title-short": "iScience",
"volume": "29",
"issue": "9",
"page": "117206",
"DOI": "10.1016/j.isci.2026.117206",
"PMID": "42643167",
"PMCID": "PMC13503129",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.isci.2026.117206",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
14
]
]
}
}

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.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: anndata, Scanpy, Pillow, 7 other tools, genetics / omics, 14 references
[2] doi:10.1093/bioinformatics/btag578 [code]
NicheDeSig: niche-aware deconvolution and adaptive signature analysis for spatial transcriptomics.
Journal: Bioinformatics (Oxford, England)
In common: anndata, Scanpy, PyTorch, 4 other tools, genetics / omics, 12 references
[3] 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: anndata, Scanpy, PyTorch, 6 other tools, genetics / omics, 8 references
[4] doi:10.1093/bib/bbag404 [code]
Navigating cell maps by deep learning integration of single-cell and spatially resolved transcriptomics.
Journal: Briefings in bioinformatics
In common: anndata, Scanpy, PyTorch, 5 other tools, genetics / omics, mouse, 8 references
[5] doi:10.1093/bioinformatics/btag515 [code]
PRISM: Prior-enhanced Inference for Spatial Transcriptomic Cell Type Mapping.
Journal: Bioinformatics (Oxford, England)
In common: anndata, Scanpy, PyTorch, 5 other tools, genetics / omics, 8 references
[6] doi:10.1093/bib/bbag298 [code]
Empowering multifaceted analysis of spatial transcriptomics data with RGAST.
Journal: Briefings in bioinformatics
In common: anndata, Scanpy, PyTorch, 6 other tools, genetics / omics, mouse, 7 references
[7] doi:10.1093/bioinformatics/btag424 [code]
SlotDeconv: spatial transcriptomics deconvolution via diversity-constrained prototype learning and spatial refinement.
Journal: Bioinformatics (Oxford, England)
In common: anndata, Scanpy, PyTorch, 6 other tools, genetics / omics, mouse, 6 references
[8] doi:10.1038/s41592-026-03211-w [code]
Spatial isoform sequencing at single-cell resolution reveals cell-type-specific spatial isoform variability in multiple brain cell types.
Journal: Nature methods
In common: Scanpy, Pillow, seaborn, 5 other tools, genetics / omics, mouse, 7 references
[9] doi:10.1038/s41592-026-03194-8 [code]
Beyond benchmarking: an expert-guided consensus approach to spatially aware clustering.
Journal: Nature methods
In common: anndata, Scanpy, Pillow, 7 other tools, genetics / omics, 5 references
[10] doi:10.1038/s41593-026-02293-1 [code]
Optics-free spatial genomics for mapping mammalian brain aging by IRISeq.
Journal: Nature neuroscience
In common: Scanpy, seaborn, scikit-learn, 4 other tools, genetics / omics, mouse, 8 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.