Last articles - STHDA : Hierarchical Clustering Essentials

Divisive Hierarchical Clustering Essentials

Wed, 06 Sep 2017 17:22:00 +0200

The divisive hierarchical clustering, also known as DIANA (DIvisive ANAlysis) is the inverse of agglomerative clustering (Chapter @ref(agglomerative-clustering)).

This article introduces the divisive clustering algorithms and provides practical examples showing how to compute divise clustering using R.

Algorithm

It starts by including all objects in a single large cluster. At each step of iteration, the most heterogeneous cluster is divided into two. The process is iterated until all objects are in their own cluster.

Recall that, divisive clustering is good at identifying large clusters while agglomerative clustering is good at identifying small clusters.

Computation

The R function diana() [cluster package] can be used to compute divisive clustering. It returns an object of class “diana” (see ?diana.object) which has also methods for the functions: print(), summary(), plot(), pltree(), as.dendrogram(), as.hclust() and cutree().

The output of DIANA can be visualized as dendrograms using the function fviz_dend() [factoextra package]. For example, the following R code shows how to computes and visualize divise clustering:

# Compute diana()
library(cluster)
res.diana <- diana(USArrests, stand = TRUE)
# Plot the dendrogram
library(factoextra)
fviz_dend(res.diana, cex = 0.5,
          k = 4, # Cut in four groups
          palette = "jco" # Color palette
          )

For interpreting dendrograms, read the “agglomerative clustering” chapter.

Heatmap - Static and Interactive: Absolute Guide

Wed, 06 Sep 2017 12:46:00 +0200

A heatmap (or heat map) is another way to visualize hierarchical clustering. It’s also called a false colored image, where data values are transformed to color scale.

Heat maps allow us to simultaneously visualize clusters of samples and features. First hierarchical clustering is done of both the rows and the columns of the data matrix. The columns/rows of the data matrix are re-ordered according to the hierarchical clustering result, putting similar observations close to each other. The blocks of ‘high’ and ‘low’ values are adjacent in the data matrix. Finally, a color scheme is applied for the visualization and the data matrix is displayed. Visualizing the data matrix in this way can help to find the variables that appear to be characteristic for each sample cluster.

Previously, we described how to visualize dendrograms. Here, we’ll demonstrate how to draw and arrange a heatmap in R.

Contents:

R Packages/functions for drawing heatmaps
Data preparation
R base heatmap: heatmap()
Enhanced heat maps: heatmap.2()
Pretty heat maps: pheatmap()
Interactive heat maps: d3heatmap()
Enhancing heatmaps using dendextend
Complex heatmap
Application to gene expression matrix
Visualizing the distribution of columns in matrix
Summary

R Packages/functions for drawing heatmaps

There are a multiple numbers of R packages and functions for drawing interactive and static heatmaps, including:

heatmap() [R base function, stats package]: Draws a simple heatmap
heatmap.2() [gplots R package]: Draws an enhanced heatmap compared to the R base function.
pheatmap() [pheatmap R package]: Draws pretty heatmaps and provides more control to change the appearance of heatmaps.
d3heatmap() [d3heatmap R package]: Draws an interactive/clickable heatmap
Heatmap() [ComplexHeatmap R/Bioconductor package]: Draws, annotates and arranges complex heatmaps (very useful for genomic data analysis)

Here, we start by describing the 5 R functions for drawing heatmaps. Next, we’ll focus on the ComplexHeatmap package, which provides a flexible solution to arrange and annotate multiple heatmaps. It allows also to visualize the association between different data from different sources.

Data preparation

We use mtcars data as a demo data set. We start by standardizing the data to make variables comparable:

df <- scale(mtcars)

R base heatmap: heatmap()

The built-in R heatmap() function [in stats package] can be used.

A simplified format is:

heatmap(x, scale = "row")

x: a numeric matrix
scale: a character indicating if the values should be centered and scaled in either the row direction or the column direction, or none. Allowed values are in c(“row”, “column”, “none”). Default is “row”.

# Default plot
heatmap(df, scale = "none")

In the plot above, high values are in red and low values are in yellow.

It’s possible to specify a color palette using the argument col, which can be defined as follow:

Using custom colors:

col<- colorRampPalette(c("red", "white", "blue"))(256)

Or, using RColorBrewer color palette:

library("RColorBrewer")
col <- colorRampPalette(brewer.pal(10, "RdYlBu"))(256)

Additionally, you can use the argument RowSideColors and ColSideColors to annotate rows and columns, respectively.

For example, in the the R code below will customize the heatmap as follow:

An RColorBrewer color palette name is used to change the appearance
The argument RowSideColors and ColSideColors are used to annotate rows and columns respectively. The expected values for these options are a vector containing color names specifying the classes for rows/columns.

# Use RColorBrewer color palette names
library("RColorBrewer")
col <- colorRampPalette(brewer.pal(10, "RdYlBu"))(256)
heatmap(df, scale = "none", col =  col, 
        RowSideColors = rep(c("blue", "pink"), each = 16),
        ColSideColors = c(rep("purple", 5), rep("orange", 6)))

Enhanced heat maps: heatmap.2()

The function heatmap.2() [in gplots package] provides many extensions to the standard R heatmap() function presented in the previous section.

# install.packages("gplots")
library("gplots")
heatmap.2(df, scale = "none", col = bluered(100), 
          trace = "none", density.info = "none")

Other arguments can be used including:

labRow, labCol
hclustfun: hclustfun=function(x) hclust(x, method=“ward”)

In the R code above, the bluered() function [in gplots package] is used to generate a smoothly varying set of colors. You can also use the following color generator functions:

colorpanel(n, low, mid, high)
- n: Desired number of color elements to be generated
- low, mid, high: Colors to use for the Lowest, middle, and highest values. mid may be omitted.
redgreen(n), greenred(n), bluered(n) and redblue(n)

