expression_heatmap

Introduction

矩阵热图和相关性热图形态相同，但回答的问题不同：相关性热图的矩阵是算出来的（两两相关系数），这里的矩阵是测出来的（表达量、丰度、强度）。

它的核心能力是双向聚类：行和列同时重排，把模式相似的挤到一起。所以图上出现的块状结构不是数据自带的顺序，是算法找出来的。

用 ComplexHeatmap 而不是 pheatmap 或 ggplot，省的是两件事：注释条，和横向拼接——同一批行可以把表达、甲基化、突变几个矩阵用 ht1 + ht2 + ht3 并排接起来，行顺序由第一个热图的聚类决定、其余自动跟着对齐。这是它名字里那个 "Complex" 的实际含义。

读的时候按这个顺序：

1. 先看色阶方向。 高低两端和中间锚点各是什么色。
2. 看块。 聚类后出现的连续色块就是"一起动"的行列组合。
3. 对照注释条。 这一步最要紧——Tumor 样本恰好聚成一块，说明矩阵里确实有和样本类型相关的信号；没聚在一起，那也是结论，而且往往更重要。
4. 看树状图的高度。 分支越早分开，两簇差别越大。
5. 看缺失格。 灰格是 NA，不是低值。

Example Data

矩阵热图要的输入是：一个数值矩阵（行 = 特征、列 = 样本），外加可选的行/列注释表。注释表靠行名和矩阵对齐，不靠顺序——但错位了也不会报错。

expression_heatmap 是一份模拟的长表：36 个基因 × 18 个样本，共 648 行，另有 6 个固定缺失值。基因分属 Cell cycle、Immune response 和 Lipid metabolism 三条模拟通路；Normal 与 Tumor 各 9 个样本。CSV 保留模拟 log-expression，热图在绘制前按基因做 z-score。

library(ComplexHeatmap)
library(circlize)
library(grid)
library(biopalette)

d <- read.csv(
  "https://assets.evanzhou.org/tessera/csv/expression_heatmap.csv",
  na.strings = "NA"
)

gene_levels <- unique(d$gene)
sample_levels <- unique(d$sample)

# 长表转矩阵。不能用 xtabs()：缺失格会被当成 0，而 0 和 NA 在表达数据里
# 不是一回事。
expr_raw <- matrix(
  NA_real_, nrow = length(gene_levels), ncol = length(sample_levels),
  dimnames = list(gene_levels, sample_levels)
)
expr_raw[cbind(match(d$gene, gene_levels), match(d$sample, sample_levels))] <-
  d$expression

# 每个基因各自标准化：颜色表示该基因在样本间相对偏高或偏低，不拿不同基因的
# 原始表达量直接比较。scale() 会保留那 6 个 NA。
expr_mat <- t(scale(t(expr_raw)))

gene_meta <- unique(d[c("gene", "pathway")])
sample_meta <- unique(d[c("sample", "type")])

row_anno_df <- data.frame(
  Pathway = gene_meta$pathway[match(rownames(expr_mat), gene_meta$gene)],
  row.names = rownames(expr_mat)
)
col_anno_df <- data.frame(
  Type = sample_meta$type[match(colnames(expr_mat), sample_meta$sample)],
  row.names = colnames(expr_mat)
)

Palettes

主体是逐基因 z-score，用发散色板 walter_white：负值、0、正值分别对应低于该基因均值、接近均值和高于均值。固定在 -2 / 0 / 2，不同热图才能使用同一种颜色语言。

注释条用定性色：Normal / Tumor 取 heat_light，三条 pathway 取 three_body。这个分工必须清楚：主体连续、注释定性，混了就分不清哪个通道在编码什么。

expr_palette <- get_palette("walter_white", type = "diverging")
expr_col_fun <- circlize::colorRamp2(
  breaks = c(-2, 0, 2),
  colors = expr_palette[c(1, 3, 5)]
)
expr_legend <- list(
  title = "Gene z-score",
  at = c(-2, 0, 2),
  labels = c("<= -2", "0", ">= 2"),
  direction = "vertical"
)

type_palette <- get_palette("heat_light", type = "qualitative")
pathway_palette <- get_palette("three_body", type = "qualitative")

# 注释色按名字绑定：聚类会重排行列，按位置给色会导致含义漂移
row_anno_colors <- list(
  Pathway = setNames(
    pathway_palette,
    c("Cell cycle", "Immune response", "Lipid metabolism")
  )
)
col_anno_colors <- list(
  Type = c(Normal = type_palette[2], Tumor = type_palette[1])
)

Recipe

1 · 基础热图

先看一眼矩阵本身有没有结构，没有注释干扰。

ComplexHeatmap 不是 ggplot 图层——Heatmap() 返回一个对象，要 draw() 才画出来。

