OSCR

Scalable discovery of spatial multicellular patterns via neighborhood-to-sequence transformation.

Code ↔ Paper

1 match 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 1 match
  1. [1] § Methods › Spatial annotation of hepatic cell types and vascular structures in the MERFISH Liver Atlas ↔ code/simu_utils.R, lines 302–329 · score 0.54 · cell populations, spatial transcriptomic, inferred, gene, tissue

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

R · 677 lines · 21 KB · no license · 1 match

  1. # Author: Zheng Li
  2. # Date: 2021-05-19
  3. # Purpose:
  4. # This script includes all functions to conduct the simulation study for BASS
  6. suppressPackageStartupMessages({
  7. library(BASS)
  8. library(Giotto)
  9. library(BayesSpace)
  10. library(Seurat)
  11. library(SC3)
  12. library(GapClust)
  13. library(mclust)
  14. library(tidyverse)
  15. library(scater)
  16. library(gtools)
  17. library(splatter)
  18. library(reticulate)
  19. library(CVXR)
  20. })
  23. #' Concatenate two strings
  24. "%&%" <- function(x, y) paste0(x, y)
  27. #' Cluster cells with Seurat
  28. seu_cluster <- function(sim_dat, C, resolutions = seq(0.1, 4, by = 0.1))
  29. {
  30. cnts <- do.call(cbind, sim_dat[[1]])
  31. colnames(cnts) <- "Cell" %&% 1:ncol(cnts)
  32. seu <- CreateSeuratObject(counts = cnts, min.cells = 1, min.features = 1)
  33. seu <- NormalizeData(seu, verbose = F)
  34. seu <- ScaleData(seu, features = rownames(seu), verbose = F)
  35. seu <- RunPCA(seu, features = rownames(seu), verbose = F)
  36. seu <- FindNeighbors(seu, dims = 1:20, verbose = F)
  37. seu <- FindClusters(seu, resolution = resolutions, verbose = F)
  38. kest <- sapply("RNA_snn_res." %&% resolutions, function(res){
  39. length(unique(seu[[res, drop = T]]))
  40. })
  41. out_res <- names(which.min(abs(kest - C)))
  42. c_est <- as.numeric(as.character(seu[[out_res, drop = T]])) + 1
  43. # evaluation
  44. ctrue <- unlist(sim_dat[[2]])
  45. ari <- eval_ARI(c_est, ctrue)
  46. F1 <- eval_F1(c_est, ctrue)
  47. MCC <- eval_MCC(c_est, ctrue)
  48. C_est <- length(unique(c_est))
  49. clust_props <- calc_clust_prop(c_est, ctrue)
  51. list(metric = c(ari, F1, MCC), C_est = C_est, clust_props = clust_props)
  52. }
  55. #' Cluster cells with SC3
  56. sc3_cluster <- function(sim_dat, C, seed = 0)
  57. {
  58. set.seed(seed)
  59. cnts <- do.call(cbind, sim_dat[[1]])
  60. x <- scater::normalizeCounts(cnts, log = TRUE)
  61. sce <- SingleCellExperiment(assays = list(counts = cnts, logcounts = x))
  62. rowData(sce)$feature_symbol <- rownames(sce)
  63. sce <- tryCatch({
  64. sc3(sce, ks = C, n_cores = 1)},
  65. error = function(e){
  66. print(e)
  67. return("error")}
  68. )
  69. if(is.character(sce)){
  70. out <- list(metric = rep("error", 9), C_est = "error",
  71. clust_props = "error")
  72. } else{
  73. if(ncol(x) > 5000) sce <- sc3_run_svm(sce, ks = C)
  74. c_est <- colData(sce)[, "sc3_" %&% C %&% "_clusters"]
  75. ctrue <- unlist(sim_dat[[2]])
  76. ari <- eval_ARI(c_est, ctrue)
  77. F1 <- eval_F1(c_est, ctrue)
  78. MCC <- eval_MCC(c_est, ctrue)
  79. C_est <- length(unique(c_est))
  80. clust_props <- calc_clust_prop(c_est, ctrue)
  81. out <- list(metric = c(ari, F1, MCC), C_est = C_est,
  82. clust_props = clust_props)
  83. }
  84. return(out)
  85. }
  88. #' Cluster cells with FICT
  89. #' FICT input files:
  90. #' 1.expression, NxJ
  91. #' 2.coordinates,
  92. fict_cluster <- function(sim_dat, xy, C, case, rep)
  93. {
  94. # The following environmental variables need to be set
  95. # to run FICT successfully
  96. PYTHONPATH <- paste(
  97. "/net/mulan/home/zlisph/BASS_pjt/9_FICT/FICT-SAMPLE/FICT/",
  98. "/net/mulan/home/zlisph/BASS_pjt/9_FICT/FICT-SAMPLE/GECT/", sep = ":")
  99. Sys.setenv(PYTHONPATH = PYTHONPATH)
  100. Sys.setenv(MPLBACKEND = 'Agg')
  101. fict <- "~/softwares/miniconda3/bin/fict"
  102. tmp_fder <- "~/BASS_pjt/2_simu/fict_tmp/case" %&% case %&% "/rep" %&% rep
  103. if(!file.exists(tmp_fder)) dir.create(tmp_fder, recursive = T)
  104. # prepare input
  105. # 1.coordinates
  106. write.table(xy, file = tmp_fder %&% "/simu.coordinates",
  107. col.names = F, quote = F)
  108. # 2.expression
  109. cnts <- do.call(cbind, sim_dat[[1]])
  110. ge <- scater::normalizeCounts(cnts, log = TRUE)
  111. ge <- t(ge) # ncell x ngene
  112. write.table(ge, file = tmp_fder %&% "/simu.expression",
  113. col.names = F, quote = F)
  115. # run FICT
  116. # Note:
  117. # by default hidden = 20, size of denoise auto-encoder
  118. cmd <- paste(fict, "-p", tmp_fder %&% "/simu", "-o",
  119. tmp_fder %&% "/out/", "--n_type", C)
  120. system(cmd, ignore.stdout = T, ignore.stderr = T)
  122. # evaluate
  123. # occasionally have error: numpy.linalg.LinAlgError: Internal Error.
  124. c_est <- tryCatch({
  125. read.table(tmp_fder %&% "/out/cluster_result.csv")
  126. }, error = function(e) e)
  127. if(!is.data.frame(c_est)){
  128. out <- list(metric = rep("error", 9), C_est = "error",
  129. clust_props = "error")
  130. } else{
  131. c_est <- c_est$V1+1
  132. ctrue <- unlist(sim_dat[[2]])
  133. ari <- eval_ARI(c_est, ctrue)
  134. F1 <- eval_F1(c_est, ctrue)
  135. MCC <- eval_MCC(c_est, ctrue)
  136. C_est <- length(unique(c_est))
  137. clust_props <- calc_clust_prop(c_est, ctrue)
  138. out <- list(metric = c(ari, F1, MCC), C_est = C_est,
  139. clust_props = clust_props)
  140. }
  141. # remove tmp files
  142. unlink(tmp_fder, recursive = T)
  143. return(out)
  144. }
  147. #' Match spatial domain labels and cell type labels
  148. #' in the estimated proportion matrix (pi_est) with
  149. #' labels in the true proportion matrix
  150. match_pi <- function(pi_est, pi_true, z_est, z_true)
  151. {
  152. C <- nrow(pi_true)
  153. R <- ncol(pi_true)
  154. # 1. Find spatial domain label correspondence
  155. perms_R <- permutations(n = R, r = R)
  156. accur <- rep(NA, nrow(perms_R))
  157. for(i in 1:nrow(perms_R))
  158. {
  159. perm <- perms_R[i, ]
  160. z_true_perm <- case_when(
  161. z_true == 1 ~ perm[1],
  162. z_true == 2 ~ perm[2],
  163. z_true == 3 ~ perm[3],
  164. z_true == 4 ~ perm[4])
  165. accur[i] <- mean(z_est == z_true_perm)
  166. }
  167. perm_final_R <- perms_R[which.max(accur), ]
  169. # 2. Match spatial domain label in pi_est
  170. pi_est_match <- pi_est[, perm_final_R]
  172. # 3. Match cell type label
  173. perms_C <- permutations(n = C, r = C)
  174. sse <- rep(NA, nrow(perms_C))
  175. for(i in 1:nrow(perms_C))
  176. {
  177. perm <- perms_C[i, ]
  178. sse[i] <- sum((pi_est_match[perm, ] - pi_true)^2)
  179. }
  180. perm_final_C <- perms_C[which.min(sse), ]
  181. pi_est_match <- pi_est_match[perm_final_C, ]
  182. out <- list(
  183. pi = pi_est_match,
  184. perm_K = perm_final_R,
  185. perm_C = perm_final_C)
  186. }
  189. #' Map spatial domain label to cell types in that domain
  190. #' 1 -> (1, 2, 3), 2 -> (2, 3, 4)
  191. #' 3 -> (3, 4, 1), 4 -> (4, 1, 2)
  192. map_z2c <- function(z)
  193. {
  194. case_when(
  195. z == 1 ~ c(1, 2, 3),
  196. z == 2 ~ c(2, 3, 4),
  197. z == 3 ~ c(3, 4, 1),
  198. z == 4 ~ c(4, 1, 2)
  199. )
  200. }
  203. #' Produce the confusion matrix
  204. #' One difficulty for constructing the confusion matrix is that
  205. #' the labeling from the clustering and truth may not match.
  206. #' We get around this difficulty by identifying the confusion matrix
  207. #' such that the number of true positives (sum of diagonal values)
  208. #' is maximized.
  209. #' Note the number of clusters in est_labels needs to be smaller than
  210. #' or equal to the number of clusters in true_labels
  211. confusion <- function(est_labels, true_labels)
  212. {
  213. C <- length(unique(true_labels))
  214. est_labels <- factor(est_labels, 1:C)
  215. A <- matrix(table(est_labels, true_labels), C, C)
  216. P <- Variable(C, C, boolean = T) # row permutation matrix
  217. Y <- Variable(C, C)
  219. problem <- Problem(Maximize(matrix_trace(Y)), list(Y == P %*% A,
  220. sum_entries(P, axis = 1) == 1, sum_entries(P, axis = 2) == 1))
  221. result <- solve(problem)
  222. P <- result$getValue(P)
  223. confusion_mtx <- P %*% A
  224. }
  227. #' Evaluate the F1-score
  228. #' Refer to paper Chicco and Jurman, BMC Genomics, 2020
  229. eval_F1 <- function(est_labels, true_labels)
  230. {
  231. C <- length(unique(true_labels))
  232. confusion_mtx <- suppressMessages(confusion(est_labels, true_labels))
  233. TP <- diag(confusion_mtx)
  234. FP <- apply(confusion_mtx, 1, sum) - diag(confusion_mtx)
  235. FN <- apply(confusion_mtx, 2, sum) - diag(confusion_mtx)
  236. F1 <- 2 * TP / (2 * TP + FP + FN)
  237. c(all_F1 = mean(F1), major_F1 = mean(F1[1:4]), rare_F1 = mean(F1[5:C]))
  238. }
  241. #' Evaluate Matthew correlation coefficient (MCC) score
  242. #' Refer to paper Chicco and Jurman, BMC Genomics, 2020
  243. eval_MCC <- function(est_labels, true_labels)
  244. {
  245. C <- length(unique(true_labels))
  246. confusion_mtx <- suppressMessages(confusion(est_labels, true_labels))
  247. TP <- diag(confusion_mtx)
  248. FP <- apply(confusion_mtx, 1, sum) - diag(confusion_mtx)
  249. FN <- apply(confusion_mtx, 2, sum) - diag(confusion_mtx)
  250. TN <- sum(confusion_mtx) - TP - FP - FN
  251. MCC <- (TP * TN - FP * FN) /
  252. sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN))
  253. MCC[is.na(MCC)] <- 0
  254. c(all_MCC = mean(MCC), major_MCC = mean(MCC[1:4]), rare_MCC = mean(MCC[5:C]))
  255. }
  258. #' Evaluate the adjusted Random Index(ARI)
  259. eval_ARI <- function(est_labels, true_labels)
  260. {
  261. major_cells <- true_labels %in% 1:4
  262. rare_cells <- !major_cells
  263. all_ari <- adjustedRandIndex(est_labels, true_labels)
  264. major_ari <- adjustedRandIndex(
  265. est_labels[major_cells], true_labels[major_cells])
  266. rare_ari <- adjustedRandIndex(
  267. est_labels[rare_cells], true_labels[rare_cells])
  268. c(all_ari = all_ari, major_ari = major_ari, rare_ari = rare_ari)
  269. }
  272. #' Reorder the labeling of cell type clusters/spatial domains to
  273. #' cell types/spatial domains 1, 2, 3, 4 and then order the remaining
  274. #' redundant clusters/domains based on their cluster/domain size from
  275. #' large to small. Finally, calculate the proportion. Note this quantity
  276. #' is specifically used for evaluation in the case of a mis-specified
  277. #' number of cell types/spatial domains.
  278. calc_clust_prop <- function(est_labels, true_labels)
  279. {
  280. C <- length(unique(true_labels))
  281. C_est <- length(unique(est_labels))
  282. C_min <- min(C, C_est)
  283. C_max <- max(C, C_est)
  284. mtx <- matrix(table(est_labels, true_labels), C_est, C)
  285. P <- Variable(C_est, C_est, boolean = T) # row permutation matrix
  286. Y <- Variable(C_est, C)
  287. problem <- suppressMessages(Problem(
  288. Maximize(matrix_trace(Y[1:C_min, 1:C_min])),
  289. list(Y == P %*% mtx, sum_entries(P, axis = 1) == 1,
  290. sum_entries(P, axis = 2) == 1)))
  291. result <- suppressMessages(solve(problem))
  292. P <- result$getValue(P)
  293. mtx <- P %*% mtx
  294. order_idx <- order(apply(mtx, 1, sum)[-c(1:C_min)], decreasing = T)
  295. mtx <- mtx[c(1:C_min, order_idx+C_min), , drop = F]
  296. props <- apply(mtx, 1, sum) / sum(mtx)
  297. names(props) <- 1:C_est
  298. paste(props, collapse = ",")
  299. }
  302. #' Generate simulated spatial transcriptomic data with splatter package
  303. #'
  304. #' Spatial domain labels are based on STARmap (20180417_BZ5_control) data
  305. #' and annotated based on the marker genes. Spatial domain labels for
  306. #' section 1 to section 9 are generated based on the annotated labels.
  307. #'
  308. #' @param starmap STARmap data for inferring simulation parameters and for
  309. #' providing realistic cell coordinates and spatial domains.
  310. #' @param scenario Simulation scenario (refer to the manuscript for details)
  311. #' @param rare_dist Rare cell types either randomly ditributed across the entire
  312. #' tissue section or located in certain spatial domians.
  313. #' @param C Number of cell types that consist of 4 major cell types and (C-4)
  314. #' rare cell types. Rare cell types together consist of 30% of the total cell
  315. #' population and randomly distribued across the entire tissue region.
  316. #' Each rare cell type consist of 30%/(C-4) of the total cell population.
  317. #' @param J Number of genes.
  318. #' @param I Number of randomly selected genes among the J genes to remain in
  319. #' the final dataset.
  320. #' @param L Number of tissue sections.
  321. #' @param batch_facLoc Batch factor location (refer to the splatter package).
  322. #' Zero if there is no batch effect.
  323. #' @param de_prop Probability that a gene will be selected to be differentially
  324. #' expressed for each cell type (refer to splatter package).
  325. #' @param de_facLoc DE factor location (refer to splatter package).
  326. #' @param de_facScale DE factor scale (refer to splatter package).
  327. #' @param sim_seed Random seed.
  328. #' @param debug Output DE genes for each cell type.
  329. simu <- function(
  330. starmap,
  331. scenario,
  332. rare_dist = c("random", "spatial"),
  333. C,
  334. J,
  335. I = NULL,
  336. L,
  337. batch_facLoc,
  338. de_prop,
  339. de_facLoc,
  340. de_facScale,
  341. sim_seed,
  342. debug = FALSE
  343. )
  344. {
  345. cnts <- starmap$cnts
  346. info <- starmap$info
  347. N <- nrow(info)
  348. init_params <- splatEstimate(cnts[, which(info$c == "eL2/3")])
  350. # reproduce randomness in cell type assignment and
  351. # count data selection.
  352. set.seed(sim_seed)
  354. # 1.simulate count data
  355. noBatch <- ifelse(batch_facLoc == 0, TRUE, FALSE)
  356. group_prob <- if(C == 4){
  357. rep(0.25, 4)
  358. } else if(C > 4){
  359. c(rep(0.7 / 4, 4), rep(1 / (C - 4), C - 4) * 0.3)
  360. }
  361. params <- setParams(
  362. init_params,
  363. batchCells = rep(3 * N, L), # 3N here represents a large number such that
  364. # we have sufficient cells of each type to be
  365. # allocated to the spatial transcriptomics data
  366. batch.rmEffect = noBatch,
  367. batch.facLoc = batch_facLoc,
  368. nGenes = J,
  369. group.prob = group_prob,
  370. out.prob = 0,
  371. de.prob = de_prop,
  372. de.facLoc = de_facLoc,
  373. de.facScale = de_facScale,
  374. seed = sim_seed)
  375. sim_groups <- splatSimulate(
  376. params = params,
  377. method = "groups",
  378. verbose = FALSE)
  379. if(!is.null(I)){
  380. # down-sample J genes to I genes
  381. keep_idx <- sort(sample(1:J, I, replace = F))
  382. sim_groups <- sim_groups[keep_idx, ]
  383. } else{
  384. I <- J
  385. }
  386. # remove cells having no expressed genes
  387. idx_zerosize <- apply(counts(sim_groups), MARGIN = 2, sum) == 0
  388. sim_groups <- sim_groups[, !idx_zerosize]
  390. # 2.parse proportion of cells types in each spatial domain
  391. if(scenario == 2){
  392. prop <- c(0.8, 0.1, 0.1)
  393. } else if(scenario == 3){
  394. prop <- c(0.5, 0.25, 0.25)
  395. } else if(scenario == 4){
  396. prop <- rep(1, 3) / 3
  397. }
  399. c_simu <- list()
  400. sim_cnt <- list()
  401. for(l in 1:L)
  402. {
  403. # 3.set spatial domain labels
  404. if(l == 1){
  405. ztrue <- info$z
  406. } else{
  407. ztrue <- info[["slice" %&% (l-1)]]
  408. }
  410. # 4.generate cell types
  411. c_simu[[l]] <- rep(NA, length(ztrue))
  412. if(scenario == 1){
  413. c_simu[[l]] <- ztrue
  414. } else if(scenario %in% c(2, 3, 4)){
  415. for(z in unique(info$z))
  416. {
  417. zi_idx <- ztrue == z
  418. c_simu[[l]][zi_idx] <- sample(map_z2c(z), sum(zi_idx), prob = prop,
  419. replace = T)
  420. }
  421. } else if(scenario == 5){ # more challenging scenario with rare cell types
  422. # assign rare cell types
  423. if(rare_dist[1] == "spatial"){
  424. rare_idx <- rep(NA, length(ztrue))
  425. nrare_each <- round(0.3 / (C - 4) * length(ztrue))
  426. for(c in 5:C)
  427. {
  428. repeat{
  429. zi_idx <- sample(unique(info$z), 1) == ztrue
  430. zi_noAlloc_idx <- zi_idx & is.na(rare_idx)
  431. if(sum(zi_noAlloc_idx) >= nrare_each) break
  432. }
  433. ci_idx <- sample(which(zi_noAlloc_idx), nrare_each, replace = F)
  434. c_simu[[l]][ci_idx] <- c
  435. rare_idx[ci_idx] <- TRUE
  436. }
  437. rare_idx[is.na(rare_idx)] <- FALSE
  438. } else if(rare_dist[1] == "random"){
  439. rare_idx <- sample(c(TRUE, FALSE), length(ztrue), prob = c(0.3, 0.7),
  440. replace = T)
  441. c_simu[[l]][rare_idx] <- sample(5:C, sum(rare_idx), replace = T)
  442. }
  443. major_idx <- !rare_idx
  445. # assign major cell types using proportion of
  446. # cell types in scenario 3
  447. for(z in unique(info$z))
  448. {
  449. zi_idx <- ztrue == z
  450. zi_major_idx <- zi_idx & major_idx
  451. c_simu[[l]][zi_major_idx] <- sample(map_z2c(z), sum(zi_major_idx),
  452. prob = c(0.5, 0.25, 0.25), replace = T)
  453. }
  454. }
  456. # 5.assign count data
  457. groups <- as.data.frame(colData(sim_groups)) %>%
  458. filter(Batch == "Batch" %&% l)
  459. sim_cnt[[l]] <- array(NA, c(I, N))
  460. for(c in 1:C)
  461. {
  462. c_size <- sum(c_simu[[l]] == c)
  463. c_cells <- groups$Cell[grepl("Group" %&% c %&% "$", groups$Group)]
  464. cells_select <- sample(as.character(c_cells), c_size, replace = F)
  465. sim_cnt[[l]][, c_simu[[l]] == c] <-
  466. as.matrix(counts(sim_groups)[, cells_select])
  467. }
  468. colnames(sim_cnt[[l]]) <- "Cell" %&% 1:N
  469. rownames(sim_cnt[[l]]) <- "Gene" %&% 1:I
  470. }
  472. # return DE genes for each group
  473. if(debug){
  474. de_genes <- list()
  475. for(c in 1:C)
  476. {
  477. de_genes[[c]] <- rownames(rowData(sim_groups))[
  478. rowData(sim_groups)[, "DEFacGroup" %&% c] != 1]
  479. }
  480. return(list(sim_cnt, c_simu, sim_seed, de_genes))
  481. }
  483. return(list(sim_cnt, c_simu, sim_seed))
  484. }
  487. #' Run BASS algorithm
  488. run_BASS <- function(sim_dat, xy, beta_method, beta, C, R,
  489. init_method)
  490. {
  491. L <- length(sim_dat[[1]])
  492. xys <- lapply(1:L, function(x) xy)
  493. sim_cnt <- sim_dat[[1]]
  494. for(l in 1:L)
  495. {
  496. colnames(sim_cnt[[l]]) <- rownames(xys[[l]])
  497. }
  499. # run algorithms
  500. BASS <- createBASSObject(sim_cnt, xys, C = C, R = R,
  501. init_method = init_method, beta_method = "SW",
  502. beta = beta, tol = 1e-4, burnin = 2000, nsample = 2000)
  503. BASS <- BASS.preprocess(BASS)
  504. BASS <- BASS.run(BASS)
  505. BASS <- BASS.postprocess(BASS)
  507. ctrue <- unlist(sim_dat[[2]])
  508. ztrue <- unlist(info[, 4:(L+3)])
  509. pi_true <- table(ctrue, ztrue)
  510. pi_true <- pi_true %*% diag(1 / apply(pi_true, 2, sum))
  511. c_est <- unlist(BASS@results$c)
  512. z_est <- unlist(BASS@results$z)
  513. pi_est <- BASS@results$pi
  515. # evaluation
  516. c_ari <- eval_ARI(c_est, ctrue)
  517. F1 <- eval_F1(c_est, ctrue)
  518. MCC <- eval_MCC(c_est, ctrue)
  519. z_ari <- adjustedRandIndex(ztrue, z_est)
  520. if(C != 4 | R != 4){ # only evaluate in the main simulation
  521. pi_est <- NA
  522. mse_pi <- NA
  523. } else{
  524. pi_est <- match_pi(pi_est, pi_true, z_est, ztrue)$pi
  525. mse_pi <- sqrt(sum((pi_est - pi_true)^2))
  526. }
  527. C_est <- length(unique(c_est))
  528. R_est <- length(unique(z_est))
  529. c_clust_prop <- calc_clust_prop(c_est, ctrue)
  530. z_clust_prop <- calc_clust_prop(z_est, ztrue)
  532. output <- list(
  533. c_ari = c_ari,
  534. c_F1 = F1,
  535. c_MCC = MCC,
  536. C_est = C_est,
  537. c_clust_prop = c_clust_prop,
  538. z_ari = z_ari,
  539. R_est = R_est,
  540. z_clust_prop = z_clust_prop,
  541. pi_est = pi_est,
  542. mse_pi = mse_pi,
  543. beta = BASS@results$beta)
  544. return(output)
  545. }
  548. #' Run HMRF algorithm
  549. run_HMRF <- function(sim_dat, xy, ztrue, R, case, rep,
  550. usePCs = F, SEgenes = NULL, dosearchSEgenes = T)
  551. {
  552. sim_cnt <- sim_dat[[1]][[1]]
  553. J <- nrow(sim_cnt)
  555. # 0.prepare
  556. my_python_path <- "/net/mulan/home/zlisph//softwares/miniconda3/bin/python"
  557. results_folder <- "/net/mulan/home/zlisph/BASS_pjt/2_simu/HMRF_tmp/case" %&%
  558. case %&% "/rep" %&% rep
  559. if(!file.exists(results_folder)) dir.create(results_folder, recursive = T)
  560. instrs <- createGiottoInstructions(python_path = my_python_path)
  562. # 1.run HMRF
  563. starmap <- createGiottoObject(
  564. raw_exprs = as.data.frame(sim_cnt),
  565. spatial_locs = xy,
  566. instructions = instrs)
  567. starmap <- filterGiotto(
  568. gobject = starmap,
  569. expression_threshold = 0.5,
  570. gene_det_in_min_cells = 20,
  571. min_det_genes_per_cell = 0)
  572. starmap <- normalizeGiotto(starmap)
  573. starmap <- addStatistics(starmap)
  574. starmap <- createSpatialNetwork(starmap, minimum_k = 2)
  576. # select SE genes using
  577. if(dosearchSEgenes){
  578. genes <- binSpect(starmap, bin_method = 'kmeans')$genes
  579. genes <- if(length(genes) < 100) genes else genes[1:100]
  580. } else{
  581. genes <- rownames((starmap@norm_scaled_expr))
  582. }
  584. if(usePCs){
  585. starmap <- Giotto::runPCA(
  586. gobject = starmap,
  587. genes_to_use = genes,
  588. ncp = 20)
  589. betas <- c(0, 1, 21)
  590. output_folder <- results_folder %&% "/result_k" %&% R %&% "_pca"
  591. HMRF_out <- doHMRF(
  592. gobject = starmap,
  593. dim_reduction_to_use = "pca",
  594. dimensions_to_use = 1:20,
  595. k = R,
  596. betas = betas,
  597. output_folder = output_folder,
  598. overwrite_output = T)
  599. } else{
  600. betas <- c(0, 2, 26)
  601. output_folder <- results_folder %&% "/result_k" %&% R
  602. HMRF_out <- doHMRF(
  603. gobject = starmap,
  604. expression_values = "scaled",
  605. spatial_genes = genes,
  606. k = R,
  607. betas = betas,
  608. output_folder = output_folder,
  609. overwrite_output = T)
  610. }
  612. # 2.evaluate results
  613. res_dir <- output_folder %&% "/result.spatial.zscore/k_" %&% R
  614. allbetas <- seq(betas[1], by = betas[2], length = betas[3])
  615. ari <- list()
  616. R_est <- list()
  617. z_clust_prop <- list()
  618. for(beta in allbetas)
  619. {
  620. beta <- beta %&% ".0"
  621. probs <- tryCatch(
  622. read.table(res_dir %&% "/ftest.beta." %&% beta %&% ".unnormprob.txt",
  623. row.names = 1),
  624. error = function(e) "error",
  625. warning = function(w) "warning")
  627. if(!is.data.frame(probs)){
  628. out <- list(ari = rep("error", length(allbetas)),
  629. R_est = rep("error", length(allbetas)),
  630. z_clust_prop = rep("error", length(allbetas)))
  631. return(out)
  632. } else{
  633. z_est <- apply(probs, MARGIN = 1, which.max)
  634. ari[[beta]] <- adjustedRandIndex(ztrue, z_est)
  635. R_est[[beta]] <- length(unique(z_est))
  636. z_clust_prop[[beta]] <- calc_clust_prop(z_est, ztrue)
  637. }
  638. }
  639. unlink(results_folder, recursive = T)
  640. out <- list(ari = ari, R_est = R_est, z_clust_prop = z_clust_prop)
  641. return(out)
  642. }
  645. #' Run BayesSpace
  646. run_BayesSpace <- function(sim_dat, xy, ztrue, R, seed = 0)
  647. {
  648. set.seed(seed)
  649. sim_cnt <- sim_dat[[1]][[1]]
  650. colnames(xy) <- c("row", "col")
  651. info <- data.frame(xy)
  652. sce <- SingleCellExperiment(assays = list(counts = sim_cnt), colData = info)
  653. sce <- spatialPreprocess(sce, n.PCs = 15, n.HVGs = nrow(sce), log.normalize = T)
  654. sce <- spatialCluster(sce, q = R, d = 15,
  655. init.method = "mclust", model = "t",
  656. nrep = 10000, burn.in = 1000)
  657. z_est <- colData(sce)$spatial.cluster
  658. ari <- adjustedRandIndex(ztrue, z_est)
  659. R_est <- length(unique(z_est))
  660. z_clust_prop <- calc_clust_prop(z_est, ztrue)
  661. out <- list(ari = ari, R_est = R_est, z_clust_prop = z_clust_prop)
  662. return(out)
  663. }
  666. #' run_SpGCN
  667. run_SpaGCN <- function(sim_dat, xy, ztrue, R, p = 0.5)
  668. {
  669. sim_cnt <- sim_dat[[1]][[1]]
  670. xy <- data.frame(xy)
  671. z_est <- run_SpaGCN_py(t(sim_cnt), xy, R, p)$refined_pred
  672. ari <- adjustedRandIndex(ztrue, z_est)
  673. R_est <- length(unique(z_est))
  674. z_clust_prop <- calc_clust_prop(z_est, ztrue)
  675. out <- list(ari = ari, R_est = R_est, z_clust_prop = z_clust_prop)
  676. return(out)
  677. }