Pretty heat maps: pheatmap()

First, install the pheatmap package: install.packages(“pheatmap”); then type this:

library("pheatmap")
pheatmap(df, cutree_rows = 4)

Arguments are available for changing the default clustering metric (“euclidean”) and method (“complete”). It’s also possible to annotate rows and columns using grouping variables.

Interactive heat maps: d3heatmap()

First, install the d3heatmap package: install.packages(“d3heatmap”); then type this:

library("d3heatmap")
d3heatmap(scale(mtcars), colors = "RdYlBu",
          k_row = 4, # Number of groups in rows
          k_col = 2 # Number of groups in columns
          )

The d3heamap() function makes it possible to:

Put the mouse on a heatmap cell of interest to view the row and the column names as well as the corresponding value.
Select an area for zooming. After zooming, click on the heatmap again to go back to the previous display

Enhancing heatmaps using dendextend

The package dendextend can be used to enhance functions from other packages. The mtcars data is used in the following sections. We’ll start by defining the order and the appearance for rows and columns using dendextend. These results are used in others functions from others packages.

The order and the appearance for rows and columns can be defined as follow:

library(dendextend)
# order for rows
Rowv  <- mtcars %>% scale %>% dist %>% hclust %>% as.dendrogram %>%
   set("branches_k_color", k = 3) %>% set("branches_lwd", 1.2) %>%
   ladderize
# Order for columns: We must transpose the data
Colv  <- mtcars %>% scale %>% t %>% dist %>% hclust %>% as.dendrogram %>%
   set("branches_k_color", k = 2, value = c("orange", "blue")) %>%
   set("branches_lwd", 1.2) %>%
   ladderize

The arguments above can be used in the functions below:

The standard heatmap() function [in stats package]:

heatmap(scale(mtcars), Rowv = Rowv, Colv = Colv,
        scale = "none")

The enhanced heatmap.2() function [in gplots package]:

library(gplots)
heatmap.2(scale(mtcars), scale = "none", col = bluered(100), 
          Rowv = Rowv, Colv = Colv,
          trace = "none", density.info = "none")

The interactive heatmap generator d3heatmap() function [in d3heatmap package]:

library("d3heatmap")
d3heatmap(scale(mtcars), colors = "RdBu",
          Rowv = Rowv, Colv = Colv)

Complex heatmap

ComplexHeatmap is an R/bioconductor package, developed by Zuguang Gu, which provides a flexible solution to arrange and annotate multiple heatmaps. It allows also to visualize the association between different data from different sources.

It can be installed as follow:

source("https://bioconductor.org/biocLite.R")
biocLite("ComplexHeatmap")

Simple heatmap

You can draw a simple heatmap as follow:

library(ComplexHeatmap)
Heatmap(df, 
        name = "mtcars", #title of legend
        column_title = "Variables", row_title = "Samples",
        row_names_gp = gpar(fontsize = 7) # Text size for row names
        )

Additional arguments:

show_row_names, show_column_names: whether to show row and column names, respectively. Default value is TRUE
show_row_hclust, show_column_hclust: logical value; whether to show row and column clusters. Default is TRUE
clustering_distance_rows, clustering_distance_columns: metric for clustering: “euclidean”, “maximum”, “manhattan”, “canberra”, “binary”, “minkowski”, “pearson”, “spearman”, “kendall”)
clustering_method_rows, clustering_method_columns: clustering methods: “ward.D”, “ward.D2”, “single”, “complete”, “average”, … (see ?hclust).

To specify a custom colors, you must use the the colorRamp2() function [circlize package], as follow:

library(circlize)
mycols <- colorRamp2(breaks = c(-2, 0, 2), 
                    colors = c("green", "white", "red"))
Heatmap(df, name = "mtcars", col = mycols)

It’s also possible to use RColorBrewer color palettes:

library("circlize")
library("RColorBrewer")
Heatmap(df, name = "mtcars",
        col = colorRamp2(c(-2, 0, 2), brewer.pal(n=3, name="RdBu")))

We can also customize the appearance of dendograms using the function color_branches() [dendextend package]:

library(dendextend)
row_dend = hclust(dist(df)) # row clustering
col_dend = hclust(dist(t(df))) # column clustering
Heatmap(df, name = "mtcars", 
        row_names_gp = gpar(fontsize = 6.5),
        cluster_rows = color_branches(row_dend, k = 4),
        cluster_columns = color_branches(col_dend, k = 2))

Splitting heatmap by rows

You can split the heatmap using either the k-means algorithm or a grouping variable.

It’s important to use the set.seed() function when performing k-means so that the results obtained can be reproduced precisely at a later time.

To split the dendrogram using k-means, type this:

# Divide into 2 groups
set.seed(2)
Heatmap(df, name = "mtcars", k = 2)

To split by a grouping variable, use the argument split. In the following example we’ll use the levels of the factor variable cyl [in mtcars data set] to split the heatmap by rows. Recall that the column cyl corresponds to the number of cylinders.

# split by a vector specifying rowgroups
Heatmap(df, name = "mtcars", split = mtcars$cyl,
        row_names_gp = gpar(fontsize = 7))

Note that, split can be also a data frame in which different combinations of levels split the rows of the heatmap.

# Split by combining multiple variables
Heatmap(df, name ="mtcars", 
        split = data.frame(cyl = mtcars$cyl, am = mtcars$am),
        row_names_gp = gpar(fontsize = 7))

It’s also possible to combine km and split:

Heatmap(df, name ="mtcars", col = mycol,
        km = 2, split =  mtcars$cyl)

