2  Heatmap Scaling

2.1 Introduction

Heatmaps are ubiquitous in omics research for visualizing large-scale data matrices. However, one of the most common mistakes is using default scaling, which allows outliers to compress all meaningful variation into an indistinguishable range of colors. This chapter shows you how to properly handle scaling and outliers in heatmaps, and reveals a critical bug in many R heatmap functions.

Learning Objectives

By the end of this chapter, you will:

• Understand how outliers compress color scales in heatmaps
• Learn three robust scaling approaches: MAD scaling, quantile capping, and log transformation
• Discover the hidden dendrogram scaling bug in base R heatmap functions
• Know why scaling matters for clustering analysis
• Be able to use massageR::heat.clust() for proper integrated scaling
• Choose appropriate color scales for different data types

2.2 The Outlier Problem

Scenario: Metabolomics data (log2-transformed)

• Most metabolites: 6 to 10 (log2 scale)
• 1 outlier metabolite: ~10 (2^10 = 1000x higher!)

What happens with default scaling?

The extreme outlier compresses the color scale for all other values!

2.3 Visual Example: The Outlier Effect

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Without Outlier Handling

• Outliers dominate the color scale
• Most data compressed into narrow range
• Group differences invisible
• Patterns lost 😱

2.4 Solution 1: Robust MAD Scaling and Cutoffs

# Step 1: Robust scaling using MAD (Median Absolute Deviation)
expr_scaled <- expr_data %>%
  mutate(across(-Sample, scale_mad))

# Step 2: Identify which values will be capped
expr_mat_scaled <- expr_scaled %>% column_to_rownames("Sample") %>% as.matrix()
capped_cells <- (expr_mat_scaled < -3) | (expr_mat_scaled > 3)

# Step 3: Cap at ±3 (meaningful after MAD scaling!)
expr_capped <- expr_scaled %>% mutate(across(-Sample, ~ pmin(pmax(.x, -3), 3)))

# Step 4: Create symmetric breaks centered at 0
max_abs <- max(abs(range(expr_capped[,-1])))
breaks  <- seq(-max_abs, max_abs, length.out = 101)

# Create asterisk markers for capped values (in original order)
asterisk_matrix <- matrix("", nrow = nrow(capped_cells), ncol = ncol(capped_cells))
asterisk_matrix[capped_cells] <- "*"

# Create final plot with asterisk markers
# display_numbers uses the original data order, clustering is applied automatically
p <- pheatmap(expr_capped %>% column_to_rownames("Sample"),
              main = "MAD-scaled + capped at ±3 (* = capped)",
              color = colorRampPalette(rev(brewer.pal(11, "RdBu")))(100),
              breaks = breaks, scale = "none",
              display_numbers = asterisk_matrix,
              number_color = "black",
              fontsize_number = 14,
              silent = TRUE)

Robust scaling approach

• Use MAD instead of SD - not affected by outliers
• Center by median (robust)
• Scale by MAD (Median Absolute Deviation)
• Then cap at ±3 MAD (~99% of normal data)
• Outliers now exceed threshold!
• Set scale = "none" (already scaled!)

2.5 Solution 2: Range scaling and Quantile-Based cut-off

# Step 1: Identify values outside 5-95 percentiles PER COLUMN
expr_mat_raw <- expr_data %>% column_to_rownames("Sample") %>% as.matrix()
capped_cells_q <- apply(expr_mat_raw, 2, function(x) {
  q_lower <- quantile(x, 0.05, na.rm = TRUE)
  q_upper <- quantile(x, 0.95, na.rm = TRUE)
  (x < q_lower) | (x > q_upper)
})

# Step 2: Cap at 5th and 95th percentiles PER COLUMN (metabolite)
expr_capped <- expr_data %>%
  mutate(across(-Sample, ~ cap_quantiles(.x, lower = 0.05, upper = 0.95)))

# Step 3: Range scaling (min-max normalization to [0,1])
expr_quantile <- expr_capped %>% mutate(across(-Sample, ~ (.x - min(.x)) / (max(.x) - min(.x))))

# Create asterisk markers for capped values (in original order)
asterisk_matrix_q <- matrix("", nrow = nrow(capped_cells_q), ncol = ncol(capped_cells_q))
asterisk_matrix_q[capped_cells_q] <- "*"

# Create final plot with asterisk markers
p <- pheatmap(expr_quantile %>% column_to_rownames("Sample"),
              main = "Capped at 5-95 percentiles + range-scaled (* = capped)",
              color = rev(viridis::magma(100)),
              display_numbers = asterisk_matrix_q,
              number_color = "white",
              fontsize_number = 14,
              silent = TRUE)

Quantile capping + range scaling

• Cap first to remove outliers per metabolite
• Then range scale to use full [0,1] color scale
• More robust to outliers than variance scaling
• Good for non-normal data
• Common quantiles: 5-95% or 2-98%

2.6 Solution 3: Log Transformation

For positive values only (e.g., counts, intensities)

# Log transform BEFORE plotting (metabolomics data is positive-only)
expr_log_sol4 <- expr_data %>%
  mutate(across(-Sample, ~ log2(.x + 1)))  # +1 to handle zeros