simu_utils.R at commit 81dfe73, no license · at the source

Overview

Authors: Fangyuan Zhao1,2,3,4, Sheng Wang1,2, Zhikang Wang3,4, Yan Cui3,4, Yi Zhao1,2, Zhiyuan Yuan3,4
ORCID iDs: Fangyuan Zhao, Sheng Wang, Zhikang Wang, Yi Zhao, Zhiyuan Yuan
  1. Institute of Computing Technology, Chinese Academy of Sciences, Beijing, China
  2. University of Chinese Academy of Sciences, Beijing, China
  3. Institute of Science and Technology for Brain-Inspired Intelligence; MOE Key Laboratory of Computational Neuroscience and Brain-Inspired Intelligence; MOE Frontiers Center for Brain Science, Fudan University, Shanghai, China
  4. Center for Medical Research and Innovation, Shanghai Pudong Hospital, Fudan University Pudong Medical Center, Fudan University, Shanghai, China
Institutions: Chinese Academy of Sciences (China); Fudan University (China); Institute of Computing Technology (China); Pudong Medical Center (China); University of Chinese Academy of Sciences (China); Frontiers Center for Brain Science of the Ministry of Education (China); Shanghai Medical College of Fudan University (China)
Journal: Communications biology, volume 9, issue 1, article 725
Dates: received 15 September 2025; accepted 13 March 2026; published online 1 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s42003-026-09923-1 · PMID 41922721 · PMCID PMC13216563 · OpenAlex W7147007029
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), Alzheimer's / dementia (population)
Methods: Connectivity, Machine learning
Keywords: Data mining, Data integration
MeSH: Computational Biology*, Data Mining*, Algorithms, Animals, Humans (* major topic)
Topic: Single-cell and spatial transcriptomics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Shanghai Science and Technology Development Foundation (23YF1403000)
Citations: cited by 1 paper (Europe PMC); 101 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.

