##########################
## Clustering Exercises ##
##########################

## Import a sample data set
## Download from GEO the Arabidopsis IAA treatment series "GSE1110" in TXT format. The direct link to the download is:
## ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/GSE1110/
## Uncompress the downloaded file.
## Import the data set into R
temp <- readLines("GSE1110_series_matrix.txt"); cat(temp[-grep("^!|^\"$", temp)], file="GSE1110clean.txt", sep="\n"); mydata <- read.delim("GSE1110clean.txt", header=T, sep="\t") # These import commands include a cleanup step to get rid of annotation lines and corrupted return signs.
rownames(mydata) <- mydata[,1]; mydata <- as.matrix(mydata[,-1]) # Assigns row indices and converts the data into a matrix object.

## Filtering
mydata <- mydata[apply(mydata>100, 1, sum)/length(mydata[1,])>0.5 & apply(log2(mydata), 1, IQR)>1.5,] # Retrieves all rows with high intensities (50% > 100) and high variability (IQR>1.5).

## Hierarchical clustering routine
y <- mydata
hr <- hclust(as.dist(1-cor(t(y), method="pearson")), method="complete"); hc <- hclust(as.dist(1-cor(y, method="spearman")), method="complete") # Generates row and column dendrograms.
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)] # Cuts the tree and creates color vector for clusters.
library(gplots); myheatcol <- redgreen(75) # Assign your favorite heatmap color scheme. Some useful examples: colorpanel(40, "darkblue", "yellow", "white"); heat.colors(75); cm.colors(75); rainbow(75); redgreen(75); library(RColorBrewer); rev(brewer.pal(9,"Blues")[-1]). Type demo.col(20) to see more color schemes.
heatmap.2(y, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=myheatcol, scale="row", density.info="none", trace="none", RowSideColors=mycolhc) # Creates heatmap for entire data set where the obtained clusters are indicated in the color bar.
x11(height=6, width=2); names(mycolhc) <- names(mycl); barplot(rep(10, max(mycl)), col=unique(mycolhc[hr$labels[hr$order]]), horiz=T, names=unique(mycl[hr$order])) # Prints color key for cluster assignments. The numbers next to the color boxes correspond to the cluster numbers in 'mycl'.
clid <- c(1,2); ysub <- y[names(mycl[mycl%in%clid]),]; hrsub <- hclust(as.dist(1-cor(t(ysub), method="pearson")), method="complete") # Select sub-cluster number (here: clid=c(1,2)) and generate corresponding dendrogram.
x11(); heatmap.2(ysub, Rowv=as.dendrogram(hrsub), Colv=as.dendrogram(hc), col=myheatcol, scale="row", density.info="none", trace="none", RowSideColors=mycolhc[mycl%in%clid]) # Create heatmap for chosen sub-cluster.
data.frame(Labels=rev(hrsub$labels[hrsub$order])) # Print out row labels in same order as shown in the heatmap.

