
---
title: "Walkthrough with the simulated dataset 8"
author: "Suoqin Jin, Lihua Zhang"
output: html_document
mainfont: Arial
vignette: >
  %\VignetteIndexEntry{Integrative analysis of single cell multi-omics data using scAI}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

  ```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  root.dir = './'
)
```


This walkthrough outlines the key steps of scAI using the simulated dataset 8. This simulated data consist of paired single-cell RNA-seq and ATAC-seq data, including five imbalanced cell clusters with five clusters in scRNA-seq data and three clusters in scATAC-seq data.

Load the required libraries
```{r message=FALSE,warning=FALSE}
library(scAI)
library(dplyr)
library(cowplot)
library(ggplot2)
```

## Load data
The algorithm takes a list of two digital data matrices as input. Genes/loci should be in rows and cells in columns. rownames and colnames should be included. Before running the scAI model, we need to normalize the data to account for library size and select highly variable features.

```{r}
load("/Users/suoqinjin/Documents/scAI/data/data_simulation8.rda")
X <- data_simulation8$data # List of data matrix
labels.RNA <- data_simulation8$labels.RNA # the true labels of cells, five clusters in scRNA-seq data
labels.ATAC <- data_simulation8$labels.ATAC # the true labels of cells, three clusters in scATAC-seq data
```

## Create a scAI object
```{r}
scAI_outs <- create_scAIobject(raw.data = X)
```
## Preprocess data
Perform quality control to remove low-quality cells and genes, and normalize the data.
Since this is a simulated data, we do not need to normalize the data. Thus we set `assay = NULL`.

```{r, results='asis'}
scAI_outs <- preprocessing(scAI_outs, assay = NULL, minFeatures = 200, minCells = 1,
                           libararyflag = F, logNormalize = F)
```
Add cell information into *pData* slot of the object
```{r}
scAI_outs <- addpData(scAI_outs, pdata = cbind(labels.RNA,labels.ATAC), pdata.name = c("labels.RNA","labels.ATAC"))
```

## Run scAI model
As depending on the random initilization the results might differ, we run scAI multiple times (e.g. nrun = 5) and output the best result. User can also output results from all runs by setting *keep_all = TRUE*. The key parameters here are the number of factors/clusters (k). The `selectK` function can aid in selecting k. A suitable k is the one at which the magnitude of cophenetic correlation begins to fall.
```{r}
scAI_outs <- run_scAI(scAI_outs, K = 5, nrun = 2)
```

## Visualize the inferred biologically relevant components
We plot the heatmap of the three learned low-rank matrices using hierarchical clustering. The ground truth labels of cells are used for validation (not necessary).
```{r, fig.width=7,fig.height = 8, fig.wide = TRUE, fig.align = "center"}
lmHeatmap(scAI_outs, color.by = "labels.RNA")
```

## Visualize cells onto the low-dimensional space
We can visualize cells onto the low-dimensional space using t-SNE, FIt-sne or UMAP.
Here, we perform comparison of the visualization of raw ATAC-seq data with the aggregated data. Cells are colored by the true labels.
```{r, fig.width=7,fig.height = 3.5,  fig.wide = TRUE, fig.align = "center"}
cell_coords.ori <- reducedDims(scAI_outs, data.use = scAI_outs@norm.data$ATAC, do.scale = F, method = "umap", return.object = F)
cell_coords.agg <- reducedDims(scAI_outs, data.use = scAI_outs@agg.data, do.scale = F, method = "umap", return.object = F)

gg1 <- cellVisualization(scAI_outs, cell_coords.ori, color.by = "labels.ATAC",show.legend = F, title = "scATAC-seq")
gg2 <- cellVisualization(scAI_outs, cell_coords.agg, color.by = "labels.ATAC", ylabel = NULL, title = "Aggregated scATAC-seq")
cowplot::plot_grid(gg1, gg2)
```

## Identify enriched features in each factor
```{r}
markers_RNA <- identifyFactorMarkers(scAI_outs, assay = 'RNA', n.top = 5)
markers_ATAC <- identifyFactorMarkers(scAI_outs, assay = 'ATAC',  n.top = 5)
```

### Ranking the features (genes/loci) and show the top markers in each factor
```{r, fig.width=8, fig.height=3,  fig.wide = TRUE, fig.align = "center"}
featureRankingPlot(scAI_outs, assay = 'RNA', feature.show = markers_RNA$markers.top$features, top.p = 0.1, ylabel = "Gene score", ncol = 5)
featureRankingPlot(scAI_outs, assay = 'ATAC', feature.show = markers_ATAC$markers.top$features, top.p = 0.1, ylabel = "Locus score", ncol = 5)

```

## Embedding cells, genes, loci and factors into 2D-dimensions using our new visualization method VscAI

```{r message=FALSE,warning=FALSE}
scAI_outs <- getEmbeddings(scAI_outs)

```

### Visualization of the embedding using VscAI
User can provide a vector of the features (e.g., key marker genes/loci) to explore the biological meaning of the cell groups and enhance the interpretation of the data. Here, we select the top two features of each factor.
```{r, fig.width=10,fig.height=4, fig.align = "center"}
genes.embed <- markers_RNA$markers.top %>% group_by(factors) %>% slice(1:2)
genes.embed <- as.character(genes.embed$features)
loci.embed <- markers_ATAC$markers.top %>% group_by(factors) %>% slice(1:2)
loci.embed <- as.character(loci.embed$features)

gg1 <- VscAIplot(scAI_outs, gene.use = genes.embed, loci.use = NULL, loci.use.names = NULL, color.by = "labels.RNA")
gg2 <- VscAIplot(scAI_outs, gene.use = NULL, loci.use = loci.embed, loci.use.names = loci.embed, color.by = "labels.ATAC")
cowplot::plot_grid(gg1, gg2)

```

## Feature plot
We can overlay the expression of features, or the cell loading values onto the low-dimensional space, e.g., VscAI, tsne, umap
```{r, fig.width=9, fig.height=5,  fig.wide = TRUE, fig.align = "center"}
featureScoreVisualization(scAI_outs, feature.scores = t(scAI_outs@fit$H), feature.use = c('factor1','factor2','factor3','factor4','factor5'),  method = "VscAI", nCol = 3, cell.size = 0.1, show.legend = T, show.legend.combined = F)
```

## Identify cell clusters
We can also identify cell clusters based on the inferred cell loading matrix using Leiden algorithm.
```{r}
scAI_outs <- identifyClusters(scAI_outs, resolution = 0.05)

```
## Visualize cells onto the low-dimensional space
We can visualize cells onto the low-dimensional space generated by t-SNE, FIt-sne or UMAP.
Here, we perform UMAP dimension reduction. Cells are colored by the clustering inferred by scAI.
```{r, fig.width=4.5,fig.height = 3.5,  fig.wide = TRUE, fig.align = "center"}
scAI_outs <- reducedDims(scAI_outs, method = "umap")
cellVisualization(scAI_outs, scAI_outs@embed$umap, color.by = "cluster")

```

