新手上路数据挖掘（一）：pheatmap热图

今天我们来学习数据可视化中非常实用的热图绘制方法。就像准备画画需要先准备好画布和颜料一样，绘制热图的第一步，就是要准备好合适的数据。

在开始之前，我们先来了解一下热图绘制需要什么样的数据：

理想情况下，我们需要一个矩阵形式的数据：

每一行代表一个特征（比如基因、产品等） 每一列代表一个样本或观测值（比如不同实验组、不同时间点等） 每个单元格的数值表示该特征在该样本中的测量值

为了让大家更容易上手，这节课我们会使用R语言中自带的airway数据集作为例子。这是一个经过整理的基因表达数据集，非常适合用来学习热图绘制。

📦 airway数据集：RNA-seq入门必备「练手数据」

🔍 数据集身份卡

• • 来源：Bioconductor官方教程常用数据集
• • 本质：RNA测序（RNA-seq）的基因表达矩阵
• • 亮点：体积小、结构清晰、自带元数据

🧬 数据内容

• • 8个样本：4人支气管上皮细胞系 × 2种处理（地塞米松/对照）
• • 64,102个基因：已预处理为原始计数（raw counts）
• • 关键变量：
  • • dex（处理组 vs 对照组）
  • • cell（4种不同细胞系）

先装一下包~

BiocManager::install("airway")

# 加载包
library(airway)
library(DESeq2)
library(pheatmap)


# 加载数据
data("airway")

如果这里运行不报错就说明装成功了。

# 创建DESeqDataSet对象
dds <- DESeqDataSet(airway, design =~ cell + dex)

# 过滤低表达基因
keep <- rowSums(counts(dds))>=
dds <- dds[keep,]

# 运行DESeq2分析
dds <- DESeq(dds)

# 获取标准化后的计数数据
vsd <- vst(dds, blind =FALSE)


# 获取数据矩阵
dat <- assay(vsd)

这时候dat的行是基因，列是样本。

图片

接下来开始绘制热图——

一. 基于样本间距离的热图

可视化样本间的相似性或差异性，快速识别数据中的潜在模式、批次效应或异常样本。 理想状态下同一处理组的样本应彼此距离更近，在热图中聚类在一起。

• • 热图中某个样本远离其所属组的其他样本，可能表明：
  • • 样本处理或测序失败（如低质量RNA）。
  • • 样本标签错误（需检查元数据）。
  • • 真实的生物学异常（如个体差异极大）。

# 计算样本间的距离矩阵
sampleDists <- dist(t(dat))## 注意这里的转置

# 转换为矩阵
sampleDistMatrix <- as.matrix(sampleDists)

# 绘制热图
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         main ="Sample distance matrix")

图片

在计算样本间距离时使用 转置（t()） 的操作是为了确保距离计算是针对 样本（列） 而非 基因（行）。这是基因表达数据分析中的关键步骤，具体原因如下：

1. 数据矩阵的默认结构

在表达矩阵（如RNA-seq的计数矩阵）中，通常：

• • 行（rows）：代表基因或转录本（特征）。
• • 列（columns）：代表样本（观察值）。

例如，assay(vsd) 返回的矩阵格式为： m × n（m个基因 × n个样本）。

2. 距离计算的对象

• • 目标：我们想比较的是 样本之间的相似性（而非基因之间的相似性）。 例如：判断“处理组”样本是否彼此更接近，而与“对照组”样本远离。
• • 问题：**dist() 函数默认对矩阵的 行（rows） 计算距离。**如果直接对 assay(vsd) 计算距离，会得到 基因间的距离（无意义）。
• • 解决：通过转置矩阵（t(assay(vsd))），将样本变为行，基因变为列，使 dist() 正确计算样本间距离。

3. 转置前后的矩阵对比

假设原始矩阵 assay(vsd) 为：

        样本1  样本2  样本3
基因A     10     20     15
基因B     50     30     40
基因C     100    80     90

• • 不转置：dist(assay(vsd)) 计算的是基因A、基因B、基因C之间的距离（无意义）。
• • 转置后：t(assay(vsd)) 变为：

        基因A  基因B  基因C
样本1     10     50    100
样本2     20     30     80
样本3     15     40     90

此时 dist(t(assay(vsd))) 计算的是样本1、样本2、样本3之间的距离（这才是我们需要的）。

4. 其他注意事项

• • 数据标准化：转置前需确保数据已归一化（如 vst 或 rlog 转换），避免高表达基因主导距离计算。
• • 距离度量选择：欧氏距离（dist）适用于标准化后的数据；若用相关性距离，需显式计算（如 1 - cor()）。
• • 热图对称性：样本间距离矩阵是对称的，因此热图的X轴和Y轴标签相同。

代码示例对比

# 错误：计算基因间距离（无意义）
gene_dists <- dist(assay(vsd))

# 正确：计算样本间距离（需转置）
sample_dists <- dist(t(assay(vsd)))

# 相关性距离示例（同样需转置）
cor_dists <- as.dist(- cor(assay(vsd), method="pearson"))# 默认对列计算，无需转置

因此，转置操作的本质是 将数据矩阵从“基因×样本”转换为“样本×基因”，确保距离函数作用于样本而非基因。这是可视化样本聚类（如热图或PCA）前的必要步骤。