If you want to use other partitioning method, rather than k-means, you can easily do it by just assigning the partitioning vector to split. In the R code below, we’ll use pam() function [cluster package]. pam() stands for Partitioning of the data into k clusters “around medoids”, a more robust version of K-means.

# install.packages("cluster")
library("cluster")
set.seed(2)
pa = pam(df, k = 3)
Heatmap(df, name = "mtcars", col = mycol,
        split = paste0("pam", pa$clustering))

Heatmap annotation

The HeatmapAnnotation class is used to define annotation on row or column. A simplified format is:

HeatmapAnnotation(df, name, col, show_legend)

df: a data.frame with column names
name: the name of the heatmap annotation
col: a list of colors which contains color mapping to columns in df

For the example below, we’ll transpose our data to have the observations in columns and the variables in rows.

df <- t(df)

Simple annotation

A vector, containing discrete or continuous values, is used to annotate rows or columns. We’ll use the qualitative variables cyl (levels = “4”, “5” and “8”) and am (levels = “0” and “1”), and the continuous variable mpg to annotate columns.

For each of these 3 variables, custom colors are defined as follow:

# Annotation data frame
annot_df <- data.frame(cyl = mtcars$cyl, am = mtcars$am, 
                       mpg = mtcars$mpg)
# Define colors for each levels of qualitative variables
# Define gradient color for continuous variable (mpg)
col = list(cyl = c("4" = "green", "6" = "gray", "8" = "darkred"),
            am = c("0" = "yellow", "1" = "orange"),
            mpg = circlize::colorRamp2(c(17, 25), 
                                       c("lightblue", "purple")) )
# Create the heatmap annotation
ha <- HeatmapAnnotation(annot_df, col = col)
# Combine the heatmap and the annotation
Heatmap(df, name = "mtcars",
        top_annotation = ha)

It’s possible to hide the annotation legend using the argument show_legend = FALSE as follow:

ha <- HeatmapAnnotation(annot_df, col = col, show_legend = FALSE)
Heatmap(df, name = "mtcars", top_annotation = ha)

Complex annotation

In this section we’ll see how to combine heatmap and some basic graphs to show the data distribution. For simple annotation graphics, the following functions can be used: anno_points(), anno_barplot(), anno_boxplot(), anno_density() and anno_histogram().

An example is shown below:

# Define some graphics to display the distribution of columns
.hist = anno_histogram(df, gp = gpar(fill = "lightblue"))
.density = anno_density(df, type = "line", gp = gpar(col = "blue"))
ha_mix_top = HeatmapAnnotation(hist = .hist, density = .density)
# Define some graphics to display the distribution of rows
.violin = anno_density(df, type = "violin", 
                       gp = gpar(fill = "lightblue"), which = "row")
.boxplot = anno_boxplot(df, which = "row")
ha_mix_right = HeatmapAnnotation(violin = .violin, bxplt = .boxplot,
                              which = "row", width = unit(4, "cm"))
# Combine annotation with heatmap
Heatmap(df, name = "mtcars", 
        column_names_gp = gpar(fontsize = 8),
        top_annotation = ha_mix_top, 
        top_annotation_height = unit(3.8, "cm")) + ha_mix_right

Combining multiple heatmaps

Multiple heatmaps can be arranged as follow:

# Heatmap 1
ht1 = Heatmap(df, name = "ht1", km = 2,
              column_names_gp = gpar(fontsize = 9))
# Heatmap 2
ht2 = Heatmap(df, name = "ht2", 
        col = circlize::colorRamp2(c(-2, 0, 2), c("green", "white", "red")),
        column_names_gp = gpar(fontsize = 9))
# Combine the two heatmaps
ht1 + ht2

You can use the option width = unit(3, “cm”)) to control the size of the heatmaps.

Note that when combining multiple heatmaps, the first heatmap is considered as the main heatmap. Some settings of the remaining heatmaps are auto-adjusted according to the setting of the main heatmap. These include: removing row clusters and titles, and adding splitting.

The draw() function can be used to customize the appearance of the final image:

draw(ht1 + ht2, 
    row_title = "Two heatmaps, row title", 
    row_title_gp = gpar(col = "red"),
    column_title = "Two heatmaps, column title", 
    column_title_side = "bottom",
    # Gap between heatmaps
    gap = unit(0.5, "cm"))

Legends can be removed using the arguments show_heatmap_legend = FALSE, show_annotation_legend = FALSE.

Application to gene expression matrix

In gene expression data, rows are genes and columns are samples. More information about genes can be attached after the expression heatmap such as gene length and type of genes.

expr <- readRDS(paste0(system.file(package = "ComplexHeatmap"),
                      "/extdata/gene_expression.rds"))
mat <- as.matrix(expr[, grep("cell", colnames(expr))])
type <- gsub("s\\d+_", "", colnames(mat))
ha = HeatmapAnnotation(df = data.frame(type = type))
Heatmap(mat, name = "expression", km = 5, top_annotation = ha, 
    top_annotation_height = unit(4, "mm"), 
    show_row_names = FALSE, show_column_names = FALSE) +
Heatmap(expr$length, name = "length", width = unit(5, "mm"),
    col = circlize::colorRamp2(c(0, 100000), c("white", "orange"))) +
Heatmap(expr$type, name = "type", width = unit(5, "mm")) +
Heatmap(expr$chr, name = "chr", width = unit(5, "mm"),
    col = circlize::rand_color(length(unique(expr$chr))))

It’s also possible to visualize genomic alterations and to integrate different molecular levels (gene expression, DNA methylation, …). Read the vignette, on Bioconductor, for further examples.

Visualizing the distribution of columns in matrix

densityHeatmap(scale(mtcars))

The dashed lines on the heatmap correspond to the five quantile numbers. The text for the five quantile levels are added in the right of the heatmap.

Summary