#| fig: basic
#| fig-width: 6
#| fig-height: 6
basic_ht <- Heatmap(
  matrix = expr_mat,
  name   = "Gene z-score",
  col    = expr_col_fun,
  na_col = "grey90",              # NA 用中性灰，别用色阶里的任何一端
  cluster_rows    = TRUE,
  cluster_columns = TRUE,
  clustering_method_rows    = "complete",   # 方法要写进图注，换一个分块就变
  clustering_method_columns = "ward.D2",
  # 36 行标签在小图里仍然拥挤。关掉之后若要标重点基因，可用
  # 单独引出来：rowAnnotation(mark = anno_mark(at = which(...), labels = genes))
  # 会在侧边用连线把标签牵到对应的行上 —— 这是关行名之后的必要配套，不是装饰
  show_row_names    = FALSE,
  show_column_names = TRUE,
  column_names_gp   = gpar(fontsize = 8),
  heatmap_legend_param = expr_legend
)

draw(basic_ht)

2 · 加注释条与分块

注释条是判断"这个块有没有意义"的唯一依据。

#| fig: annotated
#| fig-width: 8
#| fig-height: 7
row_ha <- rowAnnotation(
  df = row_anno_df, col = row_anno_colors,
  annotation_name_gp = gpar(fontsize = 8)
)

col_ha <- HeatmapAnnotation(
  df = col_anno_df, col = col_anno_colors,
  annotation_name_gp = gpar(fontsize = 8)
)

annotated_ht <- Heatmap(
  matrix = expr_mat, name = "Gene z-score", col = expr_col_fun, na_col = "grey90",
  cluster_rows = TRUE, cluster_columns = TRUE,
  clustering_method_rows = "complete", clustering_method_columns = "ward.D2",
  show_row_names = FALSE, show_column_names = TRUE,
  column_names_gp = gpar(fontsize = 8),
  left_annotation = row_ha,
  top_annotation  = col_ha,
  # 固定尺寸而不是让它填满设备 —— 否则图例一多，热图本体就被挤扁
  heatmap_width  = unit(14, "cm"),
  heatmap_height = unit(14, "cm"),
  border = "white",
  heatmap_legend_param = expr_legend,
  use_raster = TRUE,      # 大矩阵栅格化，否则 PDF 里会有几万个矢量方块
  raster_quality = 2
)

draw(annotated_ht, heatmap_legend_side = "right", annotation_legend_side = "right")

3 · 按已知分组分块

不靠聚类，直接把 row_split / column_split 传进去，聚类只在每个块内部进行。

它比聚类更强硬：假定分组是对的，只看组内。适合"我已经知道分组、想看组间差异"，不适合"我还在找有没有分组"。

#| fig: split
#| fig-width: 8
#| fig-height: 7
split_ht <- Heatmap(
  matrix = expr_mat, name = "Gene z-score", col = expr_col_fun, na_col = "grey90",
  row_split    = row_anno_df$Pathway,   # 按通路切行
  column_split = col_anno_df$Type,      # 按样本类型切列
  cluster_rows = TRUE, cluster_columns = TRUE,
  show_row_names = FALSE, column_names_gp = gpar(fontsize = 8),
  left_annotation = row_ha, top_annotation = col_ha,
  heatmap_width = unit(14, "cm"), heatmap_height = unit(14, "cm"),
  border = "white", use_raster = TRUE, raster_quality = 2,
  heatmap_legend_param = expr_legend
)

draw(split_ht, heatmap_legend_side = "right", annotation_legend_side = "right")

导出。 走"开设备 → draw() → dev.off()"。Heatmap() 返回的是对象不是 ggplot，ggsave() 存不了；它也不像 base graphics 那样画完就已经在设备上。另有两件事只在导出时才暴露：画布一大图例容易被裁，draw() 时补 padding = unit(c(2, 20, 2, 20), "mm")；而 use_raster = TRUE 的主体是按设备尺寸栅格化的，导 300 dpi 时 raster_quality = 2 会糊，要往上调。

Constraints

• 块状结构不等于生物学。 这份 toy 刻意注入了三条 pathway 的组间信号，用来检验热图能否把它找回来；它不是生物学证据。真实数据仍要结合实验设计、批次和独立验证解释。
• 注释表错位不会报错。 行名对不上时图照样画得出来，只是注释贴到了别的基因或样本上，而图看起来完全正常。
• 聚类方法会改变结论。 complete 和 ward.D2 给出的分块可能完全不同。方法要写进图注，不然这张图不可复现。
• 没标准化就聚类，会被高表达基因主导。 表达量量级差几个数量级时尤其明显，多数情况下该先按行做 z-score。
• 标准化改变了问题。 逐基因 z-score 适合比较同一基因在样本间的相对变化，但不能再从颜色比较两个基因谁的原始表达量更高。这里把色阶固定在 -2–2，超出的值使用端点色。

三个 recipe 怎么比

1. 先看一眼矩阵有没有结构 → recipe 1。没有注释干扰，纯粹看聚类能不能分出块。
2. 有已知分组 → recipe 2。注释条是判断块有没有意义的唯一依据。
3. 已经知道分组、想看组间差异 → recipe 3。它假定分组是对的，只在组内聚类。