---
title: "Feature Selection and Hierarchical Clustering of cells in Zhengmix4eq"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Example analysis}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      warning = FALSE,
                      comment = NA,
                      fig.width = 6,
                      fig.height = 4)
library(HIPPO)
library(SingleCellExperiment)
```

## Installing the Package

HIPPO is under review for Bioconductor and CRAN submission. You can download the developer version as below. Please allow up to 5 minutes to completely compile the vignette.

```{r, eval=FALSE}
devtools::install_github("tk382/HIPPO", build_vignettes = TRUE)
```

## Read the data

The data set is available in the following [link](http://imlspenticton.uzh.ch/robinson_lab/DuoClustering2018/DuoClustering2018.tar.gz), where the detailed explanation is available [here](https://github.com/markrobinsonuzh/scRNAseq_clustering_comparison). Note that the file is very large (3.3GB). We use a subset of Zheng 2017 data set saved as a toydata within this package. Also, we load a reference data that matches ENSG ID to HGNC gene names.

```{r}
data(toydata)
data(ensg_hgnc)
```

Alternatively, you can start from a matrix object.

```{r, eval = FALSE}
# X = readRDS("zhengmix4eq_counts.rds")
# toydata = SingleCellExperiment(assays = list(counts = X))
```

## Diagnostic Plot

This plot shows the zero inflation compared to the expected Poisson line. If most genes don't align with the black line, it shows that there is cell heterogeneity driving the zero inflation.

```{r, warning = FALSE}
hippo_diagnostic_plot(toydata, 
                      show_outliers = TRUE, 
                      zvalue_thresh = 2)
```

## Feature Selection and Hierarchical Clustering

HIPPO assumes that the count matrix is placed in toydata@assays@data\$counts. Some objects that we found online have the count matrix in toydata@assays\$data\$counts. In this case, HIPPO will throw an error because it cannot found a count matrix. In this case, you have to create another SingleCellExperiment object to assign the count matrix in the correct slot.

Next, you can run hippo function to do the pre-processing that simutlaneously conducts feature selection and hierarchcial clustering. There are three arguments that help you decide the stopping criterion of clustering procedure.

K is the maximum number of clusters that you want. HIPPO will return the clustering results for all k = 2, 3, ..., K, so you can overestimate the number of potential clusters. The default is 10, but users are highly recommended to adjust this.

z_threshold is the feature selection criterion. For each round of hierarchical clustering, hippo will find outlier genes where the z-value of significance is greater than the threshold. For example, if you would like to select genes with p-values less than 0.05, z_threshold would be 1.96. The default threshold is 2, but users can use their discretion to change this value.

outlier_proportion is the number of outlier genes to allow. The default is 0.01 (1%) which means the clustering procedure will automatically stop if there are less than 1% of genes remain as important features. With the example data set, the default choice has empirically worked well.

```{r, warning = FALSE}
set.seed(20200321)
toydata = hippo(toydata, K = 10, 
                z_threshold = 2, outlier_proportion = 0.00001)
```

## Dimension Reduction for Each Round of HIPPO

We offer two dimension reduction methods: umap and tsne. And we offer two separate visualization functions.

```{r, warning = FALSE}
toydata = hippo_dimension_reduction(toydata, method="umap")
hippo_umap_plot(toydata)
```

```{r}
toydata = hippo_dimension_reduction(toydata, method="tsne")
hippo_tsne_plot(toydata)
```

## Visualize the selected features at each round

This function shows how the zero-inflation decreases as HIPPO proceeds in the clustering. This function has arguments called switch_to_hgnc and ref. These aim to provide the users an option to change the gene names from ENSG IDs to HGNC symbols for ease of understanding. Many SingleCellExperiment objects have such data embedded in rowData(toydata). Users can create a data frame with ensg and hgnc columns for the genes, and HIPPO will automatically switch the row names of the count matrix from ENSG IDs to HGNC symbols. The default is set to FALSE, assuming that the row names are already HGNC symbols.

```{r, fig.width = 6, fig.height = 3}
data(ensg_hgnc)
zero_proportion_plot(toydata, 
                     switch_to_hgnc = TRUE, 
                     ref = ensg_hgnc)
hippo_feature_heatmap(toydata, k = 3, 
                      switch_to_hgnc = TRUE, 
                      ref = ensg_hgnc, 
                      top.n = 20)
```

## Differential Expression Example

We also offer a differential expression analysis tool.

This function also has an option to switch the gene names to HGNC symbols. top.n argument lets users choose how many top genes to show in the box plot. The default is 5.

The labels of boxplots are aligned with the t-SNE or UMAP plots above. When K is equal to 2, the color codes match with the cell groups as separated in the dimension reduction plot. 

```{r}
toydata = hippo_diffexp(toydata, 
                  top.n = 5, 
                  switch_to_hgnc = TRUE, 
                  ref = ensg_hgnc)
```

Each round of differential expression test results are also saved in the list of data frames.

```{r}
head(get_hippo_diffexp(toydata, 1))
```


```{r}
sessionInfo()
```