We described many functions for drawing heatmaps in R (from basic to complex heatmaps). A basic heatmap can be produced using either the R base function heatmap() or the function heatmap.2() [in the gplots package]. The pheatmap() function, in the package of the same name, creates pretty heatmaps, where ones has better control over some graphical parameters such as cell size.

The Heatmap() function [in ComplexHeatmap package] allows us to easily, draw, annotate and arrange complex heatmaps. This might be very useful in genomic fields.

Visualizing Dendrograms: Ultimate Guide

Wed, 06 Sep 2017 11:29:00 +0200

As described in previous chapters, a dendrogram is a tree-based representation of a data created using hierarchical clustering methods (Chapter @ref(agglomerative-clustering)). In this article, we provide R code for visualizing and customizing dendrograms. Additionally, we show how to save and to zoom a large dendrogram.

Contents:

Visualizing dendrograms
Case of dendrogram with large data sets
Manipulating dendrograms using dendextend
Summary

The Book:

Practical Guide to Cluster Analysis in R

We start by computing hierarchical clustering using the USArrests data sets:

# Load data
data(USArrests)
# Compute distances and hierarchical clustering
dd <- dist(scale(USArrests), method = "euclidean")
hc <- hclust(dd, method = "ward.D2")

To visualize the dendrogram, we’ll use the following R functions and packages:

fviz_dend()[in factoextra R package] to create easily a ggplot2-based beautiful dendrogram.
dendextend package to manipulate dendrograms

Before continuing, install the required package as follow:

install.packages(c("factoextra", "dendextend"))

Visualizing dendrograms

We’ll use the function fviz_dend()[in factoextra R package] to create easily a beautiful dendrogram using either the R base plot or ggplot2. It provides also an option for drawing circular dendrograms and phylogenic-like trees.

To create a basic dendrograms, type this:

library(factoextra)
fviz_dend(hc, cex = 0.5)

You can use the arguments main, sub, xlab, ylab to change plot titles as follow:

fviz_dend(hc, cex = 0.5, 
          main = "Dendrogram - ward.D2",
          xlab = "Objects", ylab = "Distance", sub = "")

To draw a horizontal dendrogram, type this:

fviz_dend(hc, cex = 0.5, horiz = TRUE)

It’s also possible to cut the tree at a given height for partitioning the data into multiple groups as described in the previous chapter: Hierarchical clustering (Chapter @ref(agglomerative-clustering)). In this case, it’s possible to color branches by groups and to add rectangle around each group.

For example:

fviz_dend(hc, k = 4, # Cut in four groups
          cex = 0.5, # label size
          k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE, # color labels by groups
          rect = TRUE, # Add rectangle around groups
          rect_border = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"), 
          rect_fill = TRUE)

To change the plot theme, use the argument ggtheme, which allowed values include ggplot2 official themes [ theme_gray(), theme_bw(), theme_minimal(), theme_classic(), theme_void()] or any other user-defined ggplot2 themes.

fviz_dend(hc, k = 4,                 # Cut in four groups
          cex = 0.5,                 # label size
          k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE,  # color labels by groups
          ggtheme = theme_gray()     # Change theme
          )

Allowed values for k_color include brewer palettes from RColorBrewer Package (e.g. “RdBu”, “Blues”, “Dark2”, “Set2”, …; ) and scientific journal palettes from ggsci R package (e.g.: “npg”, “aaas”, “lancet”, “jco”, “ucscgb”, “uchicago”, “simpsons” and “rickandmorty”).

In the R code below, we’ll change group colors using “jco” (journal of clinical oncology) color palette:

fviz_dend(hc, cex = 0.5, k = 4, # Cut in four groups
          k_colors = "jco")

If you want to draw a horizontal dendrogram with rectangle around clusters, use this:

fviz_dend(hc, k = 4, cex = 0.4, horiz = TRUE,  k_colors = "jco", 
          rect = TRUE, rect_border = "jco", rect_fill = TRUE)

Additionally, you can plot a circular dendrogram using the option type = “circular”.

fviz_dend(hc, cex = 0.5, k = 4, 
          k_colors = "jco", type = "circular")

To plot a phylogenic-like tree, use type = “phylogenic” and repel = TRUE (to avoid labels overplotting). This functionality requires the R package igraph. Make sure that it’s installed before typing the following R code.

require("igraph")
fviz_dend(hc, k = 4, k_colors = "jco",
          type = "phylogenic", repel = TRUE)

The default layout for phylogenic trees is “layout.auto”. Allowed values are one of: c(“layout.auto”, “layout_with_drl”, “layout_as_tree”, “layout.gem”, “layout.mds”, “layout_with_lgl”). To read more about these layouts, read the documentation of the igraph R package.

Let’s try phylo.layout = “layout.gem”:

require("igraph")
fviz_dend(hc, k = 4, # Cut in four groups
          k_colors = "jco",
          type = "phylogenic", repel = TRUE,
          phylo_layout = "layout.gem")

Case of dendrogram with large data sets

If you compute hierarchical clustering on a large data set, you might want to zoom in the dendrogram or to plot only a subset of the dendrogram.

Alternatively, you could also plot the dendrogram to a large page on a PDF, which can be zoomed without loss of resolution.

Zooming in the dendrogram

If you want to zoom in the first clusters, its possible to use the option xlim and ylim to limit the plot area. For example, type the code below:

fviz_dend(hc, xlim = c(1, 20), ylim = c(1, 8))

Plotting a sub-tree of dendrograms

To plot a sub-tree, we’ll follow the procedure below:

Create the whole dendrogram using fviz_dend() and save the result into an object, named dend_plot for example.
Use the R base function cut.dendrogram() to cut the dendrogram, at a given height (h), into multiple sub-trees. This returns a list with components $upper and $lower, the first is a truncated version of the original tree, also of class dendrogram, the latter a list with the branches obtained from cutting the tree, each a dendrogram.
Visualize sub-trees using fviz_dend().