Repositories

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

zhengli09/BASS-Analysis

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 81dfe7353093823eed68f403d9c4033031b73de4, 2 April 2023
Languages: R (22), JavaScript (19), Python (1)
Size: 321 files, 42 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, documentation, 7 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration
Tools: tidyverse (6 files), cowplot (3 files), Seurat (3 files), anndata (1 file), ggplot2 (1 file), NumPy (1 file), pandas (1 file), PyTorch (1 file), reticulate (1 file), Scanpy (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
43 files
At the source: github.com/zhengli09/BASS-Analysis

zhaofangyuan98/FDPMining

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: e2cccc564d3acc82c0aeffeb45b8886122dafddd, 7 March 2026
Languages: Python (7), Jupyter (2)
Size: 18 files, 9 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (requirements.txt), 2 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (7 files), Scanpy (6 files), Matplotlib (4 files), pandas (4 files), scikit-learn (3 files), anndata (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
10 files
At the source: github.com/zhaofangyuan98/FDPMining

Zenodo 18902032

License: CC-BY-4.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (7 files), Scanpy (6 files), Matplotlib (4 files), pandas (4 files), scikit-learn (3 files), anndata (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
10 files
At the source: zenodo.org/records/18902032

Code availability statement

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

Read it in the paper: doi.org/10.1038/s42003-026-09923-1.

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:

  • 3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 60 scripts, each with its path and the digest of its content;
  • 1 match 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

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:

Read it in the paper: doi.org/10.1038/s42003-026-09923-1.

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 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 2 keywords, 5 MeSH terms, 1 funder, 91 references.

Cite

This paper

Zhao, F., Wang, S., Wang, Z., Cui, Y., Zhao, Y., & Yuan, Z. (2026). Scalable discovery of spatial multicellular patterns via neighborhood-to-sequence transformation. Communications biology, 9(1), 725. https://doi.org/10.1038/s42003-026-09923-1

BibTeX

@article{zhao2026scalable,
author = {Zhao, Fangyuan and Wang, Sheng and Wang, Zhikang and Cui, Yan and Zhao, Yi and Yuan, Zhiyuan},
title = {{Scalable discovery of spatial multicellular patterns via neighborhood-to-sequence transformation}},
journal = {Communications biology},
year = {2026},
month = apr,
volume = {9},
number = {1},
pages = {725},
publisher = {Nature Publishing Group},
issn = {2399-3642},
doi = {10.1038/s42003-026-09923-1},
url = {https://doi.org/10.1038/s42003-026-09923-1},
pmid = {41922721},
pmcid = {PMC13216563}
}

RIS

TY - JOUR
AU - Zhao, Fangyuan
AU - Wang, Sheng
AU - Wang, Zhikang
AU - Cui, Yan
AU - Zhao, Yi
AU - Yuan, Zhiyuan
TI - Scalable discovery of spatial multicellular patterns via neighborhood-to-sequence transformation
T2 - Communications biology
J2 - Commun Biol
PY - 2026
DA - 2026/04/01
VL - 9
IS - 1
SP - 725
SN - 2399-3642
PB - Nature Publishing Group
DO - 10.1038/s42003-026-09923-1
UR - https://doi.org/10.1038/s42003-026-09923-1
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s42003-026-09923-1",
"type": "article-journal",
"title": "Scalable discovery of spatial multicellular patterns via neighborhood-to-sequence transformation",
"container-title": "Communications biology",
"author": [
{
"family": "Zhao",
"given": "Fangyuan"
},
{
"family": "Wang",
"given": "Sheng"
},
{
"family": "Wang",
"given": "Zhikang"
},
{
"family": "Cui",
"given": "Yan"
},
{
"family": "Zhao",
"given": "Yi"
},
{
"family": "Yuan",
"given": "Zhiyuan"
}
],
"container-title-short": "Commun Biol",
"volume": "9",
"issue": "1",
"page": "725",
"DOI": "10.1038/s42003-026-09923-1",
"PMID": "41922721",
"PMCID": "PMC13216563",
"ISSN": "2399-3642",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s42003-026-09923-1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
1
]
]
}
}

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.1038/s41592-026-03194-8 [code]
Beyond benchmarking: an expert-guided consensus approach to spatially aware clustering.
Journal: Nature methods
In common: anndata, Scanpy, Seurat, 8 other tools, 8 references
[2] 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: reticulate, anndata, Scanpy, 9 other tools, 5 references
[3] 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, Seurat, 8 other tools, 4 references
[4] doi:10.1038/s41467-026-71759-4 [code]
CellNiche represents cellular microenvironments in atlas-scale spatial omics data with contrastive learning.
Journal: Nature communications
In common: anndata, Scanpy, PyTorch, 4 other tools, 7 references
[5] doi:10.1038/s41514-026-00391-9 [code]
Region-specific transcriptional signatures of brain aging in the absence of neuropathology at the single-cell level.
Journal: npj aging
In common: anndata, Scanpy, Seurat, 8 other tools, 4 references
[6] doi:10.1038/s41467-026-71803-3 [code]
Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: reticulate, anndata, Scanpy, 9 other tools, 2 references
[7] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: reticulate, anndata, Scanpy, 7 other tools, 2 references
[8] doi:10.1038/s41467-026-68596-w [code]
Spatial cartography of human thymus enables the geopositioning of lineage transcription factors in rare mimetic thymic epithelial cells.
Journal: Nature communications
In common: anndata, Scanpy, Seurat, 8 other tools, 2 references
[9] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: reticulate, anndata, Scanpy, 8 other tools, 1 reference
[10] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: anndata, Scanpy, Seurat, 7 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.