---
title: "The FEAST User's Guide"
author:  
  - name: Kenong Su
    email: kenong.su@emory.edu
  - name: Hao Wu
    email: hao.wu@emory.edu
shorttitle: FEAST guide
bibliography: FEAST.bib
package: FEAST
abstract: >
  This vignette introduces the usage of the Bioconductor package FEAST (FEAture SelecTion for scRNA-seq data), which is specifically designed for selecting most representative genes before performing the core of clustering. It is demonstrated an improved clustering accuracy when recruting the featured genes selected by FEAST. 
vignette: >
  %\VignetteIndexEntry{The FEAST User's Guide}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\usepackage[utf8]{inputenc}
output:
  BiocStyle::html_document:
    toc: true
    toc_float:
      collapsed: true
      smooth_scroll: true
    fig_width: 5
---
<style type="text/css">
  pre:not([class]) {
    background-color: white;
  }
  .tocify-subheader > .tocify-item {
  text-indent: 2px;}
</style>


```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


\vspace{.1in}
# Installation and help

## Install FEAST
To install this package, start R (version  > "4.0") and enter:

```{r quickYo, eval = FALSE}
if(!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("FEAST")
```

## Help for FEAST
If you have any FEAST-related  questions, please post to the GitHub Issue section of FEAST at https://github.com/suke18/FEAST/issues, which will be helpful for the construction of FEAST. 

# Introduction
## Background

Cell clustering is one of the most important and commonly performed tasks in single-cell RNA sequencing (scRNA-seq) data analysis. An important step in cell clustering is to select a subset of genes (referred to as “features”), whose expression patterns will then be used for downstream clustering. A good set of features should include the ones that distinguish different cell types, and the quality of such set could have significant impact on the clustering accuracy. All existing scRNA-seq clustering tools include a feature selection step relying on some simple unsupervised feature selection methods, mostly based on the statistical moments of gene-wise expression distributions. In this work, we develop a feature selection algorithm named FEAture SelecTion (FEAST), which provides more representative features. 

## Citation
- For consensus clustering, cite @strehl2002cluster.
- For using (**FEAST**) feature selection procedure, @su2021accurate 

## Quick start
We use `Yan` dataset[@yan2013single] for the demonstration. This `Yan` dataset includes 6 cell types about human preimplantation embryos and embryonic stem cells. The users can run the **FEAST** or **FEAST_fast** function to obtain the gene rankings (from the most significant to the least significant genes) of These featured genes are deemed to better contribute to the clustering accuracy. Here, we demonstrate that some of the top 9 selected genes are very informative (or known as marker genes e.g., CYYR1 ). 
```{r quick, eval = TRUE, message=FALSE, results='hide', include = TRUE}
library(FEAST)
data(Yan)
k = length(unique(trueclass))
Y = process_Y(Y, thre = 2)
# The core function. 
ixs = FEAST(Y, k=k)
# look at the features
Ynorm = Norm_Y(Y)
par(mfrow = c(3,3))
for (i in 1:9){
  tmp_ix = ixs[i]
  tmp_gene = rownames(Ynorm)[tmp_ix]
  boxplot(as.numeric(Ynorm[tmp_ix, ])~trueclass, main = tmp_gene, xlab="", ylab="", las=2)
}
```

# Using FEAST for scRNA-seq clustering analysis
FEAST requires a count expression matrix with rows representing the genes and columns representing the cells. To preprocess the count matrix, FEAST filters out the genes depending on the dropout rates (sometimes known as zero-expression rates). Alternatively, the users can input the normalized matrix (usually log transformed of relative abundance expression matrix (scaled by average sequencing depth)). Here, we use one template `Yan` [@yan2013single] dataset. It is important to note that FEAST requires the number of clusters (**k**) as an input, which can be obtained by the prior knowledge of the cell types. 

## Load the data
The `Yan` dataset is loaded with two objects including the count expression matrix (Y) and corresponding cell type information (trueclass). Other sample datasets can be downloaded at https://drive.google.com/drive/u/0/folders/1SRT7mrX7ziJoSjuFLLkK8kjnUsJrabVM. If the users want to use one Deng dataset[@deng2014single], the users can load the data and access the phenotype information by using `colData(Deng)` function, and the count expression matrix by using `assay(Deng,"counts")` function from the SummarizedExperiment Bioconductor package. 
```{r load_data, eval = TRUE}
data(Yan)
dim(Y)
table(trueclass)
```

## Consensus clustering
To preprocess the count expression matrix (`Y`), FEAST filters out the genes based on the dropout rate (sometimes known as zero-expression rate). This consensus clustering step will output the initial clustering labels based on a modified CSPA algorithm. It will find the cells that tightly clustered together and possibly filter out the cells that are less correlated to the initial clustering centers. 
```{r consensus, eval = TRUE, message=FALSE, results='hide', include = TRUE}
Y = process_Y(Y, thre = 2) # preprocess the data if needed
con_res = Consensus(Y, k=k)
```


## Calculate the gene-level significance
After the consensus clustering step, FEAST will generate the initial clustering labels with confidence. Based on the initial clustering outcomes, FEAST further infers the gene-level significance by using F-statistics followed by the ranking of the full list of genes. 
```{r gene-level, eval = TRUE, message=FALSE, results='hide', include = TRUE}
F_res = cal_F2(Y, con_res$cluster)
ixs = order(F_res$F_scores, decreasing = TRUE) # order the features
```

## Clustering and validation
After ranking the genes by F-statistics, FEAST curates a series of top number of features (genes). Then, FEAST input these top (by default `top = 500, 1000, 2000`) features into some established clustering algorithm such as SC3 [@kiselev2017sc3], which is regarded as the most accurate scRNA-seq clustering algorithm [@duo2018systematic]. It is worth noting that one needs to confirm the **k** to be the same for every iteration. 

After the clustering steps, with the case of unknown cell types, FEAST evaluates the quality of the feature set by computing the average distance between each cell and the clustering centers, which is the same as the MSE. It is noting that FEAST uses all the features (genes) for distance calculation for the purpose of fair comparisons. The clustering centers are obtained by using previous clustering step. 

It is noted that FEAST is not limited to SC3 clustering. It also shows superior performance on other clustering methods such as TSCAN. The code for running SC3 and TSCAN are embedded in FEAST package, which can be access by **SC3_Clust** and **TSCAN_Clust**. 
```{r validation, eval = TRUE, message=FALSE, results='hide', include = TRUE}
## clustering step
tops = c(500, 1000, 2000)
cluster_res = NULL
for (top in tops){
    tmp_ixs = ixs[1:top]
    tmp_markers = rownames(Y)[tmp_ixs]
    tmp_res = TSCAN_Clust(Y, k = k, input_markers = tmp_markers)
    #tmp_res = SC3_Clust(Y, k = k, input_markers = tmp_markers)
    cluster_res[[toString(top)]] = tmp_res
}
## validation step
Ynorm = Norm_Y(Y)
mse_res = NULL
for (top in names(cluster_res)){
    tmp_res = cluster_res[[top]]
    tmp_cluster = tmp_res$cluster
    tmp_mse = cal_MSE(Ynorm = Ynorm, cluster = tmp_cluster)
    mse_res = c(mse_res, tmp_mse)
}
names(mse_res) = names(cluster_res)
```

## Compare to the real cell type labels

### Benchmark with the original SC3
After the validation step, the feature set associated with the smallest MSE value will be recommended for the further analysis. Here, we demonstrate that an improved clustering accuracy by specifying the optimal feature set. The benchmark comparison is the original setting of established clustering algorithm (e.g. SC3). 
```{r demo, eval = TRUE, message=FALSE, results='hide', include = TRUE}
original = TSCAN_Clust(Y, k=k)
id = which.min(mse_res)
eval_Cluster(original$cluster, trueclass)
eval_Cluster(cluster_res[[id]]$cluster, trueclass)
```

### Show the clustering improvement by using figures
As demonstrated in the figure below, the first panel includes the PCA illustration of the cell types. The second panel shows the clustering result by the original SC3, in which some cells (inside the gray circle) are mixed. The third panel shows the improved clustering outcomes by inputing the optimized feature set into SC3. This is an example for the `Deng` dataset. 
![](Deng_PCA.png){ width="600"}


# Quick use step-by-step
The users can run the `Consensus` function and then run the `Select_Model_short_SC3` function; however, `Consensus` step could take a while especially for the large dataset (sample size is greater than 2000). 
```{r wraper-group, eval = TRUE, message=FALSE, results='hide', include = TRUE, fig.height=4.5, fig.width=8}
data(Yan)
Y = process_Y(Y, thre = 2) # preprocess the data if needed
con_res = Consensus(Y, k=k)
mod_res = Select_Model_short_TSCAN(Y, cluster = con_res$cluster, top = c(200, 500, 1000, 2000))
# mod_res = Select_Model_short_SC3(Y, cluster = con_res$cluster, top = c(200, 500, 1000, 2000))
# to visualize the result, one needs to load ggpubr library. 
library(ggpubr)
PP = Visual_Rslt(model_cv_res = mod_res, trueclass = trueclass)
print(PP$ggobj) # show the visualization plot.
```
The result figure shows that the MSE validation process is concordant with the clustering accuracy trend. It demonstrated that using top 1000 featured genes can result in the highest clustering accuracy with improvement than the original setting. 


# Quick use by the wrapper function
To quickly obtain the ranking orders of the features, the users can directly apply `FEAST` function, which works perfectly fun on small dataset. For the large dataset, the users can consider change the dimention reduction method to *irlba* by specifying `dim_reduce="irlba"`. Moreover, the users can resort to `FEAST_fast` function for fast calculation on the large datasets. For extreme large dataset (sample size >5000), it will split the dataset equally into chucks and apply FEAST algorithm in parallel on individual splitted datasets. In this case, the users need to specify two parameters `split = TRUE`, and `batch_size = 1000` corresponding to the splitted size.  
```{r wraperfunctions, eval = TRUE, message=FALSE, results='hide', include = TRUE, fig.height=4.5, fig.width=8}
ixs = FEAST(Y, k=k)
#ixs = FEAST_fast(Y, k=k)
#ixs = FEAST_fast(Y, k=k, split=TRUE, batch_size = 1000)
```


# Session Info
```{r}
sessionInfo()
```

# Reference

