R语言ComplexHeatmap实战：样品重复性差时的热图分组聚类技巧

生物信息学研究中，热图（Heatmap）是最常用的数据可视化工具之一。它能直观展示基因表达、蛋白丰度等组学数据的整体模式，帮助研究者快速发现样本间的相似性和差异性。但在实际分析中，我们经常会遇到一个棘手问题：生物学重复样本的组内重复性较差，导致传统聚类算法无法将同组样本聚在一起。这不仅影响图表美观，更可能误导生物学结论的解读。

ComplexHeatmap是R语言中最强大的热图绘制包之一，由Zuguang Gu博士开发维护。它提供了column_split参数，允许我们基于先验分组信息强制分割热图列，实现"物理聚类"效果。这种方法特别适用于：

• 转录组测序数据中对照组与处理组差异不显著的情况
• 样本采集批次效应明显干扰真实生物学差异的场景
• 临床样本异质性较高导致传统聚类失效的癌症研究
• 需要严格按实验设计展示分组结构的论文图表制作

下面我们将通过完整案例，演示如何利用column_split解决样品重复性差的问题，并分享一系列提升热图专业度的实用技巧。

1. 数据准备与预处理

1.1 数据导入与格式转换

假设我们有一个基因表达矩阵，行代表基因，列代表样本。首先需要将数据转换为适合ComplexHeatmap处理的矩阵格式：

# 加载必要包
library(ComplexHeatmap)
library(circlize)  # 用于颜色映射

# 读取表达矩阵（示例数据）
expr_data <- read.csv("gene_expression.csv", row.names = 1, check.names = FALSE)

# 转换为矩阵格式
expr_matrix <- as.matrix(expr_data)

1.2 数据标准化处理

为了消除基因间表达量级的差异，通常需要对行（基因）进行标准化：

# 对数转换+Z-score标准化
expr_matrix <- t(scale(t(log2(expr_matrix + 1))))

# 检查数据范围
summary(as.vector(expr_matrix))

注意：标准化方法应根据数据类型选择。RNA-seq数据常用log2(CPM+1)或vst转换，而微阵列数据可能只需要log2转换。

1.3 定义样本分组信息

明确样本的分组关系是使用column_split的关键：

# 定义样本分组（假设前3个为对照，后3个为处理）
sample_groups <- factor(
  c(rep("Control", 3), rep("Treatment", 3)),
  levels = c("Control", "Treatment")  # 指定组别顺序
)

2. 基础热图绘制与问题诊断

2.1 传统聚类热图的问题

我们先绘制常规聚类热图，观察样本分组情况：

Heatmap(
  expr_matrix,
  name = "Expression",  # 图例标题
  col = colorRamp2(c(-2, 0, 2), c("blue", "white", "red")),  # 颜色映射
  show_row_names = FALSE,  # 不显示行名
  show_column_names = TRUE,  # 显示列名
  column_names_side = "top",  # 列名位置
  clustering_distance_rows = "pearson",  # 行聚类距离
  clustering_distance_columns = "pearson"  # 列聚类距离
)

如果发现Control和Treatment样本没有很好地分开，说明组内变异可能大于组间差异，传统聚类难以反映真实的生物学分组。

2.2 评估样本重复性

可以通过相关性热图定量评估样本重复性：

# 计算样本间相关性
cor_matrix <- cor(expr_matrix, method = "pearson")

# 绘制相关性热图
Heatmap(
  cor_matrix,
  name = "Pearson\ncorrelation",
  col = colorRamp2(c(0.7, 1), c("white", "red")),
  rect_gp = gpar(col = "black", lwd = 1),  # 单元格边框
  cell_fun = function(j, i, x, y, width, height, fill) {
    grid.text(sprintf("%.2f", cor_matrix[i, j]), x, y, gp = gpar(fontsize = 10))
  }
)

3. 使用column_split强制分组聚类

3.1 基本分组热图实现

当传统聚类效果不佳时，可以使用column_split按预设分组展示：

# 定义顶部注释
top_anno <- HeatmapAnnotation(
  Group = sample_groups,
  col = list(Group = c("Control" = "#1b9e77", "Treatment" = "#d95f02")),
  annotation_name_side = "left"
)

# 绘制分组热图
Heatmap(
  expr_matrix,
  name = "Expression",
  col = colorRamp2(c(-2, 0, 2), c("blue", "white", "red")),
  top_annotation = top_anno,
  column_split = sample_groups,  # 关键参数：按组分割
  show_row_names = FALSE,
  show_column_names = FALSE,
  cluster_columns = FALSE  # 关闭列聚类
)

3.2 分组内再聚类

有时我们希望在保持大分组的前提下，组内样本仍能聚类：

Heatmap(
  expr_matrix,
  name = "Expression",
  col = colorRamp2(c(-2, 0, 2), c("blue", "white", "red")),
  top_annotation = top_anno,
  column_split = sample_groups,
  cluster_column_slices = FALSE,  # 不聚类分组切片
  show_row_names = FALSE,
  show_column_names = FALSE
)

3.3 多级分组处理

对于更复杂的分组设计（如时间序列+处理条件），可以创建多级分组：

# 假设有时间因素（Day1/Day2）
time_points <- factor(rep(c("Day1", "Day2"), each = 3))

# 创建组合分组
combined_groups <- paste(time_points, sample_groups, sep = "_")

# 多级分组热图
Heatmap(
  expr_matrix,
  column_split = list(time_points, sample_groups),
  cluster_column_slices = FALSE,
  top_annotation = HeatmapAnnotation(
    Time = time_points,
    Group = sample_groups,
    col = list(
      Time = c("Day1" = "#7570b3", "Day2" = "#e7298a"),
      Group = c("Control" = "#1b9e77", "Treatment" = "#d95f02")
    )
  )
)