The R code is as follow.

Cut the dendrogram and visualize the truncated version:

# Create a plot of the whole dendrogram,
# and extract the dendrogram data
dend_plot <- fviz_dend(hc, k = 4, # Cut in four groups
          cex = 0.5, # label size
          k_colors = "jco"
          )
dend_data <- attr(dend_plot, "dendrogram") # Extract dendrogram data
# Cut the dendrogram at height h = 10
dend_cuts <- cut(dend_data, h = 10)
# Visualize the truncated version containing
# two branches
fviz_dend(dend_cuts$upper)

Plot dendrograms sub-trees:

# Plot the whole dendrogram
print(dend_plot)

# Plot subtree 1
fviz_dend(dend_cuts$lower[[1]], main = "Subtree 1")
# Plot subtree 2
fviz_dend(dend_cuts$lower[[2]], main = "Subtree 2")

You can also plot circular trees as follow:

fviz_dend(dend_cuts$lower[[2]], type = "circular")

Saving dendrogram into a large PDF page

If you have a large dendrogram, you can save it to a large PDF page, which can be zoomed without loss of resolution.

pdf("dendrogram.pdf", width=30, height=15)            # Open a PDF
p <- fviz_dend(hc, k = 4, cex = 1, k_colors = "jco" ) # Do plotting
print(p)
dev.off()                                             # Close the PDF

Manipulating dendrograms using dendextend

The package dendextend provide functions for changing easily the appearance of a dendrogram and for comparing dendrograms.

In this section we’ll use the chaining operator (%>%) to simplify our code. The chaining operator turns x %>% f(y) into f(x, y) so you can use it to rewrite multiple operations such that they can be read from left-to-right, top-to-bottom. For instance, the results of the two R codes below are equivalent.

Standard R code for creating a dendrogram:

data <- scale(USArrests)
dist.res <- dist(data)
hc <- hclust(dist.res, method = "ward.D2")
dend <- as.dendrogram(hc)
plot(dend)

R code for creating a dendrogram using chaining operator:

library(dendextend)
dend <- USArrests[1:5,] %>% # data
        scale %>% # Scale the data
        dist %>% # calculate a distance matrix, 
        hclust(method = "ward.D2") %>% # Hierarchical clustering 
        as.dendrogram # Turn the object into a dendrogram.
plot(dend)

Functions to customize dendrograms: The function set() [in dendextend package] can be used to change the parameters of a dendrogram. The format is:

set(object, what, value)

object: a dendrogram object
what: a character indicating what is the property of the tree that should be set/updated
value: a vector with the value to set in the tree (the type of the value depends on the “what”).

Possible values for the argument what include:

Value for the argument what	Description
labels	set the labels
labels_colors and labels_cex	Set the color and the size of labels, respectively
leaves_pch, leaves_cex and leaves_col	set the point type, size and color for leaves, respectively
nodes_pch, nodes_cex and nodes_col	set the point type, size and color for nodes, respectively
hang_leaves	hang the leaves
branches_k_color	color the branches
branches_col, branches_lwd , branches_lty	Set the color, the line width and the line type of branches, respectively
by_labels_branches_col, by_labels_branches_lwd and by_labels_branches_lty	Set the color, the line width and the line type of branches with specific labels, respectively
clear_branches and clear_leaves	Clear branches and leaves, respectively

Examples:

library(dendextend)
# 1. Create a customized dendrogram
mycols <- c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07")
dend <-  as.dendrogram(hc) %>%
   set("branches_lwd", 1) %>% # Branches line width
   set("branches_k_color", mycols, k = 4) %>% # Color branches by groups
   set("labels_colors", mycols, k = 4) %>%  # Color labels by groups
   set("labels_cex", 0.5) # Change label size
# 2. Create plot 
fviz_dend(dend)

Summary

We described functions and packages for visualizing and customizing dendrograms including:

fviz_dend() [in factoextra R package], which provides convenient solutions for plotting easily a beautiful dendrogram. It can be used to create rectangular and circular dendrograms, as well as, a phylogenic tree.
and the dendextend package, which provides a flexible methods to customize dendrograms.

Additionally, we described how to plot a subset of large dendrograms.

Comparing Dendrograms: Essentials

Wed, 06 Sep 2017 03:11:00 +0200

After showing how to compute hierarchical clustering (Chapter @ref(agglomerative-clustering)), we describe, here, how to compare two dendrograms using the dendextend R package.

The dendextend package provides several functions for comparing dendrograms. Here, we’ll focus on two functions:

tanglegram() for visual comparison of two dendrograms
and cor.dendlist() for computing a correlation matrix between dendrograms.

Contents:

Data preparation
Dendrograms comparison

Data preparation

We’ll use the R base USArrests data sets and we start by standardizing the variables using the function scale() as follow:

df <- scale(USArrests)

To make readable the plots, generated in the next sections, we’ll work with a small random subset of the data set. Therefore, we’ll use the function sample() to randomly select 10 observations among the 50 observations contained in the data set:

# Subset containing 10 rows
set.seed(123)
ss <- sample(1:50, 10)
df <- df[ss,]

Dendrograms comparison

We start by creating a list of two dendrograms by computing hierarchical clustering (HC) using two different linkage methods (“average” and “ward.D2”). Next, we transform the results as dendrograms and create a list to hold the two dendrograms.

library(dendextend)

# Compute distance matrix
res.dist <- dist(df, method = "euclidean")

# Compute 2 hierarchical clusterings
hc1 <- hclust(res.dist, method = "average")
hc2 <- hclust(res.dist, method = "ward.D2")

