PRECAST: simulation

Prepare the PRECASTObject with preprocessing step.

Next, we use CreatePRECASTObject() to create a PRECASTObject based on the Seurat list object seuList. This function will do three things:

• 1. Filter low-quality spots and genes, controlled by the arguments premin.features and premin.spots, respectively; the spots are retained in raw data (seuList) with at least premin.features number of nonzero-count features (genes), and the genes are retained in raw data (seuList) with at least premin.spots number of spots. To ease presentation, we denote the filtered Seurat list object as data_filter1.
• 2. Select the top 2,000 variable genes (by setting gene.number=2000) for each data batch using FindSVGs() function in DR.SC package for spatially variable genes or FindVariableFeatures() function in Seurat package for highly variable genes. Next, we prioritized genes based on the number of times they were selected as variable genes in all samples and chose the top 2,000 genes. Then denote the Seurat list object as data_filter2, where only 2,000 genes are retained.
• 3. Conduct strict quality control for data_filter2 by filtering spots and genes, controlled by the arguments postmin.features and postmin.spots, respectively; the spots are retained with at least post.features nonzero counts across genes; the features (genes) are retained with at least postmin.spots number of nonzero-count spots. Usually, no genes are filltered because these genes are variable genes.

If the argument customGenelist is not NULL, then this function only does (3) based on customGenelist gene list.

In this simulated dataset, we don’t require to select genes, thus, we set customGenelist=row.names(seuList[[1]]), representing the user-defined gene list. User can retain the raw seurat list object by setting rawData.preserve = TRUE.

## Create PRECASTObject
set.seed(2022)
PRECASTObj <-  CreatePRECASTObject(seuList, customGenelist=row.names(seuList[[1]]))

## User can retain the raw seuList by the following commond.
##  PRECASTObj <-  CreatePRECASTObject(seuList, customGenelist=row.names(seuList[[1]]), rawData.preserve = TRUE)

Add the model setting

Add adjacency matrix list and parameter setting of PRECAST. More model setting parameters can be found in model_set().

## check the number of genes/features after filtering step
PRECASTObj@seulist

## seuList is null since the default value `rawData.preserve` is FALSE.
PRECASTObj@seuList

## Add adjacency matrix list for a PRECASTObj object to prepare for PRECAST model fitting.
PRECASTObj <-  AddAdjList(PRECASTObj, platform = "Visium")

## Add a model setting in advance for a PRECASTObj object. verbose =TRUE helps outputing the information in the algorithm.
PRECASTObj <- AddParSetting(PRECASTObj, Sigma_equal=FALSE, maxIter=30, verbose=TRUE)

Fit PRECAST

For function PRECAST, users can specify the number of clusters \(K\) or set K to be an integer vector by using modified BIC(MBIC) to determine \(K\). For convenience, we give a single K here.

### Given K
set.seed(2022)
PRECASTObj <- PRECAST(PRECASTObj, K=7)

Select a best model and use ARI to check the performance of clustering

## backup the fitting results in resList
resList <- PRECASTObj@resList
# PRECASTObj@resList <- resList
PRECASTObj <- selectModel(PRECASTObj)
true_cluster <- lapply(data_simu, function(x) x$true_cluster)
str(true_cluster)
mclust::adjustedRandIndex(unlist(PRECASTObj@resList$cluster), unlist(true_cluster))

Integrate the two samples by the function IntegrateSpaData.

seuInt <- IntegrateSpaData(PRECASTObj, species='unknown')
seuInt 
## The low-dimensional embeddings obtained by PRECAST are saved in PRECAST reduction slot.