---
title: "Interactively explore and visualize Single Cell RNA seq data"
author: "Jayaram Kancherla, Kazi Tasnim Zinat, Héctor Corrada Bravo"
date: "`r Sys.Date()`"
output: BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{Explore Data using scTreeViz}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
  
```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction

[`scTreeViz`](https://github.com/HCBravoLab/scTreeViz) is a package for interactive visualization and exploration of Single Cell RNA sequencing data. `scTreeViz` provides methods for exploring hierarchical features (eg. clusters in single cell at different resolutions or taxonomic hierarchy in single cell datasets), while supporting other useful data visualization charts like heatmaps for expression and scatter plots for dimensionality reductions like UMAP or TSNE.

## Loading required packages

```{r load-packages, message=FALSE, warning=FALSE}
library(scTreeViz)
library(Seurat)
library(SC3)
library(scran)
library(scater)
library(clustree)
library(igraph)
library(scRNAseq)
```

# Preparing Datasets

The first step in using the `scTreeViz` package is to wrap datasets into `TreeViz` objects. The `TreeViz` class extends `SummarizedExperiment` and provides various methods to interactively perform various operations on the underlying hierarchy and count or expression matrices. In this section, we show various ways to generate a `TreeViz` object either from existing Single Cell packages (SingleCellExperiment or Seurat) or from a raw count matrix and cluster hierarchy.

## From `SingleCellExperiment`

A number of Single cell datasets are available as `SingleCellExperiment` objects through the `scRNAseq` package, for this usecase, we use `LunSpikeInData` dataset. In addition, we calculate the dimensionality reductions; UMAP, TSNE and PCA from the functions provided in `scater` package. 

```{r, results='hide', warning=FALSE, error=FALSE, message=FALSE}
# load dataset
sce<- ZeiselBrainData()
# Normalization
sce <- logNormCounts(sce)
# calculate umap and tsne
sce <- runUMAP(sce)
sce<- runTSNE(sce)
sce<- runPCA(sce)
```

We provide `createFromSCE` function to create a `TreeViz` object from `SingleCellExperiment` object. Here, the workflow works in two ways:

1. If no cluster information is available in the `colData` of the `SingleCellExperiment` object, we create clusters at different resolutions using the `WalkTrap` algorithm by calling an internal function `generate_walktrap_hierarchy` and use this cluster information for visualization.

```{r, warning=FALSE, error=FALSE, message=FALSE}
treeViz <- createFromSCE(sce, reduced_dim = c("UMAP","PCA","TSNE"))
plot(treeViz)
```

2. If cluster information is provided in the `colData` of the object, then the user should set the flag parameter `check_coldata` to `TRUE` and provide prefix for the columns where cluster information is stored.

```{r, eval=FALSE, warning=FALSE, error=FALSE, message=FALSE}
# Forming clusters
set.seed(1000)
for (i in  seq(10)) {
  clust.kmeans <- kmeans(reducedDim(sce, "TSNE"), centers = i)
  sce[[paste0("clust", i)]] <- factor(clust.kmeans$cluster)
}

treeViz<- createFromSCE(sce, check_coldata = TRUE, col_regex = "clust")
plot(treeViz)
```

Note: In both cases the user needs to provide the name of dimensionality reductions present in the object as a parameter.

## From `Seurat`

We use the dataset `pbmc_small` available through Seurat to create a `TreeViz` object.  

```{r, echo=TRUE, results='hide', warning=FALSE, error=FALSE, message=FALSE}
data(pbmc_small)
pbmc <- pbmc_small
```

We then preprocess the data and find clusters at different resolutions.

```{r, echo=TRUE, results='hide', warning=FALSE, error=FALSE, message=FALSE}
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- NormalizeData(pbmc)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0), print.output = 0, save.SNN = TRUE)
pbmc
```

The measurements for dimensionality reduction methods we want to visualize are also added to the object via native functions in `Seurat`. Since `PCA` is already added, we calculate `TSNE` and `UMAP`
```{r, echo=TRUE, results='hide', warning=FALSE, error=FALSE, message=FALSE}
# pbmc<- RunTSNE(pbmc)
pbmc<- RunUMAP(pbmc, dims=1:3)
Reductions(pbmc)
```

We use the `createFromSeurat` function to create a `TreeViz` object from `Seurat` object. In addition the object, we pass the name of dimensionality reductions present in the object as a paramter in vector format to indicate these measurements should be added to `treeviz` for visualization. If the mentioned reduced dimension is not present it would simply be ignored.
```{r, echo=TRUE,  warning=FALSE, error=FALSE, message=FALSE}
treeViz<- createFromSeurat(pbmc, check_metadata = TRUE, reduced_dim = c("umap","pca","tsne"))
plot(treeViz)
```

## Create TreeViz from count matrix and Cluster hierarchy

```{r, results='hide', warning=FALSE, error=FALSE, message=FALSE}
n=64
# create a hierarchy
df<- data.frame(cluster0=rep(1,n))
for(i in seq(1,5)){
  df[[paste0("cluster",i)]]<- rep(seq(1:(2**i)),each=ceiling(n/(2**i)),len=n)
}

# generate a count matrix
counts <- matrix(rpois(6400, lambda = 10), ncol=n, nrow=100)
colnames(counts)<- seq(1:64)
# create a `TreeViz` object
treeViz <- createTreeViz(df, counts)
plot(treeViz)
```


## Start the TreeViz App (using hosted app)

Start the App from the `treeViz` object we created. This adds a `facetZoom` to navigate the cluster hierarchy, a heatmap of the top `n` most variable genes from the dataset, where 'n' is selected by the user and one scatter plot for each of the reduced dimensions. 
 
```{r, eval=FALSE, echo=TRUE}
app <- startTreeviz(treeViz, top_genes = 500)
```
![Cell 416B dataset](images/cell416B.png)

![pbmc dataset](images/Treeviz.png)

Users can also use the interface to explore the same dataset using different visualizations available through Epiviz.

## Visualize gene expression across clusters

Users can also add Gene Box plots using either the frontend application, or from R session. In the following example, we visualize the 5th, 50th and 500th most variable gene as Box plots
![visualizing expression of a gene across clusters](images/Boxplot.png)

### Adding Gene Box Plots via UI

Users need to select  `Add Visualization -> Gene Box PLot` option from menu and then select the desired gene using the search pane in the appeared dialogue box

![Selecting a gene to visualize](images/dest_boxplot_dialogue.png)

### Adding Gene Box Plots via R Session

Users can also select the gene from R session by using the `plotGene` command followed by Gene name.
```{r, eval=FALSE}
app$plotGene(gene="AIF1")
```

## Stop App

After exploring the dataset, this command the websocket connection. 

```{r, eval=FALSE, echo=TRUE}
app$stop_app()
```

## Session

```{r}
sessionInfo()