# Create two dendrograms
dend1 <- as.dendrogram (hc1)
dend2 <- as.dendrogram (hc2)

# Create a list to hold dendrograms
dend_list <- dendlist(dend1, dend2)

Visual comparison of two dendrograms

To visually compare two dendrograms, we’ll use the following R functions [dendextend package]:

untangle(): finds the best layout to align dendrogram lists, using heuristic methods
tanglegram(): plots the two dendrograms, side by side, with their labels connected by lines.
entanglement(): computes the quality of the alignment of the two trees. Entanglement is a measure between 1 (full entanglement) and 0 (no entanglement). A lower entanglement coefficient corresponds to a good alignment.
Draw a tanglegram:

# Align and plot two dendrograms side by side
dendlist(dend1, dend2) %>%
  untangle(method = "step1side") %>% # Find the best alignment layout
  tanglegram()                       # Draw the two dendrograms

# Compute alignment quality. Lower value = good alignment quality
dendlist(dend1, dend2) %>%
  untangle(method = "step1side") %>% # Find the best alignment layout
  entanglement()                     # Alignment quality

## [1] 0.0384

Customized the tanglegram using many other options as follow:

dendlist(dend1, dend2) %>%
  untangle(method = "step1side") %>% 
  tanglegram(
    highlight_distinct_edges = FALSE, # Turn-off dashed lines
    common_subtrees_color_lines = FALSE, # Turn-off line colors
    common_subtrees_color_branches = TRUE # Color common branches 
    )

Note that, “unique” nodes, with a combination of labels/items not present in the other tree, are highlighted with dashed lines.

Note that, just because we can get two trees to have horizontal connecting lines, it doesn’t mean these trees are identical (or even very similar topologically).

In the following section, we’ll perform correlation analysis to measure the similarity between dendrograms.

Correlation matrix between a list of dendrogams

The function cor.dendlist() is used to compute “Baker” or “Cophenetic” correlation matrix between a list of trees. The value can range between -1 to 1. With near 0 values meaning that the two trees are not statistically similar.

# Cophenetic correlation matrix
cor.dendlist(dend_list, method = "cophenetic")

##       [,1]  [,2]
## [1,] 1.000 0.965
## [2,] 0.965 1.000

# Baker correlation matrix
cor.dendlist(dend_list, method = "baker")

##       [,1]  [,2]
## [1,] 1.000 0.962
## [2,] 0.962 1.000

The correlation between two trees can be also computed as follow:

# Cophenetic correlation coefficient
cor_cophenetic(dend1, dend2)

## [1] 0.965

# Baker correlation coefficient
cor_bakers_gamma(dend1, dend2)

## [1] 0.962

It’s also possible to compare simultaneously multiple dendrograms. A chaining operator %>% is used to run multiple function at the same time. It’s useful for simplifying the code:

# Create multiple dendrograms by chaining
dend1 <- df %>% dist %>% hclust("complete") %>% as.dendrogram
dend2 <- df %>% dist %>% hclust("single") %>% as.dendrogram
dend3 <- df %>% dist %>% hclust("average") %>% as.dendrogram
dend4 <- df %>% dist %>% hclust("centroid") %>% as.dendrogram
# Compute correlation matrix
dend_list <- dendlist("Complete" = dend1, "Single" = dend2,
                      "Average" = dend3, "Centroid" = dend4)
cors <- cor.dendlist(dend_list)
# Print correlation matrix
round(cors, 2)

##          Complete Single Average Centroid
## Complete     1.00   0.76    0.99     0.75
## Single       0.76   1.00    0.80     0.84
## Average      0.99   0.80    1.00     0.74
## Centroid     0.75   0.84    0.74     1.00

# Visualize the correlation matrix using corrplot package
library(corrplot)
corrplot(cors, "pie", "lower")

Agglomerative Clustering Essentials

Wed, 06 Sep 2017 02:44:00 +0200

The agglomerative clustering is the most common type of hierarchical clustering used to group objects in clusters based on their similarity. It’s also known as AGNES (Agglomerative Nesting). The algorithm starts by treating each object as a singleton cluster. Next, pairs of clusters are successively merged until all clusters have been merged into one big cluster containing all objects. The result is a tree-based representation of the objects, named dendrogram.

In this article we start by describing the agglomerative clustering algorithms. Next, we provide R lab sections with many examples for computing and visualizing hierarchical clustering. We continue by explaining how to interpret dendrogram. Finally, we provide R codes for cutting dendrograms into groups.

Contents:

Algorithm
Steps to agglomerative hierarchical clustering
Verify the cluster tree
Cut the dendrogram into different groups
Cluster R package
Application of hierarchical clustering to gene expression data analysis
Summary

Algorithm

Agglomerative clustering works in a “bottom-up” manner. That is, each object is initially considered as a single-element cluster (leaf). At each step of the algorithm, the two clusters that are the most similar are combined into a new bigger cluster (nodes). This procedure is iterated until all points are member of just one single big cluster (root) (see figure below).

The inverse of agglomerative clustering is divisive clustering, which is also known as DIANA (Divise Analysis) and it works in a “top-down” manner. It begins with the root, in which all objects are included in a single cluster. At each step of iteration, the most heterogeneous cluster is divided into two. The process is iterated until all objects are in their own cluster (see figure below).

Note that, agglomerative clustering is good at identifying small clusters. Divisive clustering is good at identifying large clusters. In this article, we’ll focus mainly on agglomerative hierarchical clustering.

Steps to agglomerative hierarchical clustering

We’ll follow the steps below to perform agglomerative hierarchical clustering using R software:

Preparing the data
Computing (dis)similarity information between every pair of objects in the data set.
Using linkage function to group objects into hierarchical cluster tree, based on the distance information generated at step 1. Objects/clusters that are in close proximity are linked together using the linkage function.
Determining where to cut the hierarchical tree into clusters. This creates a partition of the data.

