---

title: "Walkthrough with A549 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 ATAC-seq data of A549 cells.

This data describes lung adenocarcinoma-derived A549 cells after 0, 1 and 3 hours of 100 nM dexamethasone treatment, including scRNA-seq and scATAC-seq data of 2641 co-assayed cells.

The data is downloaded from GEO (GSM3271040 and GSM3271041) that were co-assayed using the sci-CAR protocol (Cao, J. et al. Joint profiling of chromatin accessibility and gene expression in thousands of single cells. Science 361, 1380-1385 (2018)).

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_A549.rda")

X <- data_A549$data # List of data matrix

labels <- data_A549$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 = "Time")

```

## 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 = 2, nrun = 5)

```

## Visualize the inferred biologically relevant components

We plot the heatmap of the three learned low-rank matrices using hierarchical clustering. The time information of cells are used for validation.

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

lmHeatmap(scAI_outs, color.by = "Time")

```

## 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 comparison of the visualization of raw ATAC-seq data with the aggregated data. Cells are colored by the collected time.

```{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 = "Time",show.legend = F, title = "scATAC-seq")

gg2 <- cellVisualization(scAI_outs, cell_coords.agg, color.by = "Time", 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 highlighting the important markers in each factor

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

# show the markers of GR activation

feature_genes = c('ZSWIM6','BASP1','TXNRD1','NR3C1','CKB','ABHD12','CDH16','NFKBIA','PER1','SCNN1A');

feature_loci <- c('5-17070007-17070367','5-17113334-17113956','5-17118743-17119143','5-17189512-17190061','5-17216210-17217826','5-17249185-17249758',

                  '20-25370461-25372119','16-66928981-66929756','12-6376107-6377003')

feature_loci_motif <- c('FOXA1','RELA/SP1','SP1','SMAD3','NR3C1/YY1','CEBPB/GATA3',

                        'CEBPB','SMAD3','CREB1/NR3C1')

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

featureRankingPlot(scAI_outs, assay = 'ATAC', feature.show = feature_loci, feature.show.names = feature_loci_motif, top.p = 0.5, ylabel = "Locus score")

```

## 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=10,fig.height=4, fig.align = "center"}

gg1 <- VscAIplot(scAI_outs, gene.use = feature_genes, loci.use = NULL, loci.use.names = NULL, color.by = "Time")

gg2 <- VscAIplot(scAI_outs, gene.use = NULL, loci.use = feature_loci, loci.use.names = feature_loci_motif, color.by = "Time")

cowplot::plot_grid(gg1, gg2)

```

## 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('NR3C1','CKB','ABHD12','NFKBIA'),  method = "VscAI", nCol = 4, cell.size = 0.1, show.legend = F, show.legend.combined = T)

```