---
title: 'PRECAST: DLPFC Single Sample Analysis'
author: "Wei Liu"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{PRECAST: DLPFC Single Sample Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
  
```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  tidy = TRUE,
  fig.width = 11,
  tidy.opts = list(width.cutoff = 95),
  message = FALSE,
  warning = FALSE,
  time_it = TRUE
)
all_times <- list()  # store the time for each chunk
knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      now <<- Sys.time()
    } else {
      res <- difftime(Sys.time(), now, units = "secs")
      all_times[[options$label]] <<- res
    }
  }
}))
```
This vignette introduces the PRECAST workflow for the analysis of single spatial transcriptomics dataset. The workflow consists of three steps

* Independent preprocessing and model setting
* Probabilistic embedding and clustering  using PRECAST model
* Downstream analysis (i.e. visualization of clusters and embeddings)

We demonstrate the use of PRECAST to one human dorsolateral prefrontal cortex Visium data that are [here](https://github.com/feiyoung/PRECAST/tree/main/vignettes_data), which can be downloaded to the current working path by the following command:
```{r eval=FALSE}
githubURL <- "https://github.com/feiyoung/PRECAST/blob/main/vignettes_data/dlpfc_151672.rda?raw=true"
download.file(githubURL,"dlpfc_151672.rda",mode='wb')

```

Then load to R 
```{r  eval = FALSE}
load("dlpfc_151672.rda")
```


The package can be loaded with the command:
```{r  eval = FALSE}
library(PRECAST)
library(Seurat)
```


## Compare PRECAST with DR-SC in analyzing one sample
First, we view the the spatial transcriptomics data with Visium platform.
```{r  eval = FALSE}
dlpfc_151672 ## a list including two Seurat object
```

check the meta data that must include the spatial coordinates named "row" and "col", respectively.
If the names are not, they are required to rename them.
```{r  eval = FALSE}
meta_data <- dlpfc_151672@meta.data
all(c("row", "col") %in% colnames(meta_data)) ## the names are correct!
head(meta_data[,c("row", "col")])
```

### Prepare the PRECASTObject.

Create a PRECASTObj object to prepare for PRECAST models. Here, we only show the `HVGs` method to select the 2000 highly variable genes, but users are able to choose more genes or also use the `SPARK-X` to choose the spatially variable genes.
```{r  eval = FALSE}
set.seed(2023)
library(PRECAST)
preobj <- CreatePRECASTObject(seuList = list(dlpfc_151672), selectGenesMethod="HVGs",
                              gene.number = 2000) # 
```
### Add the model setting
```{r  eval = FALSE}
## check the number of genes/features after filtering step
preobj@seulist
## Add adjacency matrix list for a PRECASTObj object to prepare for PRECAST model fitting.
PRECASTObj <-  AddAdjList(preobj, 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, coreNum = 1,
                            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$. Here, we use user-specified number of clusters. 
```{r  eval = FALSE}
### Given K
PRECASTObj <- PRECAST(PRECASTObj, K= 7)
```

Use the function `SelectModel()` to re-organize the fitted results in PRECASTObj.
```{r  eval = FALSE}
## backup the fitting results in resList
resList <- PRECASTObj@resList
PRECASTObj <- SelectModel(PRECASTObj)
ari_precast <- mclust::adjustedRandIndex(PRECASTObj@resList$cluster[[1]], PRECASTObj@seulist[[1]]$layer_guess_reordered)
```


<details>
<summary>**Other options**</summary>
Users are also able to set multiple K, then choose the  best one
```{r  eval =FALSE}
PRECASTObj2 <- AddParSetting(PRECASTObj, Sigma_equal = FALSE, coreNum = 4,
                            maxIter=30, verbose = TRUE) # set 4 cores to run in parallel.
PRECASTObj2 <- PRECAST(PRECASTObj2, K= 5:8)
## backup the fitting results in resList
resList2 <- PRECASTObj2@resList
PRECASTObj2 <- SelectModel(PRECASTObj2)
str(PRECASTObj2@resList)
mclust::adjustedRandIndex(PRECASTObj2@resList$cluster[[1]], PRECASTObj2@seulist[[1]]$layer_guess_reordered)
```

Besides, user can also use different initialization method by setting `int.model`, for example, set `int.model=NULL`; see the functions `AddParSetting()` and `model_set()` for more details.

</details>



Put the reults into a Seurat object seuInt.
```{r  eval = FALSE}
seuInt <- PRECASTObj@seulist[[1]]
seuInt@meta.data$cluster <- factor(unlist(PRECASTObj@resList$cluster))
seuInt@meta.data$batch <- 1
seuInt <- Add_embed(PRECASTObj@resList$hZ[[1]], seuInt, embed_name = 'PRECAST')
posList <- lapply(PRECASTObj@seulist, function(x) cbind(x$row, x$col))
seuInt <- Add_embed(posList[[1]], seuInt, embed_name = 'position')
Idents(seuInt) <- factor(seuInt@meta.data$cluster)

seuInt 
## The low-dimensional embeddings obtained by PRECAST are saved in PRECAST reduction slot.
```
Save the spatial and tSNE scatter plots for clusters  from PRECAST
```{r  eval = FALSE}
p_sp1 <- SpaPlot(seuInt, item='cluster', point_size = 3, combine = F)[[1]] + cowplot::theme_cowplot() + 
  ggplot2::ggtitle(paste0("PRECAST: ARI=", round(ari_precast, 2)) ) +
  ggplot2::xlab("row") + ggplot2::ylab("col")
seuInt <- AddTSNE(seuInt,n_comp = 2)
p_tsne <- dimPlot(seuInt, item='cluster')
p_tsne <- p_tsne + cowplot::theme_cowplot() + ggplot2::ggtitle("PRECAST")
```



Fit DR-SC and Plot the spatial and tSNE scatter plots for clusters 
```{r  eval = FALSE}
seu_drsc <- DR.SC::DR.SC(PRECASTObj@seulist[[1]], K=7, verbose=T)
ari_drsc <- mclust::adjustedRandIndex(seu_drsc$spatial.drsc.cluster, PRECASTObj@seulist[[1]]$layer_guess_reordered)
p_tsne_drsc <- DR.SC::drscPlot(seu_drsc)
p_tsne_drsc <- p_tsne_drsc + ggplot2::ggtitle("DR-SC")
p_sp2 <- DR.SC::spatialPlotClusters(seu_drsc)+ cowplot::theme_cowplot()  + 
  ggplot2::ggtitle(paste0("DR-SC ARI=", round(ari_drsc, 2)) ) 
```

Compare the clustering performance of PRECAST and DR-SC.

```{r  eval = FALSE, fig.width=10, fig.height=4}
drawFigs(list(p_sp1, p_sp2), layout.dim = c(1,2))

```


Compare the tSNE visualiztion performance of PRECAST and DR-SC.
```{r  eval = FALSE, fig.width=10, fig.height=4}
drawFigs(list(p_tsne, p_tsne_drsc), layout.dim = c(1,2))
```

<details>
<summary>**Session Info**</summary>
```{r}
sessionInfo()
```
</details>