We’ll describe each of these steps in the next section.

Data structure and preparation

The data should be a numeric matrix with:

rows representing observations (individuals);
and columns representing variables.

Here, we’ll use the R base USArrests data sets.

Note that, it’s generally recommended to standardize variables in the data set before performing subsequent analysis. Standardization makes variables comparable, when they are measured in different scales. For example one variable can measure the height in meter and another variable can measure the weight in kg. The R function scale() can be used for standardization, See ?scale documentation.

# Load the data
data("USArrests")
# Standardize the data
df <- scale(USArrests)
# Show the first 6 rows
head(df, nrow = 6)

##            Murder Assault UrbanPop     Rape
## Alabama    1.2426   0.783   -0.521 -0.00342
## Alaska     0.5079   1.107   -1.212  2.48420
## Arizona    0.0716   1.479    0.999  1.04288
## Arkansas   0.2323   0.231   -1.074 -0.18492
## California 0.2783   1.263    1.759  2.06782
## Colorado   0.0257   0.399    0.861  1.86497

Similarity measures

In order to decide which objects/clusters should be combined or divided, we need methods for measuring the similarity between objects.

There are many methods to calculate the (dis)similarity information, including Euclidean and manhattan distances (Chapter @ref(clustering-distance-measures)). In R software, you can use the function dist() to compute the distance between every pair of object in a data set. The results of this computation is known as a distance or dissimilarity matrix.

By default, the function dist() computes the Euclidean distance between objects; however, it’s possible to indicate other metrics using the argument method. See ?dist for more information.

For example, consider the R base data set USArrests, you can compute the distance matrix as follow:

# Compute the dissimilarity matrix
# df = the standardized data
res.dist <- dist(df, method = "euclidean")

Note that, the function () computes the distance between the rows of a data matrix using the specified distance measure method.

To see easily the distance information between objects, we reformat the results of the function dist() into a matrix using the as.matrix() function. In this matrix, value in the cell formed by the row i, the column j, represents the distance between object i and object j in the original data set. For instance, element 1,1 represents the distance between object 1 and itself (which is zero). Element 1,2 represents the distance between object 1 and object 2, and so on.

The R code below displays the first 6 rows and columns of the distance matrix:

as.matrix(res.dist)[1:6, 1:6]

##            Alabama Alaska Arizona Arkansas California Colorado
## Alabama       0.00   2.70    2.29     1.29       3.26     2.65
## Alaska        2.70   0.00    2.70     2.83       3.01     2.33
## Arizona       2.29   2.70    0.00     2.72       1.31     1.37
## Arkansas      1.29   2.83    2.72     0.00       3.76     2.83
## California    3.26   3.01    1.31     3.76       0.00     1.29
## Colorado      2.65   2.33    1.37     2.83       1.29     0.00

Linkage

The linkage function takes the distance information, returned by the function dist(), and groups pairs of objects into clusters based on their similarity. Next, these newly formed clusters are linked to each other to create bigger clusters. This process is iterated until all the objects in the original data set are linked together in a hierarchical tree.

For example, given a distance matrix “res.dist” generated by the function dist(), the R base function hclust() can be used to create the hierarchical tree.

hclust() can be used as follow:

res.hc <- hclust(d = res.dist, method = "ward.D2")

d: a dissimilarity structure as produced by the dist() function.
method: The agglomeration (linkage) method to be used for computing distance between clusters. Allowed values is one of “ward.D”, “ward.D2”, “single”, “complete”, “average”, “mcquitty”, “median” or “centroid”.

There are many cluster agglomeration methods (i.e, linkage methods). The most common linkage methods are described below.

Maximum or complete linkage: The distance between two clusters is defined as the maximum value of all pairwise distances between the elements in cluster 1 and the elements in cluster 2. It tends to produce more compact clusters.
Minimum or single linkage: The distance between two clusters is defined as the minimum value of all pairwise distances between the elements in cluster 1 and the elements in cluster 2. It tends to produce long, “loose” clusters.
Mean or average linkage: The distance between two clusters is defined as the average distance between the elements in cluster 1 and the elements in cluster 2.
Centroid linkage: The distance between two clusters is defined as the distance between the centroid for cluster 1 (a mean vector of length p variables) and the centroid for cluster 2.
Ward’s minimum variance method: It minimizes the total within-cluster variance. At each step the pair of clusters with minimum between-cluster distance are merged.

Note that, at each stage of the clustering process the two clusters, that have the smallest linkage distance, are linked together.

Complete linkage and Ward’s method are generally preferred.

Dendrogram

Dendrograms correspond to the graphical representation of the hierarchical tree generated by the function hclust(). Dendrogram can be produced in R using the base function plot(res.hc), where res.hc is the output of hclust(). Here, we’ll use the function fviz_dend()[ in factoextra R package] to produce a beautiful dendrogram.

First install factoextra by typing this: install.packages(“factoextra”); next visualize the dendrogram as follow:

# cex: label size
library("factoextra")
fviz_dend(res.hc, cex = 0.5)

In the dendrogram displayed above, each leaf corresponds to one object. As we move up the tree, objects that are similar to each other are combined into branches, which are themselves fused at a higher height.

The height of the fusion, provided on the vertical axis, indicates the (dis)similarity/distance between two objects/clusters. The higher the height of the fusion, the less similar the objects are. This height is known as the cophenetic distance between the two objects.

Note that, conclusions about the proximity of two objects can be drawn only based on the height where branches containing those two objects first are fused. We cannot use the proximity of two objects along the horizontal axis as a criteria of their similarity.

In order to identify sub-groups, we can cut the dendrogram at a certain height as described in the next sections.

Verify the cluster tree

