---
title: Cluster Analysis in R 
author: "First/last name (first.last@ucr.edu)"
date: "Last update: `r format(Sys.time(), '%d %B, %Y')`" 
output:
  html_document:
    toc: true
    toc_float:
        collapsed: true
        smooth_scroll: true
    toc_depth: 3
    fig_caption: yes
    code_folding: show
    number_sections: true

fontsize: 14pt
bibliography: bibtex.bib
weight: 13
type: docs
---

<!--
- Compile from command-line
Rscript -e "rmarkdown::render('rclustering.Rmd', c('html_document'), clean=F); knitr::knit('rclustering.Rmd', tangle=TRUE)"
-->

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
options(width=100, max.print=1000)
knitr::opts_chunk$set(
    eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")),
    cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE")))
```

```{r setup, echo=FALSE, messages=FALSE, warnings=FALSE}
suppressPackageStartupMessages({
    library(ggplot2)
    library(ape)
})
```

<div style="text-align: right"> 
Source code downloads: &nbsp; &nbsp;
[ [.Rmd](https://raw.githubusercontent.com/tgirke/GEN242//main/content/en/tutorials/rclustering/rclustering.Rmd) ] &nbsp; &nbsp; 
[ [.R](https://raw.githubusercontent.com/tgirke/GEN242//main/content/en/tutorials/rclustering/rclustering.R) ]
</div>


## Introduction

- What is Clustering?
    - Clustering is the classification of data objects into similarity groups (clusters) according to a defined distance measure. 
    - It is used in many fields, such as machine learning, data mining, pattern recognition, image analysis, genomics, systems biology, etc. 
    - Machine learning typically regards data clustering as a form of unsupervised learning.

- Why Clustering and Data Mining in R?}
    - Efficient data structures and functions for clustering
    - Reproducible and programmable
    - Comprehensive set of clustering and machine learning libraries 
    - Integration with many other data analysis tools

- Useful Links
    - [Cluster Task Views](https://cran.r-project.org/web/views/Cluster.html)
    - [Machine Learning Task Views](https://cran.r-project.org/web/views/MachineLearning.html)
    - [UCR Manual](http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual#TOC-Clustering-and-Data-Mining-in-R)

## Data Preprocessing

### Data Transformations

Choice depends on data set!
	
- Center and standardize
    1. Center: subtract from each value the mean of the corresponding vector
	2.  Standardize: devide by standard deviation
	- Result: _Mean = 0_ and _STDEV = 1_

- Center and scale with the `scale()` function
    1. Center: subtract from each value the mean of the corresponding vector
	2. Scale: divide centered vector by their _root mean square_ (_rms_):
    $$ x_{rms} = \sqrt[]{\frac{1}{n-1}\sum_{i=1}^{n}{x_{i}{^2}}} $$
    - Result: _Mean = 0_ and _STDEV = 1_

- Log transformation 
- Rank transformation: replace measured values by ranks 
- No transformation

### Distance Methods

List of most common ones!

- Euclidean distance for two profiles _X_ and _Y_:
  $$ d(X,Y) = \sqrt[]{ \sum_{i=1}^{n}{(x_{i}-y_{i})^2} }$$
    - __Disadvantages__: not scale invariant, not for negative correlations
- Maximum, Manhattan, Canberra, binary, Minowski, ...
- Correlation-based distance: _1-r_
    - Pearson correlation coefficient (PCC):
	  $$r = \frac{n\sum_{i=1}^{n}{x_{i}y_{i}} - \sum_{i=1}^{n}{x_{i}} \sum_{i=1}^{n}{y_{i}}}{ \sqrt[]{(\sum_{i=1}^{n}{x_{i}^2} - (\sum_{i=1}^{n}{x_{i})^2}) (\sum_{i=1}^{n}{y_{i}^2} - (\sum_{i=1}^{n}{y_{i})^2})} }$$
        - __Disadvantage__: outlier sensitive 
	- Spearman correlation coefficient (SCC)
	    - Same calculation as PCC but with ranked values!

There are many more distance measures

- If the distances among items are quantifiable, then clustering is possible.
- Choose the most accurate and meaningful distance measure for a given field of application.
- If uncertain then choose several distance measures and compare the results. 

### Cluster Linkage

<center><img title="cluster linkage" src="../images/linkage.png" ></center>

## Clustering Algorithms

### Hierarchical Clustering

#### Overview of algorithm 

1. Identify clusters (items) with closest distance  
2. Join them to new clusters
3. Compute distance between clusters (items)
4. Return to step 1

##### Hierarchical clustering: agglomerative Approach
<center><img title="hierarchical clustering" src="../images/hierarchical.png" ></center>

##### Hierarchical Clustering with Heatmap
<center><img title="heatmap" src="../images/heatmap.png" ></center>

- A heatmap is a color coded table. To visually identify patterns, the rows and columns of a heatmap are often sorted by hierarchical clustering trees.  
- In case of gene expression data, the row tree usually represents the genes, the column tree the treatments and the colors in the heat table represent the intensities or ratios of the underlying gene expression data set.

#### Hierarchical Clustering Approaches

1. Agglomerative approach (bottom-up)
    - R functions: `hclust()` and `agnes()`

2. Divisive approach (top-down)
    - R function: `diana()`

#### Tree Cutting to Obtain Discrete Clusters

1. Node height in tree
2. Number of clusters
3. Search tree nodes by distance cutoff


#### Examples

##### Using `hclust` and `heatmap.2`

Note, with large data sets consider using [flashClust](https://cran.r-project.org/web/packages/flashClust/index.html) which is a fast implementation of hierarchical clustering. 

```{r hclust_heatmap_example, eval=TRUE, warning=FALSE, message=FALSE}
library(gplots) 
y <- matrix(rnorm(500), 100, 5, dimnames=list(paste("g", 1:100, sep=""), paste("t", 1:5, sep=""))) 
heatmap.2(y) # Shortcut to final result
```

##### Stepwise Approach with Tree Cutting

```{r hclust_heatmap_example_setpwise, eval=TRUE, warning=FALSE, message=FALSE}
## Row- and column-wise clustering 
hr <- hclust(as.dist(1-cor(t(y), method="pearson")), method="complete")
hc <- hclust(as.dist(1-cor(y, method="spearman")), method="complete") 
## Tree cutting
mycl <- cutree(hr, h=max(hr$height)/1.5); mycolhc <- rainbow(length(unique(mycl)), start=0.1, end=0.9); mycolhc <- mycolhc[as.vector(mycl)] 
## Plot heatmap 
mycol <- colorpanel(40, "darkblue", "yellow", "white") # or try redgreen(75)
heatmap.2(y, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=mycol, scale="row", density.info="none", trace="none", RowSideColors=mycolhc) 
```

### K-Means Clustering

#### Overview of algorithm 

1. Choose the number of k clusters   
2. Randomly assign items to the k clusters
3. Calculate new centroid for each of the k clusters
4. Calculate the distance of all items to the k centroids
5. Assign items to closest centroid
6. Repeat until clusters assignments are stable

<center><img title="kmeans" src="../images/kmeans.png" ></center>
	
#### Examples

```{r kmeans_example, eval=TRUE, warning=FALSE, message=FALSE}
km <- kmeans(t(scale(t(y))), 3)
km$cluster 
```

### Fuzzy C-Means Clustering

- In contrast to strict (hard) clustering approaches, fuzzy (soft) clustering methods allow multiple cluster memberships of the clustered items [@Hathaway1996-hu]. 
- This is commonly achieved by assigning to each item a weight of belonging to each cluster. 
- Thus, items at the edge of a cluster, may be in a cluster to a lesser degree than items at the center of a cluster. 
- Typically, each item has as many coefficients (weights) as there are clusters that sum up for each item to one.

#### Examples

##### Fuzzy Clustering with `fanny`

```{r fanny_example, eval=TRUE, warning=FALSE, message=FALSE}
library(cluster) # Loads the cluster library.
fannyy <- fanny(y, k=4, metric = "euclidean", memb.exp = 1.2)
round(fannyy$membership, 2)[1:4,]
fannyy$clustering # Hard clustering result
(fannyyMA <- fannyy$membership > 0.20)[1:4,] # Soft clustering result
apply(fannyyMA, 1, which)[1:4] # First 4 clusters
```
	
### Principal Component Analysis (PCA)

Principal components analysis (PCA) is a data reduction technique that allows to simplify multidimensional data sets to 2 or 3 dimensions for plotting purposes and visual variance analysis.

#### Basic Steps

- Center (and standardize) data
- First principal component axis
    - Across centroid of data cloud
	- Distance of each point to that line is minimized, so that it crosses the maximum variation of the data cloud
- Second principal component axis 
    - Orthogonal to first principal component
	- Along maximum variation in the data
- First PCA axis becomes x-axis and second PCA axis y-axis 
- Continue process until the necessary number of principal components is obtained 

<center><img title="pca" src="../images/pca.png" ></center>

#### Example

```{r pca_example, eval=TRUE, warning=FALSE, message=FALSE}
pca <- prcomp(y, scale=T)
summary(pca) # Prints variance summary for all principal components
plot(pca$x, pch=20, col="blue", type="n") # To plot dots, drop type="n"
text(pca$x, rownames(pca$x), cex=0.8)
```
1st and 2nd principal components explain x% of variance in data.

### Multidimensional Scaling (MDS)

- Alternative dimensionality reduction approach
- Represents distances in 2D or 3D space
- Starts from distance matrix (PCA uses data points)

#### Example

The following example performs MDS analysis with `cmdscale` on the geographic distances among European cities.

```{r mds_example, eval=TRUE, warning=FALSE, message=FALSE}
loc <- cmdscale(eurodist) 
plot(loc[,1], -loc[,2], type="n", xlab="", ylab="", main="cmdscale(eurodist)")
text(loc[,1], -loc[,2], rownames(loc), cex=0.8) 
```

### Biclustering

Finds in matrix subgroups of rows and columns which are as similar as possible to each other and as different as possible to the remaining data points.

<center><img title="pca" src="../images/biclust.png"  ></center>
<center> Unclustered --------------------------> Clustered</center>

### Similarity Measures for Clusters

- Compare the numbers of identical and unique item pairs appearing in cluster sets
- Achieved by counting the number of item pairs found in both clustering sets _(a)_ as well as the pairs appearing only in the first _(b)_ or the second _(c)_ set. 
- With this a similarity coefficient, such as the Jaccard index, can be computed. The latter is defined as the size of the intersect divided by the size of the union of two sample sets: _a/(a+b+c)_. 
- In case of partitioning results, the Jaccard Index measures how frequently pairs of items are joined together in two clustering data sets and how often pairs are observed only in one set. 
- Related coefficient are the Rand Index and the Adjusted Rand Index. These indices also consider the number of pairs _(d)_ that are not joined together in any of the clusters in both sets. 

#### Example: 

##### Jaccard index for cluster sets

The following imports the `cindex()` function and computes the Jaccard Index for two sample cluster results.

```{r jaccard_index, eval=TRUE, warning=FALSE, message=FALSE}
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/clusterIndex.R") 
library(cluster); y <- matrix(rnorm(5000), 1000, 5, dimnames=list(paste("g", 1:1000, sep=""), paste("t", 1:5, sep=""))); clarax <- clara(y, 49); clV1 <- clarax$clustering; clarax <- clara(y, 50); clV2 <- clarax$clustering 
ci <- cindex(clV1=clV1, clV2=clV2, self=FALSE, minSZ=1, method="jaccard")
ci[2:3] # Returns Jaccard index and variables used to compute it 
```

##### Clustering cluster sets with Jaccard index

The following example shows how one can cluster entire cluster result sets. First, 10 sample cluster results are created with Clara using k-values from 3 to 12. The results are stored as named clustering vectors in a list object. Then a nested sapply loop is used to generate a similarity matrix of Jaccard Indices for the clustering results. After converting the result into a distance matrix, hierarchical clustering is performed with `hclust`.}

```{r jaccard_index_clustering, eval=TRUE, warning=FALSE, message=FALSE}
clVlist <- lapply(3:12, function(x) clara(y[1:30, ], k=x)$clustering); names(clVlist) <- paste("k", "=", 3:12)
d <- sapply(names(clVlist), function(x) sapply(names(clVlist), function(y) cindex(clV1=clVlist[[y]], clV2=clVlist[[x]], method="jaccard")[[3]]))
hv <- hclust(as.dist(1-d))
plot(as.dendrogram(hv), edgePar=list(col=3, lwd=4), horiz=T, main="Similarities of 10 Clara Clustering Results for k: 3-12") 
```

- Remember: there are many additional clustering algorithms.
- Additional details can be found in the Clustering Section of the [R/Bioconductor Manual](http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual\#TOC-Clustering-and-Data-Mining-in-R).

## Clustering Exercises

### Data Preprocessing

#### Scaling 

```{r scaling, eval=TRUE, warning=FALSE, message=FALSE}
## Sample data set
set.seed(1410)
y <- matrix(rnorm(50), 10, 5, dimnames=list(paste("g", 1:10, sep=""), 
            paste("t", 1:5, sep="")))
dim(y)
## Scaling
yscaled <- t(scale(t(y))) # Centers and scales y row-wise
apply(yscaled, 1, sd)
```

### Distance Matrices

#### Euclidean distance matrix

```{r euclidean_distance_matrices, eval=TRUE, warning=FALSE, message=FALSE}
dist(y[1:4,], method = "euclidean")
```

#### Correlation-based distance matrix

Correlation matrix
```{r correlation_similarity_matrices, eval=TRUE, warning=FALSE, message=FALSE}
c <- cor(t(y), method="pearson") 
as.matrix(c)[1:4,1:4]
```

Correlation-based distance matrix
```{r correlation_distance_matrices, eval=TRUE, warning=FALSE, message=FALSE}
d <- as.dist(1-c)
as.matrix(d)[1:4,1:4]
```
	
### Hierarchical Clustering with `hclust`

Hierarchical clustering with complete linkage and basic tree plotting
```{r hclust1, eval=TRUE, warning=FALSE, message=FALSE}
hr <- hclust(d, method = "complete", members=NULL)
names(hr)
par(mfrow = c(1, 2)); plot(hr, hang = 0.1); plot(hr, hang = -1) 
```

#### Tree plotting I

```{r hclust_plot_tree1, eval=TRUE, warning=FALSE, message=FALSE}
plot(as.dendrogram(hr), edgePar=list(col=3, lwd=4), horiz=T) 
```

#### Tree plotting II

The `ape` library provides more advanced features for tree plotting

```{r hclust_plot_tree2, eval=TRUE, warning=FALSE, message=FALSE}
library(ape) 
plot.phylo(as.phylo(hr), type="p", edge.col=4, edge.width=2, 
           show.node.label=TRUE, no.margin=TRUE)
```

### Tree Cutting

Accessing information in _hclust_ objects 

```{r hclust_object, eval=TRUE, warning=FALSE, message=FALSE}
hr
## Print row labels in the order they appear in the tree
hr$labels[hr$order] 
```

Tree cutting with `cutree`
```{r cutree, eval=TRUE, warning=FALSE, message=FALSE}
mycl <- cutree(hr, h=max(hr$height)/2)
mycl[hr$labels[hr$order]] 
```

### Heatmaps

#### With `heatmap.2`

All in one step: clustering and heatmap plotting
```{r heatmap2a, eval=TRUE, warning=FALSE, message=FALSE}
library(gplots)
heatmap.2(y, col=redgreen(75))
```

#### With `pheatmap`

All in one step: clustering and heatmap plotting
```{r pheatmap, eval=TRUE, warning=FALSE, message=FALSE}
library(pheatmap); library("RColorBrewer")
pheatmap(y, color=brewer.pal(9,"Blues"))
```

#### Customizing heatmaps

Customizes row and column clustering and shows tree cutting result in row color bar. Additional color schemes can be found [here](http://manuals.bioinformatics.ucr.edu/home/R\_BioCondManual\#clustering\_hc\_submap).

```{r heatmap2_custom, eval=TRUE, warning=FALSE, message=FALSE}
hc <- hclust(as.dist(1-cor(y, method="spearman")), method="complete")
mycol <- colorpanel(40, "darkblue", "yellow", "white")
heatmap.2(y, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=mycol,
          scale="row", density.info="none", trace="none", 
          RowSideColors=as.character(mycl))
```

#### Complex heatmaps

For plotting complex heatmaps, the `ComplexHeatmap` provides useful functionalities. The following creates a sample heatmap. For additional examples, visit the manual [here](https://jokergoo.github.io/ComplexHeatmap-reference/book/).

```{r complex_heatmap, eval=TRUE, warning=FALSE, message=FALSE}
library(ComplexHeatmap)
set.seed(123)
mat = matrix(rnorm(100), 10)
rownames(mat) = paste0("R", 1:10)
colnames(mat) = paste0("C", 1:10)
column_ha = HeatmapAnnotation(foo1 = runif(10), bar1 = anno_barplot(runif(10)))
row_ha = rowAnnotation(foo2 = runif(10), bar2 = anno_barplot(runif(10)))
Heatmap(mat, name = "mat", top_annotation = column_ha, right_annotation = row_ha)
```

### K-Means Clustering with PAM

Runs K-means clustering with PAM (partitioning around medoids) algorithm and shows result in color bar of hierarchical clustering result from before.

```{r kmeans2, eval=TRUE, warning=FALSE, message=FALSE}
library(cluster)
pamy <- pam(d, 4)
(kmcol <- pamy$clustering)
heatmap.2(y, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=mycol,
          scale="row", density.info="none", trace="none", 
          RowSideColors=as.character(kmcol))
```

### K-Means Fuzzy Clustering

Performs k-means fuzzy clustering
```{r kmeans_fuzzy, eval=TRUE, warning=FALSE, message=FALSE}
library(cluster)
fannyy <- fanny(d, k=4, memb.exp = 1.5)
round(fannyy$membership, 2)[1:4,]
fannyy$clustering 
## Returns multiple cluster memberships for coefficient above a certain 
## value (here >0.1)
fannyyMA <- round(fannyy$membership, 2) > 0.10 
apply(fannyyMA, 1, function(x) paste(which(x), collapse="_"))
```

### Multidimensional Scaling (MDS) 

Performs MDS analysis on the geographic distances between European cities

```{r cmdscale2, eval=TRUE, warning=FALSE, message=FALSE}
loc <- cmdscale(eurodist) 
## Plots the MDS results in 2D plot. The minus is required in this example to 
## flip the plotting orientation.
plot(loc[,1], -loc[,2], type="n", xlab="", ylab="", main="cmdscale(eurodist)")
text(loc[,1], -loc[,2], rownames(loc), cex=0.8) 
```

### Principal Component Analysis (PCA)

Performs PCA analysis after scaling the data. It returns a list with class `prcomp` that contains five components: (1) the standard deviations (sdev) of the principal components, (2) the matrix of eigenvectors (rotation), (3) the principal component data (x), (4) the centering (center) and (5) scaling (scale) used.
```{r pca2, eval=TRUE, warning=FALSE, message=FALSE}
library(scatterplot3d)
pca <- prcomp(y, scale=TRUE)
names(pca)
summary(pca) # Prints variance summary for all principal components.
scatterplot3d(pca$x[,1:3], pch=20, color="blue") 
```

### Additional Exercises

See [here](http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual#TOC-Clustering-Exercises)

## Version Information

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

## References