---
title: "Clustering methods"
author: "Pengyi Yang"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

## Learning objective
The aim in this lecture is to understand and apply various clustering algorithms for exploratory data analysis.

## Libraries to load
```{r, warning=FALSE, message=FALSE}
library(limma)
library(gplots)
library(e1071)
library(shiny)
```


## Application of clustering algorithms to Leukemia dataset
ExpressionSet containing gene expresion data from 60 bone marrow samples of patients with one
of the four main types of leukemia (ALL, AML, CLL, CML) or no-leukemia controls.

Platform: Affymetrix Human Genome U133 Plus 2.0  

Annotation: genemapperhgu133plus2 (CDF from GATExplorer)

Mapping: Gene Ensembl ID (20172 features)

Tissue: Bone Marrow

Cell type: Mononuclear cells isolated by Ficoll density centrifugation

Disease type:  
1. Acute Lymphoblastic Leukemia (ALL). Subtype: c-ALL / pre-B-ALL without t(9;22)  
2. Acute Myeloid Leukemia (AML). Subtype: Normal karyotype  
3. Chronic Lymphocytic Leukemia (CLL)  
4. Chronic Myeloid Leukemia (CML)  
5. Non-leukemia and healthy bone marrow (NoL)  

All samples were obtained from untreated patients at the time of diagnosis.
Preprocessing: The microarrays were normalized with RMA using a redefined probe mapping from
Affymetrix probesets to Ensembl genes (Ensembl IDs ENSG). This alternative Chip Definition
File (CDF) with complete unambiguous mapping of microarray probes to genes (GeneMapper)
is available at GATExplorer (http://bioinfow.dep.usal.es/xgate/mapping/mapping.php) (Risueno et
al. 2010).


```{r, warning=FALSE, message=FALSE}
# load data
library(leukemiasEset)
data(leukemiasEset)
# extract expression matrix from the "ExpressionSet object"
leukemia.mat <- exprs(leukemiasEset)
# obtain and summarize phenodata
phenodata <- pData(leukemiasEset)
```


## hierarchical clustering analysis of Leukemia dataset
```{r, warning=FALSE, message=FALSE, fig.width=12, fig.height=12}
labs <- phenodata[,"LeukemiaType"]
#
par(mfrow=c(3,1))
data.dist=dist(t(leukemia.mat))
plot(hclust(data.dist), labels=labs, main="Complete Linkage", xlab="", sub="",ylab="", cex=1.3)
plot(hclust(data.dist, method="average"), labels=labs, main="Average Linkage", xlab="", sub="",ylab="", cex=1.3)
plot(hclust(data.dist, method="single"), labels=labs,  main="Single Linkage", xlab="", sub="",ylab="", cex=1.3)

hc.out=hclust(dist(t(leukemia.mat)))
hc.clusters = cutree(hc.out, 5)
table(hc.clusters, labs)

par(mfrow=c(1,1))
plot(hc.out, labels = labs)
abline(h=139, col="red")
hc.out
```

### Identify and select differentially expressed genes from Leukemia dataset using LIMMA package

```{r, warning=FALSE, message=FALSE}
# create the design matrix
f = phenodata[colnames(leukemia.mat), 3]
design = model.matrix(~ f - 1)
colnames(design) <- c("ALL", "AML", "CLL", "CML", "NoL")
fit <- lmFit(leukemia.mat , design) # fitting data to a linear model
contrast.matrix <- makeContrasts(ALL-NoL, AML-NoL, CLL-NoL, CML-NoL, levels=design)  # defining group comparisons
fit2 <- contrasts.fit(fit, contrast.matrix) # convert the coefficients of the design matrix into these 4 contrasts which are to be tested equal to zero
fit2 <- eBayes(fit2) # With a series of related parameter estimates & standard errors, eBayes computes moderated t-statistics of differential expression by empirical Bayes shrinkage of the standard errors towards a common value.
leukemia.topDE <-leukemia.mat[order(fit2$F.p.value), ][1:500,] # Select the top 500 most differentially expressed genes
dim(leukemia.topDE)
```

### Use heatmap to visualize hierarchical clustering results
```{r, fig.width=12, fig.height=12}
patientcolors = sapply(as.character(f), switch, ALL = "red", AML =  "orange", CLL = "blue",  CML="green", NoL="gray")
heatmap.2(leukemia.topDE, col=bluered(75), ColSideColors=patientcolors, density.info="none", trace="none", na.color = "black", margins=c(8, 8), main="Clustering of top500 DE genes", dendrogram = "column")
```


# output session information
```{r, echo=FALSE}
sessionInfo()
```