After linking the objects in a data set into a hierarchical cluster tree, you might want to assess that the distances (i.e., heights) in the tree reflect the original distances accurately.

One way to measure how well the cluster tree generated by the hclust() function reflects your data is to compute the correlation between the cophenetic distances and the original distance data generated by the dist() function. If the clustering is valid, the linking of objects in the cluster tree should have a strong correlation with the distances between objects in the original distance matrix.

The closer the value of the correlation coefficient is to 1, the more accurately the clustering solution reflects your data. Values above 0.75 are felt to be good. The “average” linkage method appears to produce high values of this statistic. This may be one reason that it is so popular.

The R base function cophenetic() can be used to compute the cophenetic distances for hierarchical clustering.

# Compute cophentic distance
res.coph <- cophenetic(res.hc)
# Correlation between cophenetic distance and
# the original distance
cor(res.dist, res.coph)

## [1] 0.698

Execute the hclust() function again using the average linkage method. Next, call cophenetic() to evaluate the clustering solution.

res.hc2 <- hclust(res.dist, method = "average")
cor(res.dist, cophenetic(res.hc2))

## [1] 0.718

The correlation coefficient shows that using a different linkage method creates a tree that represents the original distances slightly better.

Cut the dendrogram into different groups

One of the problems with hierarchical clustering is that, it does not tell us how many clusters there are, or where to cut the dendrogram to form clusters.

You can cut the hierarchical tree at a given height in order to partition your data into clusters. The R base function cutree() can be used to cut a tree, generated by the hclust() function, into several groups either by specifying the desired number of groups or the cut height. It returns a vector containing the cluster number of each observation.

# Cut tree into 4 groups
grp <- cutree(res.hc, k = 4)
head(grp, n = 4)

##  Alabama   Alaska  Arizona Arkansas 
##        1        2        2        3

# Number of members in each cluster
table(grp)

## grp
##  1  2  3  4 
##  7 12 19 12

# Get the names for the members of cluster 1
rownames(df)[grp == 1]

## [1] "Alabama"        "Georgia"        "Louisiana"      "Mississippi"   
## [5] "North Carolina" "South Carolina" "Tennessee"

The result of the cuts can be visualized easily using the function fviz_dend() [in factoextra]:

# Cut in 4 groups and color by groups
fviz_dend(res.hc, k = 4, # Cut in four groups
          cex = 0.5, # label size
          k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE, # color labels by groups
          rect = TRUE # Add rectangle around groups
          )

Using the function fviz_cluster() [in factoextra], we can also visualize the result in a scatter plot. Observations are represented by points in the plot, using principal components. A frame is drawn around each cluster.

fviz_cluster(list(data = df, cluster = grp),
             palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"), 
             ellipse.type = "convex", # Concentration ellipse
             repel = TRUE, # Avoid label overplotting (slow)
             show.clust.cent = FALSE, ggtheme = theme_minimal())

Cluster R package

The R package cluster makes it easy to perform cluster analysis in R. It provides the function agnes() and diana() for computing agglomerative and divisive clustering, respectively. These functions perform all the necessary steps for you. You don’t need to execute the scale(), dist() and hclust() function separately.

The functions can be executed as follow:

library("cluster")
# Agglomerative Nesting (Hierarchical Clustering)
res.agnes <- agnes(x = USArrests, # data matrix
                   stand = TRUE, # Standardize the data
                   metric = "euclidean", # metric for distance matrix
                   method = "ward" # Linkage method
                   )
# DIvisive ANAlysis Clustering
res.diana <- diana(x = USArrests, # data matrix
                   stand = TRUE, # standardize the data
                   metric = "euclidean" # metric for distance matrix
                   )

After running agnes() and diana(), you can use the function fviz_dend()[in factoextra] to visualize the output:

fviz_dend(res.agnes, cex = 0.6, k = 4)

Application of hierarchical clustering to gene expression data analysis

In gene expression data analysis, clustering is generaly used as one of the first step to explore the data. We are interested in whether there are groups of genes or groups of samples that have similar gene expression patterns.

Several distance measures (Chapter @ref(clustering-distance-measures)) have been described for assessing the similarity or the dissimilarity between items, in order to decide which items have to be grouped together or not. These measures can be used to cluster genes or samples that are similar.

For most common clustering softwares, the default distance measure is the Euclidean distance. The most popular methods for gene expression data are to use log2(expression + 0.25), correlation distance and complete linkage clustering agglomerative-clustering.

Single and Complete linkage give the same dendrogram whether you use the raw data, the log of the data or any other transformation of the data that preserves the order because what matters is which ones have the smallest distance. The other methods are sensitive to the measurement scale.

Note that, when the data are scaled, the Euclidean distance of the z-scores is the same as correlation distance.

Pearson’s correlation is quite sensitive to outliers. When clustering genes, it is important to be aware of the possible impact of outliers. An alternative option is to use Spearman’s correlation instead of Pearson’s correlation.

In principle it is possible to cluster all the genes, although visualizing a huge dendrogram might be problematic. Usually, some type of preliminary analysis, such as differential expression analysis is used to select genes for clustering.

Selecting genes based on differential expression analysis removes genes which are likely to have only chance patterns. This should enhance the patterns found in the gene clusters.

Summary

Hierarchical clustering is a cluster analysis method, which produce a tree-based representation (i.e.: dendrogram) of a data. Objects in the dendrogram are linked together based on their similarity.

To perform hierarchical cluster analysis in R, the first step is to calculate the pairwise distance matrix using the function dist(). Next, the result of this computation is used by the hclust() function to produce the hierarchical tree. Finally, you can use the function fviz_dend() [in factoextra R package] to plot easily a beautiful dendrogram.

It’s also possible to cut the tree at a given height for partitioning the data into multiple groups (R function cutree()).