FAST: simulation

View the simulated data

First, we view the the three simulated spatial transcriptomics data with ST platform. There are 200 genes for each data batch and ~2000 spots in total

simu3 ## a list including three Seurat object with default assay: RNA

Check the content in simu3.

head(simu3[[1]])

row.names(simu3[[1]])[1:10]

Create a PRECASTObject object

We show how to create a PRECASTObject object step by step. First, we create a Seurat list object using the count matrix and meta data of each data batch. Although simu3 is a prepared Seurat list object, we re-create a same objcet seuList to show the details.

• Note: the spatial coordinates must be contained in the meta data and named as row and col, which benefits the identification of spaital coordinates by FAST

## Get the gene-by-spot read count matrices
countList <- lapply(simu3, function(x) x[["RNA"]]@counts)

## Check the spatial coordinates: Yes, they are named as "row" and "col"!
head(simu3[[1]]@meta.data)

## Get the meta data of each spot for each data batch
metadataList <- lapply(simu3, function(x) x@meta.data)


## ensure the row.names of metadata in metaList are the same as that of colnames count matrix in countList
M <- length(countList)
for(r in 1:M){
  row.names(metadataList[[r]]) <- colnames(countList[[r]])
}


## Create the Seurat list  object

seuList <- list()
for(r in 1:M){
  seuList[[r]] <- CreateSeuratObject(counts = countList[[r]], meta.data=metadataList[[r]], project = "FASTsimu")
}

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

## User can retain the raw seuList by the following commond.
##  PRECASTObj <-  CreatePRECASTObject(seuList, customGenelist=row.names(seuList[[1]]), rawData.preserve = TRUE)
## check the number of genes/features after filtering step
PRECASTObj@seulist

Fit FAST using simulated data

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

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

## Add a model setting in advance for a PRECASTObj object: verbose =TRUE helps outputing the information in the algorithm; 
PRECASTObj <- AddParSettingFAST(PRECASTObj, verbose=TRUE)
## Check the parameters
PRECASTObj@parameterList

### set q= 20 here
set.seed(2023)
PRECASTObj <- FAST(PRECASTObj, q=20)
### Check the results
str(PRECASTObj@resList)

set.seed(2023)
PRECASTObj <- FAST(PRECASTObj, q=20, fit.model='gaussian')
### Check the results
str(PRECASTObj@resList)

## Obtain the true labels
yList <- lapply(PRECASTObj@seulist, function(x) x$Group)
### Evaluate the MacR2
MacVec <- sapply(1:length(PRECASTObj@seulist), function(r) get_r2_mcfadden(PRECASTObj@resList$FAST$hV[[r]], yList[[r]]))
### output them
print(MacVec)

Embedding alignment and clustering using Harmony and Louvain

Based on the embeddings from FAST, we use Harmony to align the embeddings then followed by Louvain clustering to obtain the cluster labels. In this downstream analysis, other methods for embedding alignment and clustering can be also used. In the vignette of two sections of DLPFC Visium data, we will show another method (iSC-MEB) to jointly perform embedding alignment and spatial clustering.

PRECASTObj <- RunHarmonyLouvain(PRECASTObj, resolution = 0.4)
ARI_vec <- sapply(1:M, function(r) mclust::adjustedRandIndex(PRECASTObj@resList$Louvain$cluster[[r]], yList[[r]]))
print(ARI_vec)

Remove unwanted variations in the log-normalized expression matrices

In the following, we remove the unwanted variations in the log-normalized expression matrices to obtain a combined log-normalized expression matrix in a Seurat object. In the context of the simulated data used in this study, housekeeping genes are not present, thus, we turn to another method to remove the unwanted variations. Specifically, we leverage the batch effect embeddings estimated through Harmony to capture and mitigate unwanted variations. Additionally, we utilize the cluster labels obtained via Louvain clustering to retain the desired biological effects.

The estimated embeddings of batch effects (batchEmbed) are in the slot PRECASTObj@resList$Harmony and cluster labels (cluster) are in the slot PRECASTObj@resList$Louvain.

str(PRECASTObj@resList$Harmony)
str(PRECASTObj@resList$Louvain)

Then, we integrate the three sections by removing the unwanted variations and setting seulist_HK=NULL and Method = "HarmonyLouvain" in the function IntegrateSRTData(). After obtaining seuInt, we will see there are three embeddings: FAST, harmony and position, in the slot seuInt@reductions. FAST are the embeddings obtained by FAST model fitting and are uncorrected embeddings that may includes the unwanted effects (i.e., batch effects); harmonyare the embeddings obtained by Harmony fitting and are aligned embeddings; and position are the spatial coordinates.