二. 基于差异表达基因的热图

基于 差异表达基因（DEGs）的热图 是RNA-seq数据分析中最常用的可视化方法之一，其核心意义在于直观展示基因表达模式在不同样本组间的变化规律，同时揭示潜在的生物学趋势或异常。

# 获取差异表达结果
res <- results(dds)

# 选择显著差异表达的基因（可根据需要调整阈值）
sig_genes <- rownames(res)[which(res$padj <0.05&abs(res$log2FoldChange)>)]

# 如果sig_genes太多，可以只选择前20个
if(length(sig_genes)>){
    sig_genes <- rownames(res)[order(res$padj)][:]
}

# 提取这些基因的表达数据
sig_data <- assay(vsd)[sig_genes,]

# 绘制热图
pheatmap(sig_data,
         scale ="row",# 按行标准化
         show_rownames =TRUE,
         cluster_cols =TRUE,
         main ="Top differentially expressed genes",
         fontsize_row =)

图片

1. 核心意义

(1) 可视化差异表达模式

• • 展示目标：聚焦于显著差异的基因（如padj < 0.05，|log2FC| > 1），用颜色梯度（如红-高表达，蓝-低表达）显示基因在组间的表达变化。
• • 作用：
  • • 快速识别基因的上调或下调趋势。
  • • 验证差异分析结果的可靠性（如热图中组间差异是否与统计结果一致）。

(2) 揭示样本和基因的聚类关系

• • 样本聚类：检查样本是否按实验设计（如处理组/对照组）正确分组。 示例：所有“癌症组”样本在热图中应聚类在一起，并与“正常组”分离。
• • 基因聚类：共表达基因可能具有相似功能或共同调控机制（如通路富集分析时可关联）。

(3) 识别潜在亚群或异常样本

• • 若某个样本未按预期聚类（如对照组样本混入处理组），可能提示：
  • • 样本标签错误、批次效应或技术误差。
  • • 存在未被发现的生物学亚型（如肿瘤异质性）。

2. 热图中关键元素的解读

3. 应用场景示例

(1) 验证差异分析结果

• • 若DESeq2鉴定出“基因A在处理组上调”，热图中应显示该基因在处理组样本中明显偏红（高表达）。

(2) 功能富集线索

• • 聚类在一起的基因可能属于同一通路（如热图中一群基因在癌症组均上调，可后续做KEGG分析）。

(3) 质量控制

• • 热图中若某个样本的基因表达整体偏离同组其他样本，需检查该样本是否为离群值。

4. 绘制热图的注意事项

(1) 数据标准化

• • 必须使用归一化后的表达值（如DESeq2的 rlog 或 vst 转换数据），避免技术偏差影响颜色梯度。
• • 代码示例：# DESeq2的方差稳定变换 vsd <- vst(dds, blind=FALSE)

(2) 基因选择策略

• • 常用方法：
  • • 选择padj显著且log2FC较大的基因（如padj < 0.05 & |log2FC| > 1）。
  • • 按基因表达变化幅度排序（如取log2FC绝对值最高的50个基因）。
• • 代码示例：# 提取显著差异基因 res <- results(dds, contrast=c("dex","treated","untreated")) sig_genes <- subset(res, padj <0.05&abs(log2FoldChange)>) top_genes<- rownames(sig_genes)[order(sig_genes$padj)][:]

三. 添加样本注释的热图

在基于差异表达基因（DEGs）的热图中添加 样本注释（Sample Annotation） 可以显著增强热图的信息量和可解释性。

如何添加样本注释？

图片

注释数据需与热图中的样本顺序一致，且为因子型或字符型变量。

# 创建样本注释数据框
annotation_df<-as.data.frame(colData(vsd)[,c("dex","cell")])

# 绘制带注释的热图
pheatmap(sig_data,# 差异基因的表达矩阵
         scale ="row",# 按基因标准化（Z-score）
         annotation_col = annotation_df,# 列（样本）注释
         show_rownames =TRUE,
         cluster_cols =TRUE,
         main ="DE genes with sample annotations",
         fontsize_row =)

图片

1.关联表达模式与实验变量

• • 作用：在热图顶部或侧边添加颜色条（Color Bar），将样本的元数据（如处理组/对照组、细胞类型、批次、临床分期等）与基因表达模式直接关联。
• • 示例：
  • • 若热图中样本聚类为两大分支，通过注释可立即判断是否对应“处理组 vs 对照组”或“不同细胞系”。
  • • 发现未知协变量（如批次效应）对样本聚类的影响。

2.验证实验设计的合理性

• • 检查主变量（如药物处理）是否主导样本聚类，而其他协变量（如性别、年龄）是否干扰结果。
• • 应用场景：若“处理组”样本未聚类在一起，但注释显示它们属于同一“实验批次”，可能提示批次效应需校正。

3.识别混杂因素或亚群

• • 当样本聚类与预期不符时，注释可揭示潜在混杂因素（如不同测序平台、RNA质量）。
• • 示例：热图中某个亚群可能对应特定的临床特征（如肿瘤分级），提示新的生物学发现。

现在我们顺利生成了三张热图！虽然它们现在还带着点“新手村”的朴素画风，但别小看这第一步——数据可视化的魔法大门已经为你敞开啦！ 🔮