How to Create a Heatmap in R

Introduction

A heatmap maps a numeric value to color at each row-column position. Decide whether colors should represent raw values or standardized values, and whether rows should be reordered by clustering. This lesson uses the same small matrix for base R and ggplot2, preserves a missing cell visibly, and treats scaling as a separate analytical choice.

Tested with R 4.5.1 and ggplot2 4.0.0. Run the code blocks in order. Install the named packages before starting; inputs are supplied below.

Keep labels outside the numeric matrix

Each row is a fictional site, each column a day, and values are counts. An unknown count is not zero.

m <- matrix(c(10, 12, 14, 20, NA, 24, 30, 32, 34), nrow = 3, byrow = TRUE,
            dimnames = list(c("Site A", "Site B", "Site C"), c("Day1", "Day2", "Day3")))
stopifnot(is.numeric(m), identical(dim(m), c(3L, 3L)), sum(is.na(m)) == 1L)
print(m)
#          Day1 Day2 Day3
# Site A     10   12   14
# Site B     20   NA   24
# Site C     30   32   34

Converting a data frame containing a text ID to a matrix can coerce every value to character. Select numeric measurement columns first. Do not replace missing values or unusual observations merely to obtain a more complete-looking image.

Draw base R heatmap values without hidden scaling

heatmap defaults can scale and reorder data. Here scale="none" and disabled dendrograms make the intended raw-value display explicit.

draw_base_heatmap <- function() {
  old <- par(oma = c(0, 0, 2, 0))
  on.exit(par(old))
  heatmap(m, Rowv = NA, Colv = NA, scale = "none", revC = TRUE,
          col = hcl.colors(24, "YlGnBu"), margins = c(5, 6), cexRow = 1, cexCol = 1,
          xlab = "Day", ylab = "", main = NULL)
  mtext("Raw counts (NA is blank)", side = 3, outer = TRUE, line = 0.3, font = 2, cex = 1.1)
}
draw_base_heatmap()
stopifnot(identical(unname(m[1, ]), c(10, 12, 14)))

The base plot preserves the displayed row/column sequence and leaves the NA cell blank. It has no numeric color legend by default; the next version adds one and labels cells so the value-to-color mapping is reviewable. The heatmap documentation explains scaling and dendrogram controls.

Add a legend and an explicit missing-cell label

expand.grid varies its first field fastest, matching R’s column-major matrix flattening. That keeps the values aligned with their site and day.

library(ggplot2)
long <- expand.grid(site = rownames(m), day = colnames(m), stringsAsFactors = FALSE)
long$value <- as.vector(m)
long$site <- factor(long$site, levels = rev(rownames(m)))
long$day <- factor(long$day, levels = colnames(m))
long$label <- ifelse(is.na(long$value), "NA", as.character(long$value))
p_heatmap <- ggplot(long, aes(day, site, fill = value)) +
  geom_tile(colour = "white", linewidth = 1) +
  geom_text(aes(label = label), colour = "black", size = 4) +
  scale_fill_gradient(low = "#edf8fb", high = "#66c2a4", na.value = "grey80", name = "Count") +
  coord_equal() + labs(x = "Day", y = "Site", title = "Raw counts and one missing cell") +
  theme_minimal(base_size = 12) + theme(plot.margin = margin(8, 18, 8, 8))
print(p_heatmap)
stopifnot(long$value[as.character(long$site) == "Site C" & long$day == "Day1"] == 30)
stopifnot(nrow(ggplot_build(p_heatmap)$data[[1]]) == 9L)

The legend represents raw counts, and NA is gray rather than a count of zero. Equal-sized tiles encode equal row-column positions, not geographic area.

Standardize only when relative patterns are the question

Row z-scores compare each value with its own row mean, in row standard-deviation units. They discard between-row magnitude information.

row_z <- function(x) {
  if (!is.matrix(x) || !is.numeric(x) || !nrow(x) || !ncol(x) ||
      any(!is.finite(x) & !is.na(x))) stop("Need a nonempty finite numeric matrix, with optional NA")
  spread <- apply(x, 1, sd, na.rm = TRUE)
  if (any(!is.finite(spread) | spread == 0)) stop("Each row needs variation and two observed values")
  sweep(sweep(x, 1, rowMeans(x, na.rm = TRUE), "-"), 1, spread, "/")
}
z <- row_z(m)
stopifnot(all(abs(rowMeans(z, na.rm = TRUE)) < 1e-12), identical(is.na(z), is.na(m)))
print(round(z, 3))  # A/C: -1,0,1; B: -0.707,NA,0.707
for (bad in list(matrix(1, 2, 2), matrix(NA_real_, 2, 2), matrix(numeric(0), 0, 2), matrix("x", 2, 2)))
  stopifnot(inherits(tryCatch(row_z(bad), error = identity), "error"))

Use a diverging scale centered on zero for standardized signed values, and label the units. Clustering introduces an additional distance/linkage choice and requires a justified missing-value policy. A heatmap alone does not establish clusters, causation or statistical significance.

Primary references

• R heatmap
• ggplot2 tiles
• R scale