18-clustering.Rmd

---
output: html_document
---

# Clustering example {#clust-methods}

```{r, echo=FALSE}
library(knitr)
opts_chunk$set(fig.align = "center")
```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
library(pcaMethods)
library(pcaReduce)
library(SC3)
library(scater)
library(pheatmap)
set.seed(1234567)
```

To illustrate clustering of scRNA-seq data, we consider the `Pollen` dataset of cells from 
different human tissues [@Pollen2014-cu]. We have preprocessed the dataset and created a 
scater object in advance. We have also annotated the cells with the cell type information 
(it is the `cell_type1` column in the `phenoData` slot).

## Pollen dataset

Let's load the data and look at it:
```{r}
pollen <- readRDS("pollen/pollen.rds")
pollen
```

Let's look at the cell type annotation:
```{r}
table(pData(pollen)$cell_type1)
```

A simple PCA analysis already separates some strong cell types and provides some insights in the data structure:
```{r}
plotPCA(pollen, colour_by = "cell_type1")
```

## SC3

Let's run `SC3` clustering on the Pollen data. The advantage of the `SC3` is that it can directly take a [scater](http://bioconductor.org/packages/scater/) object (see previous chapters) as an input.

Now let's image we do not know the number of clusters _k_ (cell types). `SC3` can estimate a number of clusters for you:
```{r}
pollen <- sc3_prepare(pollen, ks = 2:5)
pollen <- sc3_estimate_k(pollen)
pollen@sc3$k_estimation
```

Interestingly, the number of cell types predicted by `SC3` is the same as the number of cell types in the Pollen data annotation.

Now we are ready to run `SC3` (we also ask it to calculate biological properties of the clusters): 
```{r}
pollen <- sc3(pollen, ks = 11, biology = TRUE)
```

`SC3` result consists of several different outputs (please look in [@Kiselev2016-bq] and [SC3 vignette](http://bioconductor.org/packages/release/bioc/vignettes/SC3/inst/doc/my-vignette.html) for more details). Here we show some of them:

Consensus matrix:
```{r, fig.height=6}
sc3_plot_consensus(pollen, k = 11, show_pdata = "cell_type1")
```

Silhouette plot:
```{r, fig.height=9}
sc3_plot_silhouette(pollen, k = 11)
```

Heatmap of the expression matrix:
```{r, fig.height=6}
sc3_plot_expression(pollen, k = 11, show_pdata = "cell_type1")
```

Identified marker genes:
```{r, fig.height=17}
sc3_plot_markers(pollen, k = 11, show_pdata = "cell_type1")
```

PCA plot with highlighted `SC3` clusters:
```{r}
plotPCA(pollen, colour_by = "sc3_11_clusters")
```

Note, that one can also run `SC3` in an interactive `Shiny` session:
```{r, eval=FALSE}
sc3_interactive(pollen)
```

This command will open `SC3` in a web browser.

* __Exercise 1__: Run `SC3` for $k$ from 9 to 13 and explore different clustering solutions in your web browser.

* __Exercise 2__: Which clusters are the most stable when $k$ is changed from 9 to 13? (Look at the "Stability" tab)

* __Exercise 3__: Check out differentially expressed genes and marker genes for the obtained clusterings. Please use $k=11$.

* __Exercise 4__: Change the marker genes threshold (the default is 0.85). Does __SC3__ find more marker genes?

## pcaReduce

`pcaReduce` operates directly on the expression matrix. It is recommended to use a gene filter and log transformation before running `pcaReduce`. We will use the default `SC3` gene filter (note that the `exprs` slot of a `scater` object is log-transformed by default).

```{r}
# use the same gene filter as in SC3
input <- exprs(pollen[fData(pollen)$sc3_gene_filter, ])
```

There are several parameters used by `pcaReduce`:
* `nbt` defines a number of `pcaReduce` runs (it is stochastic and may have different solutions after different runs)
* `q` defines number of dimensions to start clustering with. The output will contain partitions for all $k$ from 2 to q+1.
* `method` defines a method used for clustering. `S` - to perform sampling based merging, `M` - to perform merging based on largest probability.

We will run `pcaReduce` 1 time:
```{r}
# run pcaReduce 1 time creating hierarchies from 1 to 30 clusters
pca.red <- PCAreduce(t(input), nbt = 1, q = 30, method = 'S')[[1]]
```

```{r}
pData(pollen)$pcaReduce <- as.character(pca.red[,32 - 11])
plotPCA(pollen, colour_by = "pcaReduce")
```

__Exercise 5__: Run pcaReduce for $k=2$ and plot a similar PCA plot. Does it look good?

__Hint__: When running pcaReduce for different $k$s you do not need to rerun PCAreduce function, just use already calculated `pca.red` object.

__Our solution__:
```{r clust-pca-reduce2, fig.cap = "Clustering solutions of pcaReduce method for $k=2$.", echo=FALSE}
pData(pollen)$pcaReduce <- as.character(pca.red[,32 - 2])
plotPCA(pollen, colour_by = "pcaReduce")
```

__Exercise 6__: Compare the results between `SC3` and `pcaReduce` for $k=11$. What is
the main difference between the solutions provided by the two
different methods?

__Our solution__:
```{r, echo=FALSE}
pData(pollen)$pcaReduce <- as.character(pca.red[,32 - 11])
sc3_plot_expression(pollen, k = 11, show_pdata = "pcaReduce")
```


## tSNE + kmeans

[tSNE](https://lvdmaaten.github.io/tsne/) plots that we saw before (\@ref(visual-tsne)) when used the __scater__ package are made by using the [Rtsne](https://cran.r-project.org/web/packages/Rtsne/index.html) and [ggplot2](https://cran.r-project.org/web/packages/ggplot2/index.html) packages. Here we will do the same:
```{r clust-tsne, fig.cap = "tSNE map of the patient data"}
pollen <- plotTSNE(pollen, rand_seed = 1, return_SCESet = TRUE)
```

Note that all points on the plot above are black. This is different from what we saw before, when the cells were coloured based on the annotation. Here we do not have any annotation and all cells come from the same batch, therefore all dots are black.

Now we are going to apply _k_-means clustering algorithm to the cloud of points on the tSNE map. How many groups do you see in the cloud?

We will start with $k=8$:
```{r clust-tsne-kmeans2, fig.cap = "tSNE map of the patient data with 8 colored clusters, identified by the k-means clustering algorithm"}
pData(pollen)$tSNE_kmeans <- as.character(kmeans(pollen@reducedDimension, centers = 8)$clust)
plotTSNE(pollen, rand_seed = 1, colour_by = "tSNE_kmeans")
```

__Exercise 7__: Make the same plot for $k=11$.

__Exercise 8__: Compare the results between `SC3` and `tSNE+kmeans`. Can the
results be improved by changing the `perplexity` parameter?

__Our solution__:
```{r, echo=FALSE}
pData(pollen)$tSNE_kmeans <- as.character(kmeans(pollen@reducedDimension, centers = 11)$clust)
sc3_plot_expression(pollen, k = 11, show_pdata = "tSNE_kmeans")
```

As you may have noticed, both `pcaReduce` and `tSNE+kmeans` are stochastic
and give different results every time they are run. To get a better
overview of the solutions, we need to run the methods multiple times. `SC3` is also stochastic, but thanks to the consensus step, it is more robust and less likely to produce different outcomes.

## SNN-Cliq

Here we run SNN-cliq with te default parameters provided in the author's example:

```{r}
distan <- "euclidean"
par.k <- 3
par.r <- 0.7
par.m <- 0.5
# construct a graph
scRNA.seq.funcs::SNN(
    data = t(input),
    outfile = "snn-cliq.txt",
    k = par.k,
    distance = distan
)
# find clusters in the graph
snn.res <- 
    system(
        paste0(
            "python snn-cliq/Cliq.py ", 
            "-i snn-cliq.txt ",
            "-o res-snn-cliq.txt ",
            "-r ", par.r,
            " -m ", par.m
        ),
        intern = TRUE
    )
