单细胞亚群相关性分析实战：用corrplot绘制三角热图（附完整代码与避坑指南）

单细胞亚群相关性热图实战：从数据到审稿级可视化，解锁corrplot高阶玩法

在单细胞转录组数据分析的浩瀚海洋里，我们常常需要回答一个核心问题：不同细胞亚群之间，它们的基因表达模式到底有多相似？这种相似性背后，可能隐藏着发育轨迹的线索、细胞状态的连续过渡，或是功能模块的协同表达。而相关性热图，尤其是那些只展示一半信息的三角热图，就成了我们手中那把直观揭示这些关系的“解剖刀”。

但现实往往是，当你兴冲冲地跑完FindAllMarkers，拿到平均表达矩阵，计算出相关性系数后，却发现生成的图形要么色彩平淡无奇，要么布局杂乱无章，距离那些顶刊文章中简洁、信息丰富且极具美感的图表相去甚远。别担心，这篇文章就是为你准备的。我们将绕过那些基础的入门教程，直接深入corrplot包的核心，结合单细胞数据分析的实际场景，手把手带你打造出既能清晰传达科学发现，又能让审稿人眼前一亮的三角相关性热图。我会分享一套经过实战检验的代码模板，并重点剖析那些容易踩坑的参数，让你不仅“知其然”，更“知其所以然”。

1. 数据准备与相关性矩阵计算：构建热图的基石

一切精美的可视化都始于干净、规整的数据。在单细胞语境下，我们的起点通常是一个Seurat对象。假设你已经完成了基础的质控、标准化、降维和聚类，并得到了一个名为seurat_obj的对象，其中seurat_clusters字段存储了细胞的分群信息。

首先，我们需要从亚群层面提取基因表达特征。这里的关键是计算每个亚群在所有基因（或特定标记基因）上的平均表达量。我强烈建议使用差异表达基因（DEGs） 而非全部基因来计算相关性，这能有效过滤噪音，突出有生物学意义的关联。

# 加载必要的R包
library(Seurat)
library(dplyr)
library(corrplot)
library(RColorBrewer)

# 假设你的Seurat对象名为 seurat_obj
# 步骤1：识别各亚群的差异表达基因（这里以寻找正标记基因为例）
Idents(seurat_obj) <- "seurat_clusters" # 确保标识设置为聚类结果
all_markers <- FindAllMarkers(seurat_obj,
                              only.pos = TRUE,
                              min.pct = 0.25,
                              logfc.threshold = 0.25)
# 过滤掉核糖体基因等常见干扰项（可选但推荐）
all_markers <- subset(all_markers, !grepl("^RP[LS]|^MT-", gene))

# 提取在所有比较中均显著的基因（例如，调整后p值 < 0.01）
significant_genes <- all_markers %>%
  group_by(gene) %>%
  filter(any(p_val_adj < 0.01)) %>%
  pull(gene) %>%
  unique()

# 步骤2：计算亚群的平均表达矩阵
# 使用 AverageExpression 函数，针对筛选后的基因
avg_exp <- AverageExpression(seurat_obj,
                             assays = "RNA",
                             features = significant_genes,
                             group.by = "seurat_clusters")
avg_exp_matrix <- avg_exp$RNA # 得到一个 基因 x 亚群 的矩阵

# 步骤3：矩阵转置与标准化（行为基因，列为亚群）
# 为了计算亚群间的相关性，我们需要将矩阵转置，使行为亚群，列为基因
# 同时，通常对基因表达进行Z-score标准化（按行），以消除基因间表达量级差异
gene_exp_for_cor <- t(scale(t(avg_exp_matrix))) # 按基因（行）进行Z-score标准化
# 注意：scale()函数默认对列进行操作，因此我们需要先转置，标准化后再转置回来，或者直接对行操作。
# 上述代码是对行（基因）进行中心化和缩放。

# 更清晰的做法：明确对行操作
# gene_exp_for_cor <- t(apply(avg_exp_matrix, 1, scale))
# rownames(gene_exp_for_cor) <- rownames(avg_exp_matrix)
# colnames(gene_exp_for_cor) <- colnames(avg_exp_matrix)

# 此时，gene_exp_for_cor 是一个（基因 x 亚群）的矩阵，已按基因标准化。
# 计算亚群间的相关性，需要计算这个矩阵列（亚群）之间的相关性。
# 因此，我们使用 cor() 函数，并指定对列计算。
cor_matrix <- cor(gene_exp_for_cor, method = "spearman") # 使用斯皮尔曼秩相关，对离群值更稳健
# method 也可选择 "pearson"（默认，要求数据近似正态）或 "kendall"

# 步骤4：（可选但重要）计算相关性的显著性p值
# corrplot包中的 cor.mtest 函数可以方便地计算
p_matrix