Run (It is recommended to run on the server)
It is best to save the code below in an R file named scc12.R and then run it on the server.
# ---------------------------- packages
library(stadia)
# ---------------------------- data
load("./data/scc12.RData")
K <- 12
etas <- 0.15
# ---------------------------- model
d <- 35
set.seed(123)
system.time({
## set hyperparameters
hyper <- HyperParameters(scc12, dim = d, eta = etas)
## run model
scc12_stadia <- stadia(scc12, hyper, dim = d, n_cluster = K,
platform = "visium", em.maxiter = 30)
})
# ---------------------------- save result
if(!dir.exists('./result')) dir.create('./result')
save(scc12_stadia, file = "./result/scc12_stadia.RData")Visulation (It is recommended to run on the Rstudio)
First, load the dependent packages and data.
library(dplyr)
library(mclust)
library(ggplot2)
library(patchwork)
library(corrplot)
library(Seurat)
library(stringr)
library(formattable)
library(tibble)
library(purrr)
load("./data/scc12.RData")
load("./result/scc12_stadia.RData")
# self defined color
clusterCol <- c("#998099", "#79E099", "#E64B35FF", "#E7AF56", "#92EA55", "#d7ddbf","#de9f83", "#7494bc", "#7E7AD8", "#87E0D7", "#D6CED0", "#AE44E3")Merge the list of Seurat object scc12 and save the results of STADIA scc12_stadia in it.
# summary result
batch_info <- rep(paste0("P",c(2,5,9,10)), each = 3)
scc12.pp <- mapply(function(x, y) {
x$orig.ident<- y; return(x)}, scc12, batch_info)
scc12.merge <- merge(
x = scc12.pp[[1]],y = scc12.pp[-1],
add.cell.id = names(scc12.pp), project = "scc")
scc12.merge <-
scc12.merge[,scc12.merge@meta.data[,paste0('nFeature_', DefaultAssay(scc12.merge))] > 200]
scc12.merge$stadiaAnnotation <-
as.vector(scc12_stadia$c_vec)
df <- scc12.merge@meta.data %>%
mutate_at("stadiaAnnotation", as.factor)Visualize the spatial domains distribution using barplot.
dfbar <- df %>%
group_by(sample_id, .data[['stadiaAnnotation']]) %>%
summarise(cnt = n()) %>%
mutate(percentage = formattable::percent(cnt / sum(cnt))) %>%
arrange(desc(percentage), .by_group =TRUE)
dfbar$sample_id <- factor(dfbar$sample_id, levels = paste0("P", c(2, 5, 9, 10)))
ggplot( dfbar, aes( x = sample_id, weight = percentage,
fill = factor(stadiaAnnotation)))+
geom_bar( position = "stack") +
scale_fill_manual(values = clusterCol) +
theme(axis.title = element_blank(), legend.title = element_blank()) + NoLegend() 
Visualize the spatial domains identified by STADIA.
se.list <- SplitObject(scc12.merge, split.by = "slice_id")
df.list <- lapply(se.list, function(x){
df <- as.data.frame(cbind(row = x$row, col = x$col,
label = x[['stadiaAnnotation']][,1]))
df$label <- factor(df$label, levels = 1:12)
row.names(df) <- str_extract(rownames(df), "[0-9]{1,2}x[0-9]{1,2}_[0-9]{1,2}")
df <- rownames_to_column(df)
return(df)
})
df.list.1 <- df.list %>%
purrr::map2(names(.), ~mutate(.x, FXN = .y))
df_spatial <- Reduce(rbind, df.list.1)
colnames(df_spatial)[5] <- "slice_id"
sp_list <- lapply(seq_along(df.list), function(i){
obj <- df.list[[i]]
names(clusterCol) <- 1:12
plt <- ggplot(obj, aes(x = row, y = col))+
# geom_point(size = .1,shape = 19) +
geom_point(size = 1.5, shape = 21, stroke = .08,
color = "#00000000", aes(fill = label)) +
scale_fill_manual(values = clusterCol) +
theme_void() +
theme(legend.position = "none")
})
wrap_plots(sp_list[1:12], nrow = 2)