cat(paste(snn.res, collapse = "\n"))
snn.res <- read.table("res-snn-cliq.txt")
# remove files that were created during the analysis
system("rm snn-cliq.txt res-snn-cliq.txt")

pData(pollen)$SNNCliq <- as.character(snn.res[,1])
plotPCA(pollen, colour_by = "SNNCliq")
```

__Exercise 9__: Compare the results between `SC3` and `SNN-Cliq`.

__Our solution__:
```{r, echo=FALSE}
sc3_plot_expression(pollen, k = 11, show_pdata = "SNNCliq")
```

## SINCERA

As mentioned in the previous chapter [SINCERA](https://research.cchmc.org/pbge/sincera.html) is based on hierarchical clustering. One important thing to keep in mind is that it performs a gene-level z-score transformation before doing clustering:

```{r, echo=TRUE, fig.height=7, fig.width=7}
# perform gene-by-gene per-sample z-score transformation
dat <- apply(input, 1, function(y) scRNA.seq.funcs::z.transform.helper(y))
# hierarchical clustering
dd <- as.dist((1 - cor(t(dat), method = "pearson"))/2)
hc <- hclust(dd, method = "average")
```

If the number of cluster is not known [SINCERA](https://research.cchmc.org/pbge/sincera.html) can identify __k__ as the minimum height of the hierarchical tree that generates no more than a specified number of singleton clusters (clusters containing only 1 cell)
```{r, echo=TRUE}
num.singleton <- 0
kk <- 1
for (i in 2:dim(dat)[2]) {
    clusters <- cutree(hc, k = i)
    clustersizes <- as.data.frame(table(clusters))
    singleton.clusters <- which(clustersizes$Freq < 2)
    if (length(singleton.clusters) <= num.singleton) {
        kk <- i
    } else {
        break;
    }
}
cat(kk)
```

Let's now visualize the SINCERA results as a heatmap:
```{r clust-sincera, fig.cap = "Clustering solutions of SINCERA method using $k=3$"}
pheatmap(
    t(dat),
    cluster_cols = hc,
    cutree_cols = 14,
    kmeans_k = 100,
    show_rownames = FALSE
)
```

__Exercise 10__: Compare the results between `SC3` and `SNN-Cliq`.

__Our solution__:
```{r, echo=FALSE}
pData(pollen)$SINCERA <- as.character(cutree(hc, k = i))
sc3_plot_expression(pollen, k = 11, show_pdata = "SINCERA")
```

__Exercise 11__: Is using the singleton cluster criteria for finding __k__ a good idea?

## SEURAT

Here we follow an [example](http://satijalab.org/seurat/get_started.html) created by the authors of `SEURAT` (8,500 Pancreas cells). We mostly use default values in various function calls, for more details please consult the documentation and the authors:

```{r, message=FALSE, warning=FALSE, eval=FALSE}
library(Seurat)
pollen_seurat <- new("seurat", raw.data = get_exprs(pollen, exprs_values = "tpm"))
pollen_seurat <- Setup(pollen_seurat, project = "Pollen")
pollen_seurat <- MeanVarPlot(pollen_seurat)
pollen_seurat <- RegressOut(pollen_seurat, latent.vars = c("nUMI"), 
                            genes.regress = pollen_seurat@var.genes)
pollen_seurat <- PCAFast(pollen_seurat)
pollen_seurat <- RunTSNE(pollen_seurat)
pollen_seurat <- FindClusters(pollen_seurat)
TSNEPlot(pollen_seurat, do.label = T)
```

__Exercise 12__: Compare the results between `SC3` and `SEURAT`.

__Our solution__:
```{r, echo=TRUE, eval=FALSE}
pData(pollen)$SEURAT <- as.character(pollen_seurat@ident)
sc3_plot_expression(pollen, k = 11, show_pdata = "SEURAT")
```


Seurat can also find marker genes, e.g. marker genes for cluster 2:
```{r, eval=FALSE}
markers <- FindMarkers(pollen_seurat, 2)
FeaturePlot(pollen_seurat, 
            head(rownames(markers)), 
            cols.use = c("lightgrey", "blue"), 
            nCol = 3)
```

__Exercise 13__: Compare marker genes provided by `SEURAT` and `SC3`.