DEG聚类分析热图怎么看？

对于一般的统计分析，基于傻瓜式 *** 作的SPSS（PASW）软件已经足够，但在涉及个性化要求很高的复杂数据处理时，SPSS就开始显得力不从心，这时必须依赖功能更为强大的SAS等软件。以前在自己的科研过程中分析数据多用SPSS、SAS等。在统计遗传和基因组学领域，SAS可以处理很多问题，但与此同时，SAS实现复杂问题过于麻烦，很多问题SAS也不是首选。后来开始运用R环境中的各种免费统计包，特别是Bioconductor的系列分析包，我发觉非常适合生命科学领域的研究者。R有很多优点：
（1）免费，不需要去寻找破解版，不用担心版权问题，使用非常方便；
（2）功能非常强大，单个包的功能比较有限，但多个包组合起来使用则功能无比强大，远胜于SPSS、SAS等；
（3）源代码开放，稍作修改后就能满足个性化的复杂统计分析，满足个性化需求是R的最大特点之一；
（4）程序阅读容易，再加上参考学习资料很多，上手比较容易，提高也不是很难，根据个人经验，要比SAS高级阶段的进阶容易许多；
（5）国际同行高度认同R，我发现很多专用软件都开发了软件的R版，今后R将是数据分析的主流发展方向。
R软件的安装、基本使用等初级教程就不谈了，随便在官方网站找个学习资料就搞定了。“R系列”专辑拟推出中级、高级分析教程。今天推出基因表达谱芯片的聚类分析专题。
本专题示例芯片数据来自GEO数据库中检索号为GSE11787的Affymetrix芯片的CEL文件，共6个CEL文件，3个正常对照组，3个HPS刺激组，为免疫器官脾脏的表达数据。
（一）原始数据的读入、RNA降解评估和标准化
> pd <-readAnnotatedDataFrame("Targettxt",header=TRUE,rownames=1,asis=TRUE)
>rawAffyData <- ReadAffy(filenames=pData(pd)$FileName, phenoData=pd)
> summary(exprs(rawAffyData))
> deg <- AffyRNAdeg(rawAffyData)
> plotAffyRNAdeg(deg,col=c(1,2,3,4,5,6))

> eset <-rma(rawAffyData)
> summary(exprs(eset))

> op <-par(mfrow=c(1,2))
>cols <- brewerpal(6, "Set3")
>boxplot(rawAffyData,col=cols,names=1:6, main ="unnormalizeddata")
>boxplot(dataframe(exprs(eset)) ,names=1:6, main ="normalizationdata", col="blue", border="brown")
>par(op)

（二）聚类分析
原始数据读入，经AffyBatch目标转成ExpressionSet目标后，为提高后续分析（如差异表达基因的检测）的统计功效，往往需要进一步经过Detection CallFilter和IQR filter等过滤（“基因芯片数据的特异性过滤与非特异性过滤”将在另一专题里专门讨论）。
需要说明的是，常规做法是先筛选出差异表达基因，然后只用差异表达基因进行聚类分析（本示例直接用了过滤后的数据集，聚类图的效果稍差一点）。
（1）样本聚类
>dd <-dist2(log2(exprs(eset2)))
>diag(dd) <- 0
>ddrow <- asdendrogram(hclust(asdist(dd)))
>roword <- orderdendrogram(ddrow)
>library("latticeExtra")
>legend <- list(top = list(fun = dendrogramGrob,
args = list(x = ddrow, side = "top")))
>lp <- levelplot(dd[roword, roword],
scales = list(x = list(rot = 90)),
xlab = "", ylab = "", legend = legend)
>plot(lp)

（2）二维聚类
>source(">

这是热图，行是gene，列是RNA-seq样本。这个数据已经通过两种方式进行了修饰，因此我们可以从中获得一些见解。

将行/基因按相似性进行分组。这些基因在样本2中转录最多（在样本4中转录最少)。这些基因在样本1中转录最多(在样本4中最少)。这些基因在样本2中转录最多（在样本3中转录最少)。“聚类”不是偶然的，而是由于一个计算机程序试图把“相似的”东西放在一起。

