---
title: "Introduction to heatmaply"
date: "`r Sys.Date()`"
author: "Tal Galili"
output:
  html_document:
    self_contained: yes
    toc: true    
---
<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Introduction to heatmaply}
-->


```{r, echo = FALSE, message = FALSE}
library(heatmaply)
library(knitr)
knitr::opts_chunk$set(
   # cache = TRUE,
   dpi = 60,
  comment = "#>",
  tidy = FALSE)

# http://stackoverflow.com/questions/24091735/why-pandoc-does-not-retrieve-the-image-file
# < ! -- rmarkdown v1 -->

```


**Author**: [Tal Galili](http://www.r-statistics.com/) ( Tal.Galili@gmail.com )



Introduction
============

A heatmap is a popular graphical method for visualizing high-dimensional data, in which a table of numbers are encoded as a grid of colored cells. The rows and columns of the matrix are ordered to highlight patterns and are often accompanied by dendrograms. Heatmaps are used in many fields for visualizing observations, correlations, missing values patterns, and more.

Interactive heatmaps allow the inspection of specific value by hovering the mouse over a cell, as well as zooming into a region of the heatmap by dragging a rectangle around the relevant area.

This work is based on the ggplot2 and plotly.js engine. It produces similar heatmaps as d3heatmap, with the advantage of speed (plotly.js is able to handle larger size matrix), and the ability to zoom from the dendrogram.


Installation
============

To install the stable version on CRAN:

```{r, eval = FALSE}
install.packages('heatmaply')
```

To install the GitHub version:

```{r, eval = FALSE}
# You'll need devtools
install.packages.2 <- function (pkg) if (!require(pkg)) install.packages(pkg);
install.packages.2('devtools')
# make sure you have Rtools installed first! if not, then run:
#install.packages('installr'); install.Rtools()

devtools::install_github("ropensci/plotly") 
devtools::install_github('talgalili/heatmaply')

```


And then you may load the package using:

```{r}
library("heatmaply")
```

Usage
======


Default
-------------

```{r}
library(heatmaply)
heatmaply(mtcars)
```

Because the labels are somewhat long, we need to manually fix the margins (hopefully this will be fixed in future versions of plotly). This also allows the inclusion of xlab/ylab/main texts.

```{r}
heatmaply(mtcars, xlab = "Features", ylab = "Cars", 
		main = "An example of title and xlab/ylab",
		margins = c(60,100,40,20))
```

We can use this with correlation. Notice the use of limits to set the range of the colors, and how we color the branches:

```{r}
heatmaply(cor(mtcars), margins = c(40, 40),
          k_col = 2, k_row = 2,
          limits = c(-1,1))
# Better to use:
# heatmaply_cor(cor(mtcars), margins = c(40, 40),
#           k_col = 2, k_row = 2)

```



Data transformation (scaling, normalize, and percentize)
===========

The variables in mtcars includes values that reflect different types of measurement, each with its own range (and meaning) of values. In such a case, it is best to transform the data so to have all the variables have comparable values.

scale
--------

If we would assume all variables come from some normal distribution, then scaling (i.e.: subtract the mean and divide by the standard deviation) would bring them all close to the standard normal distribution. In such a case, each value would reflect the distance from the mean in units of standard deviation. The "scale" parameter in heatmaply supports column and row scaling, and can be used as follows:

```{r}
heatmaply(mtcars, xlab = "Features", ylab = "Cars", 
          scale = "column",
		main = "Data transformation using 'scale'",
		margins = c(60,100,40,20))
```

normalize
--------



When variables in the data comes from possibly different (and non-normal) distributions, other transformations may be in order. For example, scaling on a binary variable with many 0's and just a few 1's will lead to a column with a very extreme value that would cause the color legend to be very skewed by that variable, leading to a masking of the distribution of the rest of the variables.

Another possibility is to use the `normalize` function to brings data to the 0 to 1 scale by subtracting the minimum and dividing by the maximum of all observations. This preserves the shape of each variable's distribution while making them easily comparable on the same "scale". Using the function on mtcars easily reveals columns with only two (am, vs) or three (gear, cyl) variables compared with variables that have a higher resolution of possible values:


```{r}
heatmaply(normalize(mtcars), xlab = "Features", ylab = "Cars", 
		main = "Data transformation using 'scale'",
		margins = c(60,100,40,20))
```

percentize
--------

An alternative to `normalize` is the `percentize` function. This is similar to ranking the variables, but instead of keeping the rank values, divide them by the maximal rank. This is done by using the `ecdf` of the variables on their own values, bringing each value to its  empirical percentile. The benefit of the percentize function is that each value has a relatively clear interpretation, it is the percent of observations that got that value or below it.

```{r}
heatmaply(percentize(mtcars), xlab = "Features", ylab = "Cars", 
		main = "Data transformation using 'scale'",
		margins = c(60,100,40,20))
```



is.na10 (missing values)
---------------

Reviewing missing values can easily be done using the `is.na10` function. When using it with heatmaply, it is often helpful to use grid_gap = 1.

```{r}

library(heatmaply)

heatmaply(is.na10(airquality), grid_gap = 1,
            k_col =3, k_row = 3,
            margins = c(55, 30),
            colors = c("grey80", "grey20"))
# Better to use:
# heatmaply_na(airquality, k_col =3, k_row = 3,
            # margins = c(55, 30))


# warning - using grid_color cannot handle a large matrix!
# airquality[1:10,] %>% is.na10 %>% 
#   heatmaply(color = c("white","black"), grid_color = "grey",
#             k_col =3, k_row = 3,
#             margins = c(40, 50)) 
# airquality %>% is.na10 %>% 
#   heatmaply(color = c("grey80", "grey20"), # grid_color = "grey",
#             k_col =3, k_row = 3,
#             margins = c(40, 50)) 
# 

```




Changing color palettes
-------------------

We can use colors other than the default `viridis`. For example, we may want to use other color palettes in order to get divergent colors for the correlations (these will, sadly, be less useful for colorblind people):

```{r}
# divergent_viridis_magma <- c(rev(viridis(100, begin = 0.3)), magma(100, begin = 0.3))
# rwb <- colorRampPalette(colors = c("darkred", "white", "darkgreen"))
# library(RColorBrewer)
# # display.brewer.pal(11, "BrBG")
# BrBG <- colorRampPalette(brewer.pal(11, "BrBG"))
# Spectral <- colorRampPalette(brewer.pal(11, "Spectral"))

heatmaply(cor(mtcars), margins = c(40, 40),
          k_col = 2, k_row = 2,
          colors = BrBG,
          limits = c(-1,1))
# Better to use:
# heatmaply_cor(cor(mtcars), margins = c(40, 40),
#           k_col = 2, k_row = 2)

```


Another example for using colors:


```{r}
heatmaply(percentize(mtcars), margins = c(40, 130),
          colors = heat.colors(100))
```

Or even more customized colors using `scale_fill_gradient_fun`:

```{r}
heatmaply(percentize(mtcars), margins = c(40, 130), 
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient2(low = "blue", high = "red", midpoint = 200, limits = c(0, 500)))

```



Customized dendrograms and side annotation
===================


Various seriation options
-------------------

heatmaply uses the `seriation` package to find an optimal ordering of rows and columns. Optimal means to optimize the Hamiltonian path length that is restricted by the dendrogram structure. This, in other words, means to rotate the branches so that the sum of distances between each adjacent leaf (label) will be minimized. This is related to a restricted version of the travelling salesman problem. 

The default options is "OLO" (Optimal leaf ordering) which optimizes the above criterion (in O(n^4)). Another option is "GW" ([Gruvaeus and Wainer](https://www.researchgate.net/publication/230266994_Two_additions_to_hierarchical_cluster_analysis)) which aims for the same goal but uses a potentially faster heuristic. The option "mean" gives the output we would get by default from heatmap functions in other packages such as `gplots::heatmap.2`. The option "none" gives us the dendrograms without any rotation that is based on the data matrix. 

```{r}
# The default of heatmaply:
heatmaply(percentize(mtcars)[1:10,], margins = c(40, 130),
          seriate = "OLO")
```

```{r}
# Similar to OLO but less optimal (since it is a heuristic)
heatmaply(percentize(mtcars)[1:10,], margin = c(40, 130),
          seriate = "GW")
```

```{r}
# the default by gplots::heatmaply.2
heatmaply(percentize(mtcars)[1:10,], margins = c(40, 130),
          seriate = "mean")
```

```{r}
# the default output from hclust
heatmaply(percentize(mtcars)[1:10,],  margins = c(40, 130),
          seriate = "none")
```

This works heavily relies on the seriation package (their [vignette](https://CRAN.R-project.org/package=seriation/vignettes/seriation.pdf) is well worth the read), and also lightly on the dendextend package (see [vignette](https://CRAN.R-project.org/package=dendextend/vignettes/introduction.html) )


Customized dendrograms using dendextend
-----------------------------

A user can supply their own dendrograms for the rows/columns of the heatmaply using the `Rowv` and the `Colv` parameters:

```{r}

x  <- as.matrix(datasets::mtcars)

# now let's spice up the dendrograms a bit:
library(dendextend)

row_dend  <- x %>% dist %>% hclust %>% as.dendrogram %>%
   set("branches_k_color", k = 3) %>% set("branches_lwd", c(1,3)) %>%
   ladderize
#    rotate_DendSer(ser_weight = dist(x))
col_dend  <- x %>% t %>% dist %>% hclust %>% as.dendrogram %>%
   set("branches_k_color", k = 2) %>% set("branches_lwd", c(1,2)) %>%
   ladderize
#    rotate_DendSer(ser_weight = dist(t(x)))

heatmaply(percentize(x), Rowv = row_dend, Colv = col_dend)


```



Replicating the dendrogram ordering of heatmap.2
-----------------------------------------------

The following example shows how to get the same result in heatmaply as with heatmap.2:

```{r}
x  <- as.matrix(datasets::mtcars)
gplots::heatmap.2(x, trace = "none", col = viridis(100), key = FALSE)
```

With heatmaply, the only difference is the side of the row dendrogram. This might be possible to modify in future versions of heatmaply:

```{r}
heatmaply::heatmaply(x, seriate = "mean")
```



Adding annotation based on additional factors using RowSideColors
-------------------------------

With heatmap.2

```{r}
# Example for using RowSideColors

x  <- as.matrix(datasets::mtcars)
rc <- colorspace::rainbow_hcl(nrow(x))

library(gplots)
library(viridis)
heatmap.2(x, trace = "none", col = viridis(100),
          RowSideColors=rc, key = FALSE)

```

With heatmaply

```{r}
heatmaply(x, seriate = "mean",
          RowSideColors=rc)

# heatmaply(x, seriate = "mean",
#           RowSideColors=factor(rc))


```

A more sophisticated heatmap (the hover at the top row doesn't work due to an issue with plotly. We hope this would get resolved in the future):

```{r}
heatmaply(x[,-c(8,9)], seriate = "mean",
          col_side_colors = c(rep(0,5), rep(1,4)),
          row_side_colors = x[,8:9])

```




Credit
========

This package is thanks to the amazing work done by MANY people in the open source community. Beyond the many people working on the pipeline of R, thanks should go to the plotly team, and especially to Carson Sievert and others working on the R package of plotly. Also, many of the design elements were inspired by the work done on heatmap, heatmap.2 and d3heatmap, so special thanks goes to the R core team, Gregory R. Warnes, and Joe Cheng from RStudio. The dendrogram side of the package is based on the work in dendextend, in which special thanks should go to Andrie de Vries for his original work on bringing dendrograms to ggplot2 (which evolved into the richer ggdend objects, as implemented in dendextend).


Contact
============


You are welcome to:

* submit suggestions and bug-reports at: <https://github.com/talgalili/heatmaply/issues>
* send a pull request on: <https://github.com/talgalili/heatmaply/>
* compose a friendly e-mail to: <tal.galili@gmail.com>


Latest news
===============

You can see the most recent changes to the package in the [NEWS.md file](https://github.com/talgalili/heatmaply/blob/master/NEWS.md)





Session info
=============

```{r, cache=FALSE}
sessionInfo()
```
