BNPSpace Tutorial

Example - Mouse ofactory Bulb (MOB) ST data

Import the required R packages and other useful functions

suppressMessages(library(SingleCellExperiment))
suppressMessages(library(Matrix))
suppressMessages(library(flexclust))
suppressMessages(library(scater))
suppressMessages(library(scran))
source("R/utils.R")

Load data

The current ST data is stored in a SingleCellExperiment object.

set.seed(1)
load("application/MOB/data/MOB_sce.RData")
sce

## class: SingleCellExperiment 
## dim: 16034 282 
## metadata(0):
## assays(2): counts logcounts
## rownames(16034): Plekha1 Cnih3 ... St14 Gdf5
## rowData names(1): is.HVG
## colnames(282): 1 2 ... 281 282
## colData names(4): row col label sizeFactor
## reducedDimNames(1): PCA
## mainExpName: NULL
## altExpNames(0):

Filtering

We filter spots with library size less than 100, and genes with small reads.

sce = spatialFilter(sce, cutoff_max = 10)
sce

## class: SingleCellExperiment 
## dim: 1117 278 
## metadata(0):
## assays(2): counts logcounts
## rownames(1117): Plekha1 Ctsb ... Cyp2a5 Bpifb9a
## rowData names(1): is.HVG
## colnames(278): 1 2 ... 281 282
## colData names(4): row col label sizeFactor
## reducedDimNames(1): PCA
## mainExpName: NULL
## altExpNames(0):

Construct adjacent matrix

The neighbors is a \(278 \times 4\) matrix, representing the neighbors of each spot. If the number of neighbors less than \(4\), the index is \(0\).

### construct spots network
Adj = find_neighbors(sce, "ST", "lattice")
neighbors = find_neighbor_index(Adj, "ST")
head(neighbors)

##      [,1] [,2] [,3] [,4]
## [1,]    2   19  150  188
## [2,]    1    4  152  187
## [3,]    4    6  148  190
## [4,]    2    3  144  189
## [5,]    6    8  130  191
## [6,]    3    5  132  192

Run BNPSpace model

If the parameter n_clusters is not NULL and given a positive integer value, it is the finite mixture model. The default value is NULL, the MRF-constrained MFM model. K_init is the initial number of clusters and d is the hyperparameter in the MRF model.

source("R/main.R")
timestart = Sys.time()
result1 = NSCFS(sce, neighbors, n_clusters = NULL, d = 1, K_init = 10, n_iters = 300, seed = 1)

## Warning in runMCMC(group_t - 1, D_t, gamma_t, D_0_t, R_t, pi_t, f = d, a, :
## When called from R, the RNG seed has to be set at the R level via set.seed()

timeend = Sys.time()
timerun = timeend - timestart
print(timerun)

## Time difference of 1.774684 mins

pred_label = result1$pred_label

Results

The clustering result

ARI = randIndex(table(result1$pred_label[!is.na(colData(sce)$label)], colData(sce)$label[!is.na(colData(sce)$label)]))
library(BayesSpace)

## 
## Attaching package: 'BayesSpace'

## The following object is masked _by_ '.GlobalEnv':
## 
##     spatialPreprocess

metadata(sce)$BayesSpace.data <- list()
metadata(sce)$BayesSpace.data$platform <- "ST"
metadata(sce)$BayesSpace.data$is.enhanced <- FALSE
p <- clusterPlot(sce, label= pred_label, palette=NULL, size=0.05) +
    #scale_fill_viridis_d(option = "A", labels = 1:clustNumMat[j,2]) +
    labs(title=paste0("BNPSpace: ARI=", round(ARI, 3)), fill = "Domains") + 
    theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
          legend.key.height = unit(0.5, 'cm'), #change legend key height
          legend.key.width = unit(0.5, 'cm'), #change legend key width
          legend.title = element_text(size=12), #change legend title font size
          legend.text = element_text(size=10),#change legend text font size
          panel.border = element_rect(colour = "black", fill=NA, size=1),
          plot.title = element_text(size = 16, face = "bold", hjust = 0.5))

## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

p