没有聚类数据会像这样，数据看起来混乱很难去解释。

没有聚类和缩放，将会变成这样，注意到一个基因是高转录的，它是异常值，以至于无法看到其他基因的表达。

如果我们在第一个热图中使用全局缩放会怎样

这一异常值极大地扭曲了缩放，以至于不能看到其他基因。同时，注意到聚类的变化和基因有一个新的顺序。缩放可以影响两件事：

6个样本的RNA-seq的read，

Z值缩放公式 ：

不管原始数据的变化如何，除以 标准偏差 就可以确保数据范围得到缩小。为什么我们要缩小数据的范围，因为我们只能辨别有限颜色的深浅。范围越大，色度的差异就越微妙。通过对数据进行缩放，我们使用的色度更少，更容易看到：“样本1比样本2有更多的转录……”

如果有一个异常值，那个 标准差 将会变非常大，也就是Z值的分母变大，接近于零的值会被压缩到很多，用几个色度很难将它们分开。

当我们使用 异常值 对数据集进行“全局缩放”时，我们看到其中一个基因明显高度表达，但我们看不出其他基因有什么不同。

聚类主要有两种类型：

层次聚类(hierarchical clustering)

热图是很方便的数据可视化方法，很多情景下简单、实用。比如在拿到RNA-seq表达矩阵后大概看一个表达情况，而且可以做个聚类大体了解基因表达模式（当然样本多的时候可以选择WGCNA）。

得到的表就是在输入表的基础上排序加cluster信息。

导出基因顺序参考文章： >

热图是一个以颜色变化来显示数据的可视化矩阵，Toussaint Loua在1873年就曾使用过热图来绘制对巴黎各区的社会学统计。我们就拿这张简单朴素的热图来讲一下热图怎么看。

首先映入我们眼帘的是有的地方是黑的，有的地方是白的（颜色），每一块颜色都有对应的XY轴。言下之意，对象X的属性Y的值是用颜色表征的。颜色的聚集代表相应对象X的属性Y具有相似性（模式，pattern）。本质上它是表现一个数值矩阵，图上每一个小方格都是一个数值，按一条预设好的色彩变化尺（称为色键，Color Key），给每个数值分配颜色。

有时候我们还能看到对象X或者属性Y的聚类结果也绘制在热图的旁边，但是这就不属于热图的部分了，因为他已经不热了（热，就是有的地方冷，有的地方热）。

广泛的应用就是用热图来可视化表达量。我们想象一下一个9个样本50个基因的表达谱，人类一眼看过去就是一堆数字，而表达量数值大小映射到颜色的深浅上，看起来就很清楚了。

很多时候，为了同一个基因在不同样本中的表达量有可比性，需要对表达量取对数，或取Z-score，把数据标准化到一个水平上。

计算两个矩阵的相关性，可以得到两两的相关性，这时，用热图的颜色来表示相关性可以看出哪些配对相关性较高。

这是一张典型的seurat做的热图，可以清楚地看出不同分群有着不同的表达模式。这里的每一个色块是一个细胞某基因的表达量。cluster可以看做是细胞的聚类，Y轴的基因，我们看到也是聚类了的（很可能是手动的，每一类基因作者都给出了注释）。所以这张热图的关键是什么？ 细胞和基因及其顺序。选择合适的细胞和基因（一般是每个群的差异高表达基因）后，为什么我做的图是一团黑？很可能是因为：

人们经常需要根据差异分析的结果来探索基因列表的排序，如SC3的策略。差异基因的计算采用非参数Kruskal-Wallis检验。SC3提供了调整p值< 001的所有差异表达基因的列表，并绘制了p值最低的50个基因的基因表达谱。值得注意的是，聚类后的差异表达计算可能会在p值的分布中引入偏差，因此我们建议仅使用p值对基因进行排序。

