Skip to content

Commit e5754d2

Browse files
committed
add poi script
1 parent 7008d17 commit e5754d2

File tree

1 file changed

+74
-0
lines changed

1 file changed

+74
-0
lines changed
Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
suppressMessages(library(SingleCellExperiment, quietly = TRUE))
2+
suppressMessages(library(scDesign3, quietly = TRUE))
3+
4+
## VIASH START
5+
par <- list(
6+
input = "resources_test/spatialsimbench_mobnew/MOBNEW.rds",
7+
base = "domain",
8+
output = "simulated_dataset.h5ad",
9+
family = "poi",
10+
usebam = FALSE
11+
)
12+
meta <- list(
13+
name = "scdesign3"
14+
)
15+
## VIASH END
16+
reticulate::py_config()
17+
cat("Read input files\n")
18+
input <- anndata::read_h5ad(par$input)
19+
20+
sce <- SingleCellExperiment(
21+
list(counts = Matrix::t(input$layers[["counts"]])),
22+
colData = input$obs
23+
)
24+
25+
cat("scDesign3 Simulation Start\n")
26+
27+
if (par$base == "tissue") {
28+
sce@colData$cell_type <- "cell_type"
29+
} else if (par$base == "domain") {
30+
sce@colData$cell_type <- as.character(sce@colData$spatial_cluster)
31+
} else {
32+
stop("wrong base parameter")
33+
}
34+
35+
sce_simu <- scdesign3(
36+
sce = sce,
37+
assay_use = "counts",
38+
celltype = "cell_type",
39+
pseudotime = NULL,
40+
spatial = c("row", "col"), # spatial location
41+
other_covariates = NULL,
42+
mu_formula = "s(row, col, bs = 'gp', k = 100)",
43+
sigma_formula = "1",
44+
family_use = par$family, # could be change another distribution
45+
n_cores = 1,
46+
usebam = par$usebam,
47+
corr_formula = "1",
48+
copula = "gaussian",
49+
DT = TRUE,
50+
pseudo_obs = FALSE,
51+
return_model = FALSE,
52+
nonzerovar = FALSE,
53+
parallelization = "mcmapply"
54+
)
55+
56+
cat("Generating output file\n")
57+
new_obs <- sce_simu$new_covariate[c("row", "col")]
58+
59+
output <- anndata::AnnData(
60+
layers = list(
61+
counts = Matrix::t(sce_simu$new_count)
62+
),
63+
obs = new_obs,
64+
var = input$var,
65+
uns = c(
66+
input$uns,
67+
list(
68+
method_id = meta$name
69+
)
70+
)
71+
)
72+
73+
cat("Write output files\n")
74+
output$write_h5ad(par$output, compression = "gzip")

0 commit comments

Comments
 (0)