---
title: "SC3 package manual"
author: "Vladimir Kiselev"
date: "`r Sys.Date()`"
output:
    BiocStyle::html_document:
        toc: true
vignette: >
  %\VignetteIndexEntry{SC3 package manual}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---


```{r knitr-options, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
opts_chunk$set(fig.align = 'center', fig.width = 6, fig.height = 5, dev = 'png')
```

# Introduction

Single-Cell Consensus Clustering (`SC3`) is a tool for unsupervised clustering of scRNA-seq data. `SC3` achieves high accuracy and robustness by consistently integrating different clustering solutions through a consensus approach. An interactive graphical implementation makes `SC3` accessible to a wide audience of users. In addition, `SC3` also aids biological interpretation by identifying marker genes, differentially expressed genes and outlier cells. A manuscript describing `SC3` in details is published in [Nature Methods](http://dx.doi.org/10.1038/nmeth.4236).

# Quality Control, Normalisation and `scater`

`SC3` is a purely clustering tool and it does not provide functions for the sequencing quality control (QC) or normalisation. On the contrary it is expected that these preprocessing steps are performed by a user in advance. To encourage the preprocessing, `SC3` is built on top of the Bioconductor’s `scater` package. To our knowledge the `scater` is the most comprehensive toolkit for the QC and normalisation analysis of the single-cell RNA-Seq data.

The basic `scater` data container is an `SCESet` object. `SC3` implements several methods that allow one to perform clustering of the expression data contained in the `SCESet` object. All results of `SC3` calculations are written to the `sc3` slot of the `SCESet` object.

# Quick Start

## `SC3` Input

If you already have an `SCESet` object created and QCed using `scater` then proceed to the next chapter.

If you have a matrix containing expression data that was QCed and normalised by some other tool, then we first need to form an `SCESet` object containing the data. For illustrative purposes we will use an example expression matrix provided with `SC3`. This matrix (`treutein`) represents __FPKM__ gene expression of 80 cells derived from the distal lung epithelium of mice. The authors ([Treutlein et al.](http://www.nature.com/nature/journal/v509/n7500/full/nature13173.htm)) had computationally identified 5 clusters in the data. The rows in the `treutlein` dataset correspond to genes and columns correspond to cells. Column names correspond to clusters identified by the authors.

```{r, message=FALSE, warning=FALSE}
library(scater)
library(SC3)
treutlein[1:3, 1:3]
```

It is easy to create an `SCESet` object from `treutlein` expression matrix. We will follow the [`scater`'s manual](https://www.bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/vignette.html):
```{r}
# cell annotation
ann <- data.frame(cell_type1 = colnames(treutlein))
pd <- new("AnnotatedDataFrame", data = ann)
# cell expression
tmp <- treutlein
colnames(tmp) <- rownames(ann)
# SCESEt object
sceset <- newSCESet(fpkmData = tmp, phenoData = pd, logExprsOffset = 1)
```

It is also essential for `SC3` that the QC metrics is computed for the created object:
```{r}
is_exprs(sceset) <- exprs(sceset) > 0
sceset <- calculateQCMetrics(sceset)
```

The `treutlein_cell_info` dataframe contains just `cell_type1` column which correspond to the cell labels provided by authors of the original publication. Note that in general it can also contain more information about the cells, such as plate, run, well, date etc.

After the `SCESet` object is created and QC is run, `scater` allows a user to quickly visualize and assess the data, for example using a PCA plot:
```{r}
plotPCA(sceset, colour_by = "cell_type1")
```

## Run SC3

If you would like to explore clustering of your data in the range of `k`s (the number of clusters) from 2 to 4, you just need to run the main `sc3` method and define the range of `k`s using the `ks` parameter (here we also ask `SC3` to calculate biological features based on the identified cell clusters):
```{r}
# Note that n_cores = 1 is required for compilation of this vignette.
# Please remove this parameter when running on your computer:
# sceset <- sc3(sceset, ks = 2:4, biology = TRUE)
sceset <- sc3(sceset, ks = 2:4, biology = TRUE, n_cores = 1)
```

To quickly and easily explore the `SC3` solutions using an interactive Shiny application use the following method:
```{r, eval=FALSE}
sc3_interactive(sceset)
```

Visual exploration can provide a reasonable estimate of the number of clusters `k`. Once a preferable `k` is chosen it is also possible to export the results into an Excel file:
```{r eval=FALSE}
sc3_export_results_xls(sceset)
```

This will write all results to `sc3_results.xls` file. The name of the file can be controlled by the `filename` parameter.

## phenoData

`SC3` writes all its results obtained for cells to the `phenoData` slot of the `SCESet` object by adding additional columns to it. This slot also contains all other cell features calculated by the `scater` package either automatically during the `SCESet` object creation or during the `calculateQCMetrics` call. One can identify the `SC3` results using the `"sc3_"` prefix:

```{r}
p_data <- pData(sceset)
head(p_data[ , grep("sc3_", colnames(p_data))])
```

Additionally, having `SC3` results stored in the same slot makes it possible to highlight them during any of the `scater`'s plotting function call, for example:
```{r}
plotPCA(
    sceset, 
    colour_by = "sc3_3_clusters", 
    size_by = "sc3_3_log2_outlier_score"
)
```

## featureData

`SC3` writes all its results obtained for features (genes/transcripts) to the `featureData` slot of the `SCESet` object by adding additional columns to it. This slot also contains all other feature values calculated by the `scater` package either automatically during the `SCESet` object creation or during the `calculateQCMetrics` call. One can identify the `SC3` results using the `"sc3_"` prefix:

```{r}
f_data <- fData(sceset)
head(f_data[ , grep("sc3_", colnames(f_data))])
```

Because the biological features were also calculated for each `k`, one can find ajusted p-values for both differential expression and marker genes, as well as the area under the ROC curve values (see `?sc3_calcl_biology` for more information).

Again, having `SC3` results stored in the same slot makes it possible to highlight them during any of the `scater`'s plotting function call, for example:
```{r}
plotFeatureData(
    sceset, 
    aes(
        x = sc3_3_markers_clusts, 
        y = sc3_3_markers_auroc, 
        colour = sc3_3_markers_padj
    )
)
```

# Number of Сells

The default settings of `SC3` allow to cluster (using a single `k`) a dataset of 2,000 cells in about 20-30 minutes. 

For datasets with more than 2,000 cells `SC3` automatically adjusts some of its parameters (see below). This allows to cluster a dataset of 5,000 cells in about 20-30 minutes. The parameters can also be manually adjusted for datasets with any number of cells.

For datasets with more than 5,000 cells `SC3` utilizes a hybrid approach that combines unsupervised and supervised clusterings (see below). Namely, `SC3` selects a subset of cells uniformly at random, and obtains clusters from this subset. Subsequently, the inferred labels are used to train a Support Vector Machine (SVM), which is employed to assign labels to the remaining cells. Training cells can also be manually selected by providing their indeces.

# Plot Functions

`SC3` also provides methods for plotting all figures from the interactive session.

## Consensus Matrix

The consensus matrix is a <em>N</em> by <em>N</em> matrix, where <em>N</em> is the number of cells in the input dataset. It represents similarity between the cells based on the averaging of clustering results from all  combinations of clustering parameters. Similarity 0 (blue) means that the two cells are always assigned to different clusters. In contrast, similarity 1 (red) means that the two cells are always assigned to the same cluster. The consensus matrix is clustered by hierarchical clustering and has a diagonal-block structure. Intuitively, the perfect  clustering is achieved when all diagonal blocks are completely red and all off-diagonal elements are completely blue.
```{r, fig.height=6}
sc3_plot_consensus(sceset, k = 3)
```

It is also possible to annotate cells (columns of the consensus matrix) with any column of the `phenoData` slot of the `SCESet` object.
```{r, fig.height=6, fig.width=8}
sc3_plot_consensus(
    sceset, k = 3, 
    show_pdata = c(
        "cell_type1", 
        "log10_total_features",
        "sc3_3_clusters", 
        "sc3_3_log2_outlier_score"
    )
)
```

## Silhouette Plot

A silhouette is a quantitative measure of the diagonality of the consensus matrix. An average silhouette width (shown at the bottom left of the silhouette plot) varies from 0 to 1, where 1 represents a perfectly block-diagonal consensus matrix and 0 represents a situation where there is no block-diagonal structure. The best clustering is achieved when the average silhouette width is close to 1.
```{r}
sc3_plot_silhouette(sceset, k = 3)
```

## Expression Matrix

The expression panel represents the original input expression matrix (cells in columns and genes in rows) after cell and gene filters. Genes are clustered by kmeans with k = 100 (dendrogram on the left) and the heatmap represents the expression levels of the gene cluster centers after log2-scaling.
```{r, fig.height=6}
sc3_plot_expression(sceset, k = 3)
```

It is also possible to annotate cells (columns of the expression matrix) with any column of the `phenoData` slot of the `SCESet` object.
```{r, fig.height=6, fig.width=8}
sc3_plot_expression(
    sceset, k = 3, 
    show_pdata = c(
        "cell_type1", 
        "log10_total_features",
        "sc3_3_clusters", 
        "sc3_3_log2_outlier_score"
    )
)
```

## Cluster Stability

Stability index shows how stable each cluster is accross the selected range of `k`s. The stability index varies between 0 and 1, where 1 means that the same cluster appears in every solution for different `k`.
```{r, fig.height=3}
sc3_plot_cluster_stability(sceset, k = 3)
```

## DE genes

Differential expression is calculated using the non-parametric Kruskal-Wallis test. A significant p-value indicates that gene expression in at least one cluster stochastically dominates one other cluster. SC3 provides a list of all differentially expressed genes with adjusted p-values < 0.01 and plots gene expression profiles of the 50 genes with the lowest p-values. Note that the calculation of differential expression after clustering can introduce a bias in the distribution of p-values, and thus we advise to use the p-values for ranking the genes only.
```{r, fig.height=9}
sc3_plot_de_genes(sceset, k = 3)
```

It is also possible to annotate cells (columns of the matrix containing DE genes) with any column of the `phenoData` slot of the `SCESet` object.
```{r, fig.height=9, fig.width=8}
sc3_plot_de_genes(
    sceset, k = 3, 
    show_pdata = c(
        "cell_type1", 
        "log10_total_features",
        "sc3_3_clusters", 
        "sc3_3_log2_outlier_score"
    )
)
```

## Marker Genes

To find marker genes, for each gene a binary classifier is constructed based on the mean cluster expression values. The classifier prediction is then calculated using the gene expression ranks. The area under the receiver operating characteristic (ROC) curve is used to quantify the accuracy of the prediction. A p-value is assigned to each gene by using the Wilcoxon signed rank test. By default the genes with the area under the ROC curve (AUROC) > 0.85 and with the p-value < 0.01 are selected and the top 10 marker genes of each cluster are visualized in this heatmap.
```{r, fig.height=6}
sc3_plot_markers(sceset, k = 3)
```

It is also possible to annotate cells (columns of the matrix containing marker genes) with any column of the `phenoData` slot of the `SCESet` object.
```{r, fig.height=6, fig.width=8}
sc3_plot_markers(
    sceset, k = 3, 
    show_pdata = c(
        "cell_type1", 
        "log10_total_features",
        "sc3_3_clusters", 
        "sc3_3_log2_outlier_score"
    )
)
```

# SC3 in Detail

The main `sc3` method explained above is a wrapper that calls several other `SC3` methods in the following order:

* `sc3_prepare`
* _(optional)_ `sc3_estimate_k`
* `sc3_calc_dists`
* `sc3_calc_transfs`
* `sc3_kmeans`
* `sc3_calc_consens`
* _(optional)_ `sc3_calc_biology`

Let us go through each of them independently.

## `sc3_prepare`

We start with `sc3_prepare`. This method prepares an object of `SCESet` class for `SC3` clustering. This method also defines all parameters needed for clustering and stores them in the `sc3` slot. The parameters have their own defaults but can be manually changed. For more information on the parameters please use `?sc3_prepare`. 

```{r}
# Note that n_cores = 1 is required for compilation of this vignette.
# Please remove this parameter when running on your computer:
# sceset <- sc3_prepare(sceset, ks = 2:4)
sceset <- sc3_prepare(sceset, ks = 2:4, n_cores = 1)
str(sceset@sc3)
```

## _(optional)_ `sc3_estimate_k`

When the `SCESet` object is prepared for clustering, `SC3` can also estimate the optimal number of clusters `k` in the dataset. `SC3` utilizes the Tracy-Widom theory on random matrices to estimate `k`. `sc3_estimate_k` method creates and populates the following items of the `sc3` slot:

* `k_estimation` - contains the estimated value of `k`.
```{r}
sceset <- sc3_estimate_k(sceset)
str(sceset@sc3)
```

## `sc3_calc_dists`

Now we are ready to perform the clustering itself. First `SC3` calculates distances between the cells. Method `sc3_calc_dists` calculates the distances, creates and populates the following items of the `sc3` slot:

* `distances` - contains a list of distance matrices corresponding to Euclidean, Pearson and Spearman distances.
```{r}
sceset <- sc3_calc_dists(sceset)
names(sceset@sc3$distances)
```

## `sc3_calc_transfs`

Next the distance matrices are transformed using PCA and graph Laplacian. Method `sc3_calc_transfs` calculates transforamtions of the distance matrices contained in 
the `distances` item of the `sc3` slot. It then creates and populates the following items of the `sc3` slot:

* `transformations` - contains a list of transformations of the distance matrices corresponding to PCA and graph Laplacian transformations.
```{r}
sceset <- sc3_calc_transfs(sceset)
names(sceset@sc3$transformations)
```

It also removes the previously calculated `distances` item from the `sc3` slot:
```{r}
sceset@sc3$distances
```

## `sc3_kmeans`

kmeans should then be performed on the transformed distance matrices contained in the `transformations` item of the `sc3` slot. Method `sc3_kmeans` creates and populates the following items of the `sc3` slot:

* `kmeans` - contains a list of kmeans clusterings.

By default the `nstart` parameter passed to `kmeans` defined in `sc3_prepare` method, is set 1000 and written to `kmeans_nstart` item of the `sc3` slot. If the number of cells in the dataset is more than 2,000, this parameter is set to 50. A user can also manually define this parameter by changing the value of the `kmeans_nstart` item of the `sc3` slot.
```{r}
sceset <- sc3_kmeans(sceset)
names(sceset@sc3$kmeans)
```

## `sc3_calc_consens`

In this step `SC3` will provide you with a clustering solution. Let's first check that there are no `SC3` related columns in the `phenoData` slot:
```{r}
p_data <- pData(sceset)
head(p_data[ , grep("sc3_", colnames(p_data))])
```

When calculating consensus for each value of `k` `SC3` averages the clustering results of `kmeans` using a consensus approach. Method `sc3_calc_consens` calculates consensus matrices based on the clustering solutions contained in the `kmeans` item of the `sc3` slot. It then creates and populates the following items of the `sc3` slot:

* `consensus` - for each value of `k` it contains: a consensus matrix, an `hclust` object, corresponding to hierarchical clustering of the consensus matrix and the Silhouette indeces of the clusters.
```{r}
sceset <- sc3_calc_consens(sceset)
names(sceset@sc3$consensus)
names(sceset@sc3$consensus$`3`)
```

It also removes the previously calculated `kmeans` item from the `sc3` slot:
```{r}
sceset@sc3$kmeans
```

As mentioned before all the clustering results (cell-related information) are written to the `phenoData` slot of the `SCESet` object:
```{r}
p_data <- pData(sceset)
head(p_data[ , grep("sc3_", colnames(p_data))])
```

We can see that `SC3` calculated clusters for `k = 2, 3` and `4` and wrote them to the `phenoData` slot of the `SCESet` object.

## _(optional)_ `sc3_calc_biology`

`SC3` can also calculates DE genes, marker genes and cell outliers based on the calculated consensus clusterings. Similary to the clustering solutions, method `sc3_calc_biology` writes the results for the cell outliers (cell-related information) to the `phenoData` slot of the `SCESet` object. In contrast, DE and marker genes results (gene-related information) is are written to the `featureData` slot. In addition `biology` item of the `sc3` slot is set to `TRUE`.
```{r}
sceset <- sc3_calc_biology(sceset)
```

### Cell Outliers

Now we can see that cell outlier scores have been calculated for each value of `k`:
```{r}
p_data <- pData(sceset)
head(p_data[ , grep("sc3_", colnames(p_data))])
```

For more information on how the cell outliers are calculated please see `?get_outl_cells`.

### DE and marker genes

We can also see that DE and marker genes characteristics (adjusted p-values and area under the ROC curve) have been calculated for each value of `k`
```{r}
f_data <- fData(sceset)
head(f_data[ , grep("sc3_", colnames(f_data))])
```

For more information on how the DE and marker genes are calculated please see `?get_de_genes` and `?get_marker_genes`.

# Hybrid `SVM` Approach

For datasets with more than 5,000 cells `SC3` automatically utilizes a hybrid approach that combines unsupervised and supervised clusterings. Namely, `SC3` selects a subset of cells uniformly at random (5,000), and obtains clusters from this subset. The inferred labels can be used to train a Support Vector Machine (`SVM`), which is employed to assign labels to the remaining cells.

The hybrid approach can also be triggered by defining either the `svm_num_cells` parameter (the number of training cells, which is different from 5,000) or `svm_train_inds` parameter (training cells are manually selected by providing their indexes).

Let us first save the `SC3` results for `k = 3` obtained without using the hybrid approach:
```{r}
no_svm_labels <- pData(sceset)$sc3_3_clusters
```

Now let us trigger the hybrid approach by asking for 50 training cells:
```{r}
# Note that n_cores = 1 is required for compilation of this vignette.
# Please remove this parameter when running on your computer:
# sceset <- sc3(sceset, ks = 2:4, svm.num.cells = 50)
sceset <- sc3(sceset, ks = 2:4, biology = TRUE, svm_num_cells = 50, n_cores = 1)
```

Note that when `SVM` is used all results (including marker genes, DE genes and cell outliers) correspond to the training cells only (50 cells), and values of all other cells are set to `NA`:
```{r}
p_data <- pData(sceset)
head(p_data[ , grep("sc3_", colnames(p_data))])
```

Now we can run the `SVM` and predict labels of all the other cells:
```{r, message=FALSE, warning=FALSE}
sceset <- sc3_run_svm(sceset)
p_data <- pData(sceset)
head(p_data[ , grep("sc3_", colnames(p_data))])
```

Note that the cell outlier scores (and also DE and marker genes values) were not updated and they still contain `NA` values for non-training cells. To recalculate biological characteristics using the labels predicted by `SVM` one need to clear the `svm_train_inds` item in the `sc3` slot and rerun the `sc3_calc_biology` method:
```{r}
sceset@sc3$svm_train_inds <- NULL
sceset <- sc3_calc_biology(sceset)
p_data <- pData(sceset)
head(p_data[ , grep("sc3_", colnames(p_data))])
```

Now the biological characteristics are calculated for all cells (including those predicted by the `SVM`)

```{r}
svm_labels <- pData(sceset)$sc3_3_clusters
```

Now we can compare the labels using the adjusted rand index (`ARI`):
```{r}
if (require("mclust")) {
  adjustedRandIndex(no_svm_labels, svm_labels)
}
```

`ARI` is less than `1`, which means that `SVM` results are different from the non-`SVM` results, however `ARI` is still pretty close to `1` meaning that the solutions are very similar.