---
title: "k-means and c-means clustering and their use in network reconstruction"
author: "Pengyi Yang"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

# Topic 1. k-means based clustering algorithms

### Apply k-means clustering on iris data
```{r, fig.width=10, fig.height=10}
dim(iris)
head(iris)
table(iris$Species)
class(iris)

# set up plots
par(mfrow=c(2,2))
# plot original data
speciesColors = sapply(as.character(iris$Species), switch, setosa = "red", versicolor = "blue",  virginica="green")
data.mat <- cbind(iris$Petal.Length, iris$Sepal.Width, iris$Species)
colnames(data.mat) <- c("Petal.Length", "Sepal.Width", "Species")
plot(data.mat[,-3], col=speciesColors, main="orignal data with class information")

# apply k-means with k=2
set.seed(1)
km.out2 <- kmeans(data.mat[,-3], centers=2)
plot(data.mat[,-3], col=(km.out2$cluster+1), main="k-means clustering results with k=2", xlab="f1", ylab="f2", pch=20, cex=2)

# apply k-means with k=3
set.seed(1)
km.out3 <- kmeans(data.mat[,-3], centers=3)
plot(data.mat[,-3], col=(km.out3$cluster+1), main="k-means clustering results with k=3", xlab="f1", ylab="f2", pch=20, cex=2)

# apply k-means with k=4
set.seed(1)
km.out4 <- kmeans(data.mat[,-3], centers=4)
plot(data.mat[,-3], col=(km.out4$cluster+1), main="k-means clustering results with k=4", xlab="f1", ylab="f2", pch=20, cex=2)
```

### Cluster statistics
```{r}
# what is the output from k-means clustering?
km.out2

# check between and within cluster sum of squares.
km.out2$betweenss
km.out2$tot.withinss

km.out3$betweenss
km.out3$tot.withinss

km.out4$betweenss
km.out4$tot.withinss
```


# Topic 2. Application of c-means clustering.

notice that we use the **membership** score to control the size of the point to plot. k-means clustering algorithm will not create membership score and will create a hard clustering whereas c-means gives membership score which can be used as a confidence meansure of clustering for each sample. 

```{r fig.width=10, fig.height=10}
library(e1071)
# set up plots
par(mfrow=c(2,2))
# apply c-means with k=1
set.seed(1)
cm.out2 <- cmeans(data.mat[,-3], centers=2, m = 2)
plot(data.mat[,-3], col=(cm.out2$cluster+1), main="c-means clustering results with k=2", xlab="f1", ylab="f2", pch=20, cex=apply(cm.out2$membership+0.5, 1, max)^2)

# apply c-means with k=3
set.seed(1)
cm.out3 <- cmeans(data.mat[,-3], centers=3)
plot(data.mat[,-3], col=(cm.out3$cluster+1), main="c-means clustering results with k=3", xlab="f1", ylab="f2", pch=20, cex=apply(cm.out2$membership+0.5, 1, max)^2)

# apply c-means with k=4
set.seed(1)
cm.out4 <- cmeans(data.mat[,-3], centers=4)
plot(data.mat[,-3], col=(cm.out4$cluster+1), main="c-means clustering results with k=4", xlab="f1", ylab="f2", pch=20, cex=apply(cm.out2$membership+0.5, 1, max)^2)
```


## Clustering results validation
### data structure based metrics
The internal measures include the connectivity, and Silhouette Width, and Dunn Index. The connectivity indicates the degree of connectedness of the clusters, as determined by the k-nearest neighbors. The neighbSize argument specifies the number of
neighbors to use. The connectivity has a value between 0 and infinity and should be minimized.
Both the Silhouette Width and the Dunn Index combine measures of compactness and
separation of the clusters. The Silhouette Width is the average of each observation's Silhouette value. The Silhouette value measures the degree of confidence in a particular clustering assignment and lies in the interval [-1,1], with well-clustered observations having values near 1 and poorly clustered observations having values near -1. 

```{r, warning=FALSE, fig.width=10, fig.height=10}
set.seed(1)
library(clValid)
intern <- clValid(adipocyte.scaled, nClust=5:20, validation="internal", clMethods="kmeans", maxitems = 2000)
summary(intern)
optimalScores(intern)
# set up plots
par(mfrow=c(2,2))
plot(intern)


library(clValid) # clValid is a package that containing various cluster validation methods
intern <- clValid(data.mat[,-3], nClust=2:10, validation="internal", clMethods=c("hierarchical","kmeans", "pam"))
summary(intern)
optimalScores(intern)
# set up plots
par(mfrow=c(2,2))
plot(intern)
```


### stability based metrics
The stability measures evaluate the stability of a clustering result by comparing it with the clusters obtained by removing one column at a time. These measures include the average proportion of non-overlap (APN), the average distance (AD), the average distance between means (ADM), and the figure of merit (FOM). The APN, AD, and ADM are all based on the cross-classification table of the original clustering with the clustering based on the removal of one column. The APN measures the average proportion of observations not placed in the same cluster under both cases, while the AD measures the average distance between observations placed in the same cluster under both cases and the ADM measures the average distance between cluster centers for observations placed in the same cluster under both cases. The FOM measures the average intra-cluster variance of the deleted column, where the clustering is based on the remaining (undeleted) columns. In all cases the average is taken over all the deleted columns, and all measures should be minimized.

```{r, warning=FALSE, fig.width=10, fig.height=10}
par(mfrow=c(2,2))
stab <- clValid(iris[,-5], nClust=2:10, validation="stability", clMethods=c("hierarchical", "kmeans", "pam"))
optimalScores(stab)
plot(stab)
```





# Topic 3. Use k-means to clustering data and infer networks

```{r}
# (1)
library(ClueR)
data(adipocyte)

# (2)
adipocyte.exp <- adipocyte[,-ncol(adipocyte)]
idx <- which(rowSums(abs(adipocyte.exp) > 1) > 0)
adipocyte.reg <- adipocyte.exp[idx,]
# number of genes before filtering
nrow(adipocyte.exp)
# number of genes after filtering
nrow(adipocyte.reg)
# number of time points in the data
ncol(adipocyte.reg)

# (3)
standardize <- function(mat) {
  means <- apply(mat, 1, mean)
  stds <- apply(mat, 1, sd)
  tmp <- sweep(mat, 1, means, FUN="-")
  mat.stand <- sweep(tmp, 1, stds, FUN="/")
  return(mat.stand)
}

adipocyte.scaled <- standardize(adipocyte.reg)

# (4) 
set.seed(1)
km <- kmeans(adipocyte.scaled, centers = 13)

# (5)
km.gene.clustered <- names(which(km$cluster == km$cluster["ISL1"]))

# (6)
library(igraph)
mat <- cor(t(adipocyte.scaled[km.gene.clustered,]))
mat[mat<0.99]=0

# Make an Igraph object from this matrix:
network <- graph_from_adjacency_matrix(mat, weighted=T, mode="undirected", diag=F)
V(network)$label.cex <- 0.5
l <- layout_on_sphere(network)
plot(network, layout=l, cex=1, vertex.size=3)
```