p <- pheatmap(expr_log_sol4 %>% column_to_rownames("Sample"),
              main = "Log2 transformed metabolite intensities",
              color = rev(viridis::magma(100)),
              silent = TRUE)

For count/intensity data

• Compresses wide ranges
• Add +1 to handle zeros
• Common for RNA-seq, proteomics
• Use log2, log10, or ln

2.7 The Dendrogram Scaling Trap

Hidden Technical Issue

Critical R bug: Functions like heatmap(), heatmap.2(), and heatplot() have a dangerous inconsistency:

• The scale parameter affects color visualization
• But NOT dendrogram calculation!

Result: Dendrograms cluster on unscaled data while colors show scaled data!

P.S: pheatmap() seems to apply scaling and cropping before clustering!

2.8 Why Scaling Matters for Clustering

The Problem: High-variance features dominate correlations

Without scaling, features with large values dominate sample correlations!

corrplot 0.95 loaded

2.9 Why Scaling Matters for Clustering (2)

Key point: Without scaling, high-variance metabolites completely dominate the correlation calculation between samples!

Scaling ensures all metabolites contribute equally to sample clustering.

2.10 Using massageR::heat.clust

Better Approach: heat.clust

The massageR package provides heat.clust() which handles scaling and dendrogram calculation correctly in one step!

Key advantages:

• Scales data and calculates dendrograms together
• Controls exactly where limits are applied (data and/or dendrograms)
• Returns pre-computed dendrograms
• Works seamlessly with pheatmap

2.11 heat.clust with pheatmap

library(massageR)

# Convert tibble to matrix for heat.clust
expr_matrix <- expr_data %>% column_to_rownames("Sample") %>% as.matrix()

# Use heat.clust with robust MAD scaling
z <- heat.clust(expr_matrix,
                scaledim = "column",           # Scale by column
                zlim = c(-3, 3),               # Cap at ±3 MAD
                zlim_select = c("dend", "outdata"),  # Apply to both
                reorder = c(),                 # Reorder dendrograms off for consistency
                distfun = function(x) dist(x),
                hclustfun = function(x) hclust(x, method = "complete"),
                scalefun = scale_mad)          # Use MAD scaling instead of default

max_abs <- max(abs(range(z$data)))
breaks  <- seq(-max_abs, max_abs, length.out = 101)

# Use with pheatmap
p <- pheatmap(z$data,
              cluster_rows = as.hclust(z$Rowv),
              cluster_cols = as.hclust(z$Colv),
              scale = "none",
              color = colorRampPalette(rev(brewer.pal(11, "RdBu")))(100),
              breaks = breaks,
              main = "heat.clust + pheatmap: Properly scaled!",
              silent = TRUE)

One-step workflow

• Scales data automatically
• Calculates dendrograms on scaled data
• Caps at specified zlim
• Returns everything needed for pheatmap
• Ensures consistency throughout

2.12 Comparison: Before and After

Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine

2.13 Best Practices & Recommendations

Workflow

1. Inspect data distribution before making heatmap
2. Consider log transformation
3. Scale data BEFORE passing to heatmap function
   • Use MAD (robust) if data has outliers: (x - median(x)) / mad(x)
   • Use SD if data is clean: scale(x)
4. Cap extremes
   • Cap at ±3 MAD for robust scaling with outliers
   • Cap using quantiles (5-95%) if using range scaling
5. Calculate dendrograms on the same scaled and capped data
6. Consider massageR::heat.clust() for automatic proper scaling workflow
7. Use appropriate palette:
   • Often centering data highlights contrasts
   • Diverging color scale for data that has been centered (red-white-blue)
   • Sequential color scale for one-directional data (viridis, magma)

When NOT to Cap

• If outliers are biologically meaningful (rare events)
• Small datasets where each value matters
• When you want to highlight extreme values

2.14 Summary

Key Takeaways

1. Outliers compress color scales - default scaling makes all non-outlier data invisible
2. Use robust MAD scaling - (x - median(x)) / mad(x) is not affected by outliers
3. Cap at ±3 MAD - retains ~99% of normal data while handling extremes
4. Alternative: quantile capping - cap at 5-95 percentiles before range scaling
5. Log transform positive data - metabolomics, RNA-seq, proteomics intensities
6. Beware the dendrogram bug - base R heatmap functions don’t scale before clustering!
7. Use massageR::heat.clust() - handles scaling and dendrograms correctly in one step
8. Match color scales to data:
   • Diverging (RdBu) for centered data
   • Sequential (magma/viridis) for positive-only data

2.15 Exercises

Try It Yourself

1. Create a heatmap with simulated data containing outliers (use rnorm() plus a few extreme values)
2. Compare three approaches:
   • Default scaling with pheatmap(..., scale = "row")
   • MAD scaling with scale_mad() function and capping at ±3
   • Quantile capping at 5-95 percentiles
3. Check if your dendrograms change between methods
4. Try massageR::heat.clust() and verify it produces consistent results

2.16 Further Reading

• massageR package documentation - for heat.clust() function
• Why you should use MAD over SD - robust statistics
• ComplexHeatmap book - advanced heatmap customization
• Wilkinson & Friendly (2009). “The History of the Cluster Heat Map”. The American Statistician - origins and best practices