4. 高级美化与注释技巧

4.1 自定义分组标识

使用anno_block可以创建更醒目的分组标识：

top_anno <- HeatmapAnnotation(
  Group = anno_block(
    gp = gpar(fill = c("#1b9e77", "#d95f02")),
    labels = levels(sample_groups),
    labels_gp = gpar(col = "white", fontface = "bold")
  )
)

Heatmap(
  expr_matrix,
  top_annotation = top_anno,
  column_split = sample_groups,
  # 其他参数...
)

4.2 添加行/列注释

丰富的注释可以帮助解读热图：

# 假设有一些基因是差异表达基因
deg_indices <- 1:20  # 示例：前20个基因是DEGs

row_anno <- rowAnnotation(
  DEG = anno_mark(
    at = deg_indices,
    labels = rownames(expr_matrix)[deg_indices],
    side = "left"
  )
)

# 绘制带注释的热图
Heatmap(
  expr_matrix,
  left_annotation = row_anno,
  # 其他参数...
)

4.3 热图组合与布局

ComplexHeatmap支持将多个热图组合展示：

# 假设有另一个数据矩阵（如甲基化数据）
meth_matrix <- matrix(rnorm(nrow(expr_matrix) * ncol(expr_matrix)), nrow = nrow(expr_matrix))

# 绘制组合热图
ht1 <- Heatmap(
  expr_matrix,
  name = "Expression",
  column_split = sample_groups,
  show_column_names = FALSE
)

ht2 <- Heatmap(
  meth_matrix,
  name = "Methylation",
  column_split = sample_groups,
  show_column_names = FALSE
)

ht1 + ht2  # 并排组合

4.4 导出高质量图片

为满足期刊发表要求，需要调整输出参数：

pdf("publication_heatmap.pdf", width = 8, height = 10)
draw(ht_list)  # ht_list为你的热图对象
dev.off()

# 或者保存为高分辨率PNG
png("heatmap.png", width = 2400, height = 3000, res = 300)
draw(ht_list)
dev.off()

提示：对于大型热图（>1000行），建议先保存为PDF或SVG格式，再在矢量图形编辑软件中调整细节。

5. 常见问题与解决方案

5.1 分组顺序调整

如果需要改变分组的显示顺序：

# 重新定义因子水平顺序
sample_groups <- factor(sample_groups, levels = c("Treatment", "Control"))

# 热图将按新顺序显示分组

5.2 处理缺失值

表达矩阵中的NA值需要特殊处理：

# 用行均值填充NA
expr_matrix[is.na(expr_matrix)] <- rowMeans(expr_matrix, na.rm = TRUE)[row(expr_matrix)[is.na(expr_matrix)]]

# 或者直接移除含NA的行
expr_matrix <- expr_matrix[complete.cases(expr_matrix), ]

5.3 大型矩阵优化

对于包含数千基因的表达矩阵，可以采取以下优化措施：

# 只显示差异显著基因
sig_genes <- rownames(expr_matrix)[p.adjust(p_values, method = "fdr") < 0.05]
expr_subset <- expr_matrix[sig_genes, ]

# 或者按方差筛选
gene_vars <- apply(expr_matrix, 1, var)
expr_subset <- expr_matrix[gene_vars > quantile(gene_vars, 0.9), ]

# 绘制子集热图
Heatmap(expr_subset, column_split = sample_groups)

5.4 颜色方案选择

选择合适的颜色方案对数据解读至关重要：

数据类型	推荐颜色方案	适用场景
差异表达	蓝-白-红	展示上下调趋势
相关性矩阵	白-红	强调高相关性区域
离散分类数据	定性调色板（如Set3）	区分不同类别
连续型元数据	渐变色（如viridis）	展示梯度变化

# 使用viridis颜色方案
library(viridis)
Heatmap(
  expr_matrix,
  col = viridis(100),
  column_split = sample_groups
)

6. 实际案例分析

6.1 转录组数据分析

在一个真实的RNA-seq项目中，对照组和处理组各4个生物学重复，但PCA分析显示组内变异较大。使用column_split后，成功在热图中呈现了实验设计结构，同时通过行聚类揭示了基因表达模式。

关键代码调整：

# 调整行聚类参数
Heatmap(
  expr_matrix,
  clustering_distance_rows = "spearman",
  clustering_method_rows = "ward.D2",
  row_dend_width = unit(3, "cm"),
  column_split = sample_groups
)

6.2 临床样本分型研究

在癌症亚型分析中，50个临床样本的异质性导致传统聚类无法清晰展示已知亚型。通过结合临床信息和column_split，热图成功反映了病理分型与分子特征的对应关系。

特殊处理：

# 添加样本临床信息注释
clinical_anno <- HeatmapAnnotation(
  Subtype = clinical_data$Subtype,
  Grade = clinical_data$Grade,
  col = list(
    Subtype = c("LumA" = "pink", "LumB" = "red", "Basal" = "blue"),
    Grade = c("I" = "green", "II" = "yellow", "III" = "orange")
  )
)

Heatmap(
  expr_matrix,
  top_annotation = clinical_anno,
  column_split = clinical_data$Subtype
)

6.3 时间序列实验可视化

一个包含三个时间点（0h, 12h, 24h）的药物处理实验，每个时间点有2个重复。使用多级column_split清晰展示了时间效应和处理效应。

代码实现：

# 定义多级分组
time_group <- paste(experiment_design$Time, experiment_design$Treatment, sep = "_")

Heatmap(
  expr_matrix,
  column_split = list(experiment_design$Time, experiment_design$Treatment),
  cluster_column_slices = FALSE,
  column_title = c("Time", "Treatment", ""),
  column_title_side = "bottom"
)