seulist_HK <- NULL
seuInt <- IntegrateSRTData(PRECASTObj, seulist_HK=seulist_HK, Method = "HarmonyLouvain", seuList_raw=NULL, covariates_use=NULL, verbose=TRUE)
seuInt

Visualization

First, user can choose a beautiful color schema using chooseColors() in the R package PRECAST.

cols_cluster <- chooseColors(palettes_name = "Nature 10", n_colors = 7, plot_colors = TRUE)

Then, we plot the spatial scatter plot for clusters using the function SpaPlot() in the R package PRECAST.

p12 <- SpaPlot(seuInt, item = "cluster", batch = NULL, point_size = 1, cols = cols_cluster, combine = TRUE)
p12

Users can re-plot the above figures for specific need by returning a ggplot list object. For example, we plot the spatial heatmap using a common legend by using the function drawFigs() in the R package PRECAST.

pList <- SpaPlot(seuInt, item = "cluster", title_name= 'Section',batch = NULL, point_size = 1, cols = cols_cluster, combine = FALSE)
drawFigs(pList, layout.dim = c(1, 3), common.legend = TRUE, legend.position = "right", align = "hv")

We use the function AddUMAP() in the R package PRECAST to obtain the three-dimensional components of UMAP using the aligned embeddings in the reduction harmony.

seuInt <- AddUMAP(seuInt, n_comp=3, reduction = 'harmony', assay = 'RNA')
seuInt

We plot the spatial tNSE RGB plot to illustrate the performance in extracting features.

p13 <- SpaPlot(seuInt, batch = NULL, item = "RGB_UMAP", point_size = 1, combine = FALSE, text_size = 15)
drawFigs(p13, layout.dim = c(1, 3), common.legend = TRUE, legend.position = "right", align = "hv")

We use the function AddUTSNE() in the R package PRECAST to obtain the two-dimensional components of UMAP using the aligned embeddings in the reduction harmony, and plot the tSNE plot based on the extracted features to check the performance of integration.

seuInt <- AddTSNE(seuInt, n_comp = 2, reduction = 'harmony', assay = 'RNA')
p1 <- dimPlot(seuInt, item = "cluster", point_size = 0.5, font_family = "serif", cols = cols_cluster,
    border_col = "gray10",  legend_pos = "right")  # Times New Roman
p2 <- dimPlot(seuInt, item = "batch", point_size = 0.5, font_family = "serif", legend_pos = "right")

drawFigs(list(p1, p2), layout.dim = c(1, 2), legend.position = "right")

Combined differential epxression analysis

Finally, we condut the combined differential expression analysis using the integrated log-normalized expression matrix saved in the seuInt object. The function FindAllMarkers() in the Seurat R package is ued to achieve this analysis.

dat_deg <- FindAllMarkers(seuInt)
library(dplyr)
n <- 5
dat_deg %>%
    group_by(cluster) %>%
    top_n(n = n, wt = avg_log2FC) -> top10
seuInt <- ScaleData(seuInt)

Plot dot plot of normalized expressions for each spatial domain identified by using the FAST embeddings.

col_here <- c("#F2E6AB", "#9C0141") 
library(ggplot2)
p1 <- DotPlot(seuInt, features=unname(top10$gene), cols=col_here, #  idents = ident_here,
              col.min = -1, col.max = 1) + coord_flip()+ theme(axis.text.y = element_text(face = "italic"))+
  ylab("Domain") + xlab(NULL) + theme(axis.text.x = element_text(size=12, angle = 25, hjust = 1, family='serif'),
                                      axis.text.y = element_text(size=12, face= "italic", family='serif'))
p1

sessionInfo()
#> R version 4.2.1 (2022-06-23 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22621)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=C                               
#> [2] LC_CTYPE=Chinese (Simplified)_China.utf8   
#> [3] LC_MONETARY=Chinese (Simplified)_China.utf8
#> [4] LC_NUMERIC=C                               
#> [5] LC_TIME=Chinese (Simplified)_China.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.33   R6_2.5.1        jsonlite_1.8.7  evaluate_0.21  
#>  [5] cachem_1.0.8    rlang_1.1.1     cli_3.4.1       rstudioapi_0.14
#>  [9] jquerylib_0.1.4 bslib_0.5.0     rmarkdown_2.23  tools_4.2.1    
#> [13] xfun_0.39       yaml_2.3.7      fastmap_1.1.1   compiler_4.2.1 
#> [17] htmltools_0.5.5 knitr_1.43      sass_0.4.7