这类图无疑反映了 某geneList在某cluster的表达情况 。如果巧了，这个geneList是某个细胞类型的marker基因，或者是某个功能的主要集合，热图有助于细胞群功能和类型的鉴定。热图很好地将对象（X，一般是我们的细胞）与它的属性（Y，一般是我们的基因）联系起来。

在monocle2 中我们还看到一种热图将基因的表达情况与细胞发育轨迹结合到一起。可视化所有明显依赖于分支的基因的变化（如果愿意也可以自己定义geneList）。这张热图同时显示了两种命运的变化，它还要求选择分支点（branch_point ）。列是伪时间中的点，行是基因，伪时间的开始在热图的中间。当你从热图的中间读到右边的时候，你正在跟随一个伪时间谱系。当你读到左边时，另一个。这些基因是分层聚类的，因此您可以可视化具有类似的依赖于序列的基因模块

Answer: In Loupe Cell Browser (version 200), the heatmap is a compact display of a subset of differentially expressed genes per cluster Specifically, the gene list is the union of the top 120/N upregulated genes for each cluster ranked by log2 fold-change (N=total number of clusters)The gene names are on the plot when you export the heatmap However, as of version 200, there is no 1-click function to export the associated information for the subset of heatmap genes

提到相关性，我们很容易注意到WGCNA(weighted correlation network analysis，加权基因共表达网络分析）, 用于提取与性状或临床特征相关的基因模块，解析与表达量相关生物学过程。这是除了富集分析之外另一个寻找好的geneList的方法。这里的颜色不再是表达量的度量而是相似性的度量。

人们针对单细胞发展了相应的数据结构如seurat的S4类，monocle的CDS，SingleCellExperiment的sce,scanpy的anndata等，可见单细胞的故事远大于一张二维的表达谱。那么一张热图往往也不能完全的说明问题，于是我们希望能够灵活地 *** 纵热图来讲更多的故事。于是，我们发现ComplexHeatmap这个R包真的是热图神器。

数据可视化的过程就是一段探索意义的旅程，给每一种颜色、每一种形状、每一种聚集和离散找到一种生物学意义。这让我想起海子的《面朝大海，春暖花开》：

ComplexHeatmap
R数据可视化3：热图
如何画热图
10秒钟-完美掌握-热图（heatmap）绘制 - 所有人都可以！

热图是作为生信中一种常用的数据表现方法，可以简单地聚合大量数据，同时会使用一种相对应的渐变的颜色条带来展示数据。这种方法可以很直观地展现空间数据的疏密程度或频率高低，但由于太直观了，所以热图在数据表现的准确性并不能保证，也就是体现出来一个大概程度。一般情况下的使用场景是展示各种基因或者是RNA在不同的样本中的表达量，然后从整体来观察其表达模式。

在前面的章节中，我们介绍了如何使用 ggplot2 绘制热图

ggplot2 绘制热图的方式很多，如 geom_raster 、 geom_tile 等

但通常仅仅绘制热图是不够的，还需要对数据进行聚类，即绘制聚类热图。

例如，最常用的就是将差异基因的表达值绘制聚类热图，来查看基因在不同样本中的表达差异情况，或者比较不同聚类分组之间的差异。

绘制聚类热图的包有很多，我们主要介绍 pheatmap 和 ComplexHeatmap

假设我们有如下数据

要绘制简单的热图，可以使用内置的 heatmap 函数

更改颜色，并为列添加列样本的分类颜色条

内置函数提供的样式较少，无法对某些图形属性进行设置。

所以下面我们使用 pheatmap 包来绘制热图

pheatmap 对图形属性提供了更精细的控制

这样看起来怪怪的，应该是基因的表达量差异，所以对行进行标准化

嗯，一下子就顺眼多了，实验组和对照组的基因表达量差别明显

默认情况下，会对数据的行列分别进行层次聚类，如果我们想在进行层次聚类之前，先对行特征，也就是基因进行 k-means 聚类，我们可以

先将基因聚为 3 类，再进行层次聚类

如果只想对其中行列中的一个进行聚类，可以使用 cluster_rows 和 cluster_cols 参数，取消对行或列的距离

默认的距离度量为欧氏距离，也可以分别为行列指定不同的距离度量，例如

也可以使用 clustering_method 参数来指定不同的聚类方法，支持以下几种方法：

图例的设置很简单，即通过 legend_breaks 参数设置断点， legend_labels 参数设置断点处的标签

如果不想显示图例，直接设置 legend = F 就行

设置边框颜色

删除边框

默认情况下，单元格的长度和宽度会根据的大小自动调整，如果想固定单元格的大小，可以使用 cellwidth 和 cellheight 两个参数

如果我们想在单元格中显示对于的数值，可以设置 display_numbers = TRUE

对显示的数值进行格式化

或者，为 display_numbers 参数传递一个矩阵

例如，根据表达值是否大于 100 来显示不同的标记

在不对数据进行聚类的情况下，可以对行列进行自定义划分为不同的块

或者只对行或列进行分块

总之，只能对未聚类的行或列进行分块

或者，根据层次聚类的结果，对数据进行分块

使用 main 参数来设置图像的标题

可以使用 show_colnames 和 show_rownames 不显示标签

分别设置标签的大小，同时设置列标签的倾斜角度，可选的角度有 270、0、45、90、315

也可以使用 fontsize 参数统一行列标签的大小

也可以自定义行列标签

我们可以分别为行和列构建分组信息，例如对于行是基因，可以将其分为癌基因和抑癌基因等，而列为样本可以分为癌症和配对正常样本，同时样本对应的患者应该会有年龄性别等信息

例如

我们可以将这些信息以颜色条的方式添加到图中

隐藏图例

我们可以回去 pheatmap 函数返回的对象的信息

可以看到，返回对象 p 中包含 4 个变量，我们可以根据 tree_row 和 tree_col 提取出对应的行列顺序

提取这些信息有助于我们对数据进行分组，用于后续分析

参数列表

数据：
>我们之前在 SC3 : 单细胞转录组聚类分析R包 中介绍过SC3聚类算法，对于细胞异质性未知的样本，聚成几个类我们不知道，这时候可以多尝试几种聚类策略。SC3这个时候就派上用场了。

SC3在可视化上面为我们贡献了几张热图，那么他们都是什么意思呢？

共识矩阵是一个N×N矩阵，其中N是输入数据集中的细胞数。根据所有聚类参数组合的聚类结果的平均值，表示细胞之间的相似性。similarity 0(蓝色)表示这两个细胞总是被分配到不同的群。相反，similarity 1(红色)表示这两个细胞总是被分配到相同的群。共识矩阵采用层次聚类方法进行聚类，具有对角块结构。直观地说，当所有对角块都是完全红色，所有非对角元素都是完全蓝色时，就可以实现完美的聚类。

表达热图表示基因以kmeans进行聚类，k = 100(左侧为树状图)，heatmap表示基因簇中心在log2- scale后的表达水平。

差异基因的计算采用非参数Kruskal-Wallis检验。SC3提供了调整p值< 001的所有差异表达基因的列表，并绘制了p值最低的50个基因的基因表达谱。值得注意的是，聚类后的差异表达计算可能会在p值的分布中引入偏差，因此我们建议仅使用p值对基因进行排序。

为了找到标记基因，每个基因都建立一个基于平均聚类表达值的二元分类器。然后利用基因表达序列计算分类器预测，ROC曲线下的面积用来量化预测的准确性。利用Wilcoxon符号秩检验为每个基因分配一个p值。默认选择ROC曲线下面积(AUROC)为> 085且p值< 001的基因，在此heatmap中显示每个簇的前10个标记基因。

---