---

title: "Walkthrough with mESC dataset"

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 will perform integrative analysis of paired single cell RNA-seq and DNA methylation data of mouse embryonic development.

This data describes mouse embryonic stem cells that are cultured in "2i" and ''serum" conditions, including 77 cells profiled by parallel single cell methylation and transcriptome sequencing technique scM&T-seq.

Load the required libraries

```{r message=FALSE,warning=FALSE}

library(scAI)

library(dplyr)

library(cowplot)

```

## 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. Our preprocessed data matrices after normalization and feature selection were provided for this walkthrough.

```{r}

load("/Users/suoqinjin/Documents/scAI/data/data_mESC.rda")

X <- data_mESC$data # List of data matrix

labels <- data_mESC$labels # the collected time of cells, which is used for validation

```

## 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 preprocessed 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 = labels, pdata.name = "Conditions")

```

## 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 component/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 = 3, nrun = 5,alpha = 0.01,lambda = 1000, gamma = 100000)

```

## 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 = 1)

levels(scAI_outs@identity)

```

## 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 either the published cell labels or the clustering inferred by scAI.

```{r, fig.width=7,fig.height = 3,  fig.wide = TRUE, fig.align = "center"}

scAI_outs <- reducedDims(scAI_outs, method = "umap",do.scale = F)

gg1 <- cellVisualization(scAI_outs, scAI_outs@embed$umap, color.by = "Conditions",show.legend = T, title = "Conditions")

gg2 <- cellVisualization(scAI_outs, scAI_outs@embed$umap, color.by = "cluster", ylabel = NULL, title = "scAI clusters")

cowplot::plot_grid(gg1, gg2)

```

Here, we perform comparison of the visualization of raw DNA-seq data with the aggregated data. Cells are colored by the collected time.

```{r, fig.width=9,fig.height = 3,  fig.wide = TRUE, fig.align = "center"}

cell_coords.RNA <- reducedDims(scAI_outs, data.use = scAI_outs@norm.data$RNA, do.scale = T, method = "pca", return.object = F)

cell_coords.DNA <- reducedDims(scAI_outs, data.use = scAI_outs@norm.data$DNA, do.scale = T, method = "pca", return.object = F)

cell_coords.DNAagg <- reducedDims(scAI_outs, data.use = scAI_outs@agg.data, do.scale = T, method = "pca", return.object = F)

gg1 <- cellVisualization(scAI_outs, cell_coords.RNA, color.by = "cluster",  show.legend = F, title = "scRNA-seq",xlabel = "PCA1", ylabel = "PCA2")

gg2 <- cellVisualization(scAI_outs, cell_coords.DNA, color.by = "cluster",show.legend = F, xlabel = "PCA1",ylabel = NULL,title = "scDNA-seq")

gg3 <- cellVisualization(scAI_outs, cell_coords.DNAagg, color.by = "cluster", xlabel = "PCA1",ylabel = NULL, title = "Aggregated scDNA-seq")

cowplot::plot_grid(gg1, gg2, gg3, ncol = 3)

```

## 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=2.5,  fig.wide = TRUE, fig.align = "center"}

featureScoreVisualization(scAI_outs, feature.scores = t(scAI_outs@fit$H), feature.use = c('factor1','factor2','factor3'),  method = "umap", nCol = 3, cell.size = 0.1, show.legend = T, show.legend.combined = F)

```

### Ranking the features (genes/loci) and highlighting the important markers in each factor

```{r, fig.width=6, fig.height=3,  fig.wide = TRUE, fig.align = "center"}

# show the markers of GR activation

feature_genes = c('Zfp42','Esrrb','Morc1','Fbxo15','Jam2','Klf4','Tcl1','Tbx3',

                  'Tex19.1','Krt8','Cald1','Anxa5','Tagln','Ahnak','Dsp','Anxa3','Krt19','Fgf5');

featureRankingPlot(scAI_outs, assay = 'RNA', feature.show = feature_genes, top.p = 0.5, ylabel = "Gene score")

```

### Identify differentially expressed features of each cell cluster

```{r}

markers.RNA.cluster <- identifyClusterMarkers(scAI_outs, assay = "RNA")

markers.DNA.cluster <- identifyClusterMarkers(scAI_outs, assay = 'DNA')

```

### Generate a heatmap to show the expression patterns of top features in cell clusters

```{r, fig.width=10,fig.height = 4, fig.align = "center"}

n.top = 10

markers.RNA.clusterTop <- markers.RNA.cluster %>% group_by(clusters) %>% top_n(n.top, logFC) %>% slice(1:n.top)

featureHeatmap(scAI_outs, assay = "RNA", feature.use = markers.RNA.clusterTop$features, group.by = "cluster")

markers.DNA.clusterTop <- markers.DNA.cluster %>% group_by(clusters) %>% top_n(n.top, logFC) %>% slice(1:n.top)

featureHeatmap(scAI_outs, assay = "DNA", feature.use = markers.DNA.clusterTop$features, group.by = "cluster")

```

## 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 embedded cells, genes, loci and factors

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. 

```{r, fig.width=5.5,fig.height=4, fig.align = "center"}

VscAIplot(scAI_outs, gene.use = feature_genes, loci.use = NULL, loci.use.names = NULL, color.by = "cluster")

```

## Feature plot

We can overlay the expression of features onto the low-dimensional space, e.g., VscAI, tsne, umap

```{r, fig.width=9, fig.height=2.5,  fig.wide = TRUE, fig.align = "center"}

featureVisualization(scAI_outs, assay = "RNA", feature.use = c('Tcl1','Krt19','Dsp','Fgf5'),  method = "umap", nCol = 4, cell.size = 0.1, show.legend = F, show.legend.combined = F)

```