## Hierarchical clustering
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/my.colorFct.R") # Import an alternative color scheme for the heatmap function.
mydatascale <- t(scale(t(mydata))) # Centers and scales data.
hr <- hclust(as.dist(1-cor(t(mydatascale), method="pearson")), method="complete") # Cluster rows by Pearson correlation.
hc <- hclust(as.dist(1-cor(mydatascale, method="spearman")), method="complete") # Clusters columns by Spearman correlation.
heatmap(mydata, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=my.colorFct(), scale="row") # Plot the data table as heatmap and the cluster results as dendrograms.
mycl <- cutree(hr, h=max(hr$height)/1.5) # Cut the tree at specific height and color the corresponding clusters in the heatmap color bar.
mycolhc <- sample(rainbow(256))
mycolhc <- mycolhc[as.vector(mycl)]
heatmap(mydata, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", RowSideColors=mycolhc) 

## Obtain significant clusters by pvclust bootstrap analysis
library(pvclust) # Loads the required pvclust package.
pv <- pvclust(scale(t(mydata)), method.dist="correlation", method.hclust="complete", nboot=10) # Perform the hierarchical cluster analysis. Due to time resrictions, we are using here only 10 bootstrap repetitions. Usually, one should use at least 1000 repetitions.
plot(pv, hang=-1); pvrect(pv, alpha=0.95) # Plots result as a dendrogram where the significant clusters are highlighted with red rectangles.
clsig <- unlist(pvpick(pv, alpha=0.95, pv="au", type="geq", max.only=TRUE)$clusters) # Retrieve members of significant clusters.
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/dendroCol.R") # Import tree coloring function.
dend_colored <- dendrapply(as.dendrogram(pv$hclust), dendroCol, keys=clsig, xPar="edgePar", bgr="black", fgr="red", pch=20) # Create dendrogram object where the significant clusters are labeled in red.
heatmap(mydata, Rowv=dend_colored, Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", RowSideColors=mycolhc) # Plot the heatmap from above, but with the significant clusters in red and the cluster bins from the tree cutting step in the color bar.
x11(height=12); heatmap.2(mydata, Rowv=dend_colored, Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", trace="none", RowSideColors=mycolhc) # Plot heatmap with heatmap.2() function which scales better for many entries.
mydatasort <- mydata[pv$hclust$labels[pv$hclust$order], hc$labels[hc$order]] # Sort rows in data table by 'dend_colored' and its colums by 'hc'.
x11(height=16, width=12); par(mfrow=c(1,2)); plot(dend_colored, horiz=T, yaxt="n"); image(scale(t(mydatasort)), col=my.colorFct(), xaxt="n",yaxt="n") # Plot heatmap with bootstrap tree in larger format using instead of heatmap the image function.
pdf("pvclust.pdf", height=21, width=6); plot(dend_colored, horiz=T, yaxt="n"); dev.off(); pdf("heatmap.pdf", height=20, width=6); image(scale(t(mydatasort)), col=my.colorFct(), xaxt="n",yaxt="n"); dev.off() # Save graphical results to two PDF files: 'pvclust.pdf' and'heatmap.pdf'.

## Compare PAM (K-means) with hierarchical clustering
library(cluster) # Loads required library.
mydist <- t(scale(t(mydata))) # Center and scale data.
mydist <- as.dist(1-cor(t(mydist), method="pearson")) # Generates distance matrix using Pearson correlation as distance method.
pamy <- pam(mydist, max(mycl)) # Clusters distance matrix into as many clusters as obtained by tree cutting step (6).
mycolkm <- sample(rainbow(256)); mycolkm <- mycolkm[as.vector(pamy$clustering)]; heatmap(mydata, Rowv=dend_colored, Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", RowSideColors=mycolkm) # Compare PAM clustering results with hierarchical clustering by labeling it in heatmap color bar.
pdf("pam.pdf", height=20, width=20); heatmap(mydata, Rowv=dend_colored, Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", RowSideColors=mycolkm); dev.off() # Save graphical results to PDF file: 'pvclust.pdf'.

## Compare SOM with hierarchical clustering
library(som) # Loads required library.
y <- t(scale(t(mydata))) # Center and scale data.
y.som <- som(y, xdim = 2, ydim = 3, topol = "hexa", neigh = "gaussian") # Performs SOM clustering.
plot(y.som) # Plots results.
pdf("som.pdf"); plot(y.som); dev.off() # Save plot to PDF: 'som.pdf'.
somclid <- as.numeric(paste(y.som$visual[,1], y.som$visual[,2], sep=""))+1 # Returns SOM cluster assignment in order of input data.
mycolsom <- sample(rainbow(256)); mycolsom <- mycolsom[somclid]; heatmap(mydata, Rowv=dend_colored, Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", RowSideColors=mycolsom) # Compare SAM clustering results with hierarchical clustering by labeling it in heatmap color bar.
pdf("somhc.pdf", height=20, width=20); heatmap(mydata, Rowv=dend_colored, Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", RowSideColors=mycolsom); dev.off() # Save graphical results to PDF file: 'somhc.pdf'.

## Compare PCA with SOM
pca <- prcomp(mydata, scale=T) # Performs principal component analysis after scaling the data.
summary(pca) # Prints variance summary for all principal components.
library(scatterplot3d) # Loads 3D library.
scatterplot3d(pca$x[,1:3], pch=20, color=mycolsom) # Plots PCA result in 3D. The SOM clusters are highlighted in their color.
pdf("pcasom.pdf"); scatterplot3d(pca$x[,1:3], pch=20, color=mycolsom); dev.off() # Saves PCA plot in PDF format: 'pcasom.pdf'.

## Compare MDS with HC, SOM and K-means
loc <- cmdscale(mydist, k = 3) # Performs MDS analysis and returns results for three dimensions.
x11(height=8, width=8, pointsize=12); par(mfrow=c(2,2)) # Sets plotting parameters.
plot(loc[,1:2], pch=20, col=mycolsom, main="MDS vs SOM 2D") # Plots MDS-SOM comparison in 2D. The SOM clusters are highlighted in their color.
scatterplot3d(loc, pch=20, color=mycolsom, main="MDS vs SOM 3D") # Plots MDS-SOM comparison in 3D.
scatterplot3d(loc, pch=20, color=mycolhc, main="MDS vs HC 3D") # Plots MDS-HC comparison.
scatterplot3d(loc, pch=20, color=mycolkm, main="MDS vs KM 3D") # Plots MDS-KM comparison.

## Fuzzy clustering
library(cluster) # Loads cluster library.
fannyy <- fanny(mydist, k= max(mycl), memb.exp = 1.5); round(fannyy$membership, 2); fannyy$clustering # Performs fuzzy clustering with as many coefficients as clusters were obtained by tree cutting step in HC. The hard clustering results are provided in the 'clustering' slot.
fannyyMA <- round(fannyy$membership, 2) > 0.3; apply(fannyyMA, 1, which) # Returns multiple cluster memberships for coefficient above a certain value (here >0.3).