ComplexHeatmap包热图绘制进阶版

大家好，我是风。在上一次推文中我们已经使用ComplexHeatmap绘制了简单热图，也就是我们平时在论文中看到的普通热图，并且对简单热图的colors和title的设置技巧进行了展示。对于一个可玩性非常高的R包，难道只有这点能耐吗？当然不是！今天我们继续来看看，ComplexHeatmap还能对我们的热图玩出什么花来呢？
传送门：手把手教你复现一篇CNS级别的美图！附代码，建议收藏！
加载包
# 加载R包library(ComplexHeatmap)library(circlize)# 构建一个跟上次推文一样的矩阵set.seed(123)nr1 = 4; nr2 = 8; nr3 = 6; nr = nr1 + nr2 + nr3nc1 = 6; nc2 = 8; nc3 = 10; nc = nc1 + nc2 + nc3mat <- cbind(rbind(matrix(rnorm(nr1*nc1, mean = 1, sd = 0.5), nr = nr1), matrix(rnorm(nr2*nc1, mean = 0, sd = 0.5), nr = nr2), matrix(rnorm(nr3*nc1, mean = 0, sd = 0.5), nr = nr3)), rbind(matrix(rnorm(nr1*nc2, mean = 0, sd = 0.5), nr = nr1), matrix(rnorm(nr2*nc2, mean = 1, sd = 0.5), nr = nr2), matrix(rnorm(nr3*nc2, mean = 0, sd = 0.5), nr = nr3)), rbind(matrix(rnorm(nr1*nc3, mean = 0.5, sd = 0.5), nr = nr1), matrix(rnorm(nr2*nc3, mean = 0.5, sd = 0.5), nr = nr2),         matrix(rnorm(nr3*nc3, mean = 1,   sd = 0.5), nr = nr3)))mat <- mat[sample(nr, nr), sample(nc, nc)] rownames(mat) = paste0("gene", seq_len(nr))colnames(mat) = paste0("Sample", seq_len(nc))head(mat)## Sample1 Sample2 Sample3 Sample4 Sample5 Sample6## gene1 0.90474160 -0.35229823 0.5016096 1.26769942 0.8251229 0.16215217## gene2 0.90882972 0.79157121 1.0726316 0.01299521 0.1391978 0.46833693## gene3 0.28074668 0.02987497 0.7052595 1.21514235 0.1747267 0.20949120## gene4 0.02729558 0.75810969 0.5333504 -0.49637424 -0.5261114 0.56724357## gene5 -0.32552445 1.03264652 1.1249573 0.66695147 0.4490584 1.04236865## gene6 0.58403269 -0.47373731 0.5452483 0.86824798 -0.1976372 -0.03565404## Sample7 Sample8 Sample9 Sample10 Sample11 Sample12## gene1 -0.2869867 0.68032622 -0.1629658 0.8254537 0.7821773 -0.49625358## gene2 1.2814948 0.38998256 -0.3473535 1.3508922 1.1183375 2.05005447## gene3 -0.6423579 -0.31395304 0.2175907 -0.2973086 0.4322058 -0.25803192## gene4 0.8127096 -0.01427338 1.0844780 0.2426662 0.8783874 1.38452112## gene5 2.6205200 0.75823530 -0.2333277 1.3439584 0.8517620 0.85980233## gene6 -0.3203530 1.05534136 0.7771690 0.4594983 0.2550648 -0.02778098## Sample13 Sample14 Sample15 Sample16 Sample17 Sample18## gene1 -0.0895258 -0.35520328 0.1072694 0.96322199 -0.39245223 -0.1878014## gene2 1.3770269 -0.77437640 0.9829664 0.23854381 -0.53589561 1.3003544## gene3 -0.5686518 -0.51321045 -0.0451598 0.82272880 -0.02251386 0.2427300## gene4 0.8376570 0.10797078 0.4520019 0.81648036 0.02650211 0.4072600## gene5 1.9986067 -0.50928769 0.7708173 -0.09421702 -1.15458444 1.2715970## gene6 -0.2112484 0.01669142 -0.3259750 0.54460361 0.89101254 0.3699738## Sample19 Sample20 Sample21 Sample22 Sample23 Sample24## gene1 1.08736320 0.7132199 -0.1853300 -0.14238650 0.6407669 1.3266288## gene2 0.14237891 0.4471643 0.4475628 -0.31251963 0.7057150 0.8120937## gene3 1.09511516 0.5852612 0.1926402 0.51278568 1.6361334 1.1339175## gene4 -0.02625664 0.9556956 0.3443201 0.07668656 1.7857291 0.5280084## gene5 0.60599022 0.8304101 -0.1902355 -0.14753574 1.0123366 0.2140750## gene6 0.81537706 1.2737905 1.8575325 0.88491126 0.5670193 0.5372756
ComplexHeatmap cluster
看到cluster相信有的学员就开始兴奋了，这不就是机器学习中的聚类嘛！
没错哈，所谓的聚类，简单来讲，就是通过一定的算法把具有相同特征的一堆东西放在一起作为一个簇。聚类可能是热图可视化的关键组成部分,我们看到的热图一般都会使用行聚类或者列聚类，以突出热图中各种特征之间的差异性，从而找到某种规律。比如在差异分析中使用行聚类从而区分哪一区域的分子具有相似的上下调特征。在ComplexHeatmap包中，分层聚类具有非常大的灵活性，包括且不限于基于person、euclidean距离等，我们一起看看：
# 使用参数控制是否进行聚类Heatmap(mat, name = "mat", cluster_rows = TRUE, cluster_columns = FALSE) # 列不聚类，行聚类
# 使用参数控制是否展示聚类树(本质还是进行了聚类，只是不展示聚类树)Heatmap(mat, name = "mat", show_row_dend = TRUE, show_column_dend = FALSE)
# 使用参数调整聚类树方向，并调整聚类树的大小Heatmap(mat, name = "mat", row_dend_side = "right", column_dend_side = "bottom",column_dend_height = unit(2, "cm"),row_dend_width = unit(4, "cm"))
计算距离矩阵
一般情况下，热图的聚类都会分为两步：
1.计算距离矩阵
2.应用聚类方法
在ComplexHeatmap中有三种方法来计算距离，分别是：
（1）将距离指定为预定义选项。即使用dist（）函数以及定义“ pearson”，“ spearman”或“ kendall”，相关距离定义为1-cor（x，y，method），cor函数用来计算相关性，所以相关距离应该是1和相关系数的差值。这种方法允许数据中存在NA值。
（2）一个自定义函数，用于计算矩阵的距离。该函数应仅包含一个参数。但是在列上的聚类将会由矩阵将自动转置。
（3）一个自定义函数，用于计算两个向量之间的距离。该函数应仅包含两个参数。因为这种方法是由两个嵌套的for循环实现的，所以运行速度会比较慢。
接下来我们看一下每种方法都该怎么设置：
# 第一种方法Heatmap(mat, name = "mat", clustering_distance_rows = "pearson",column_title = "第一种方法中的pearson")
# 第二种方法Heatmap(mat, name = "mat", clustering_distance_rows = function(m) dist(m), # 这里的m指代数据mat column_title = "自定义距离矩阵的函数")
# 第三种方法Heatmap(mat, name = "mat", clustering_distance_rows = function(x, y) 1 - cor(x, y),column_title = "自定义函数计算两个向量之间的距离")
基于这些特征，我们可以基于成对距离将鲁棒的聚类应用于离群值。（如果不知道什么叫鲁棒性，请自行百度，如果百度了还不懂，那就等后面有模型类的推文再来一并讲解鲁棒性、稳定性等内容*^_^*）
# 定义矩阵和距离计算函数mat_with_outliers = matfor(i in 1:10) mat_with_outliers[i, i] = 1000# 设置离群值为1000robust_dist = function(x, y) {qx = quantile(x, c(0.1, 0.9)) # 用quantile取百分位比，此处为10%和90% qy = quantile(y, c(0.1, 0.9)) l = x > qx[1] & x < qx[2] & y > qy[1] & y < qy[2]x = x[l]y = y[l]sqrt(sum((x - y)^2)) # sqrt函数取平方根}
我们可以比较使用或不使用鲁棒距离法的两个热图:
# 设置颜色映射，使离群值不影响整体颜色Heatmap(mat_with_outliers, name = "mat", col = colorRamp2(c(-2, 0, 2), c("navy", "white", "red")), column_title = "不使用鲁棒距离法")
Heatmap(mat_with_outliers, name = "mat", col = colorRamp2(c(-2, 0, 2), c("navy", "white", "red")), clustering_distance_rows = robust_dist, clustering_distance_columns = robust_dist, column_title = "使用鲁棒距离法")
由上面可以看出我们可以对矩阵中的距离计算进行自定义，这给了我们很大的操作空间。如果你不知道可以怎么定义函数，也可以使用上面提到的方法，这些方法都是ComplexHeatmap包推荐的方法。当然如果你的矩阵本身就具有某种特征，比如本来就是一个字符矩阵，我们按照stringdist包中的方法对字符矩阵进行聚类：
mat_letters = matrix(sample(LETTERS[1:5], 400, replace = TRUE), 20)# 设置聚类函数dist_letters = function(x, y) { x = strtoi(charToRaw(paste(x, collapse = "")), base = 16) y = strtoi(charToRaw(paste(y, collapse = "")), base = 16)sqrt(sum((x - y)^2))}# 绘制热图Heatmap(mat_letters, name = "letters", col = structure(2:6, names = LETTERS[1:5]), clustering_distance_rows = dist_letters, clustering_distance_columns = dist_letters, cell_fun = function(j, i, x, y, w, h, col) { # 添加文本 grid.text(mat_letters[i, j], x, y) })
应用聚类方法
前面说了，一般热图的聚类都会分为两步：1：计算距离矩阵；2：应用聚类方法。
那么计算完距离之后，下一步就是应用聚类方法了，我们可以使用clustering_method_rows 和 clustering_method_columns 参数进行设置：
Heatmap(mat, name = "mat", clustering_method_rows = "single")
如果已经有一个直接返回集群对象的集群对象或函数，则可以忽略距离设置，并将集群对象或集群函数设置为cluster_rows或cluster_columns。如果它是一个聚类函数，那么唯一的参数应该是矩阵，并且它应该返回一个hclust或dendrogram对象，或者一个可以强制转换为dendrogram类的对象。这种方式可以让我们展示特定聚类方法下的结果，比如，使用来自cluster package的方法进行聚类:
library(cluster)Heatmap(mat,name = "mat", cluster_rows = diana(mat),cluster_columns = agnes(t(mat)), column_title = "cluster package 1")
# 如果cluster_columns被设置为一个函数，则不需要转置矩阵（这种方式类似于自定义函数计算距离）Heatmap(mat,name = "mat", cluster_rows = diana,cluster_columns = agnes, column_title = "cluster package 2_1")
# 上面的代码也可以写成这样Heatmap(mat, name = "mat", cluster_rows = function(m) as.dendrogram(diana(m)), cluster_columns = function(m) as.dendrogram(agnes(m)),column_title = "clusterpackage 2_2")
# 对比上面两个代码，可以发现在cluster_rows的function(m)中的m为数据mat本身，而cluster_columns的function(m)中的m为转置的数据mat，即t(mat)，也可以使用hclust法fh <- function(x) fastcluster::hclust(dist(x))Heatmap(mat, name = "mat", cluster_rows = fh, cluster_columns = fh, column_title = "hclust method")
如果你的行或者列本身就具有特定的聚类标签，比如T分期、stage分期甚至是自定义聚类算法的的cluster，那也可以在聚类树上突出展示这种标签，这样就不需要再多一个图片注释，也可以直接把特定的聚类标签结果反应到热图上。由于示例数据并没有这种特定的聚类标签，因此我们使用算法进行聚类，这里使用IEEE2006年ICDM评选出的数据挖掘的十大算法中排名第二的算法——kmeans算法对样本进行聚类。需要注意的是，在R中，kmeans函数是对行进行聚类，因此我们需要先转置矩阵让矩阵变成行为样本、列为基因，然后再绘制热图：
kmeans_cluster <- kmeans(t(mat), centers = 3, iter.max = 100)kmeans_cluster # 查看kemans结果## K-means clustering with 3 clusters of sizes 10, 8, 6## ## Cluster means:## gene1 gene2 gene3 gene4 gene5 gene6## 1 0.90125516 0.5597888 0.80812370 0.4476205 0.5335800 0.57930227## 2 -0.03974881 1.2003372 -0.16617673 0.6827711 1.3675396 -0.07065958## 3 -0.09300193 -0.1887666 0.01222321 0.2709474 -0.2461226 0.91377635## gene7 gene8 gene9 gene10 gene11 gene12## 1 0.60940007 1.05831848 0.544611273 0.7482195 0.8477304 0.5605398## 2 1.00110216 -0.01285905 -0.009637146 0.9588394 -0.2776678 -0.1445118## 3 -0.01764198 0.32505837 0.979207624 0.1016129 0.1014296 0.8927298## gene13 gene14 gene15 gene16 gene17 gene18## 1 0.764856253 0.98886581 1.17271779 0.3978471 0.5042117 0.5047975## 2 0.008724514 0.08392861 0.09047515 0.8944287 1.1496964 0.8071096## 3 1.196934720 -0.11990428 -0.04449938 0.2487632 0.1568983 -0.1110327## ## Clustering vector:## Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 Sample7 Sample8 ## 1 2 1 1 1 2 2 3 ## Sample9 Sample10 Sample11 Sample12 Sample13 Sample14 Sample15 Sample16 ## 3 2 1 2 2 3 2 1 ## Sample17 Sample18 Sample19 Sample20 Sample21 Sample22 Sample23 Sample24 ## 3 2 1 1 3 3 1 1 ## ## Within cluster sum of squares by cluster:## [1] 43.29525 29.30922 19.22125## (between_SS / total_SS = 47.2 %)## ## Available components:## ## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"## [6] "betweenss" "size" "iter" "ifault"# 对上面kmeans的聚类结果进行解释：第一行表示各个cluster下样本的数量分别是8、10和6个。接着是聚类的均值，即每个基因在各个聚类的中心点。然后是聚类向量，表明每个样本所属的cluster。Within cluster sum of squares by cluster表示每个cluster内部的距离平方和，表示该簇的紧密程度,其中between_SS /total_SS这一项表示组间距离的平方和占整体距离平方和的结果，这个值越接近1越好。最后的Available components表示运行结果返回的对象包含的组成部分group <- kmeans_cluster$clusterHeatmap(mat, name = "mat", cluster_columns = cluster_within_group(mat, group), column_title = "kmeans cluster method")
# 当然也可以像pheatmap一样再打上颜色标签ha <- HeatmapAnnotation(bar = group, col = list(bar = c("1" = "red", "2" = "blue", "3" = "purple")))Heatmap(mat, name = "mat", cluster_columns = cluster_within_group(mat, group), top_annotation = ha, column_title = "kmeans cluster method")
前面我们说的都是在聚类上面做文章，但是聚类树除了cluster，还有dendrograms，接下来我们来看看怎么让我们热图的dendrograms变得与众不同：

可以对聚类树设置颜色映射，一样使用cluster_rows和cluster_columns完成：
library(dendextend)# 根据聚类结果设置颜色row_dend <- as.dendrogram(hclust(dist(mat)))row_dend # 查看聚类结果## 'dendrogram' with 2 branches and 18 members total, at height 6.556313row_dend <- color_branches(row_dend, k = 2) # 根据聚类结果设置k值（颜色）Heatmap(mat, name = "mat", cluster_rows = row_dend)
# 使用row_dend_gp和column_dend_gp设置全局颜色Heatmap(mat, name = "mat",cluster_rows = row_dend, column_dend_gp = gpar(col = c("blue","red")))
使用热图进行可视化，最大的意义应该在于聚类，把相同特性的因素聚类在一起，把不同特性的因素进行区分，展示出因素与因素之间的相似性和差异性，这种因素可以是某种临床特征，也可以是某种表型。在Heatmap()函数中， 如果将row_dend_reorder和column_dend_reorder的值设置为逻辑值的TRUE或者FALSE，则控制热图是否应用树状图重排序；如果将这两个参数设置为数值向量，则它们还控制重新排序的权重（将被映射到reorder.dendrogram（）的wts参数）。
所以我们可以将cluster_rows/cluster_columns设置为逻辑值或聚类函数或者是聚类对象，我们一起看一下：
# 使用参数进行重聚类Heatmap(mat, name = "mat", row_dend_reorder = FALSE, column_dend_reorder = FALSE, column_title = "no reordering")
Heatmap(mat, name = "mat", row_dend_reorder = TRUE, column_dend_reorder = TRUE, column_title = "apply reordering")
# 使用dendsort包进行聚类# install.packages("dendsort")library(dendsort)dend <- dendsort(hclust(dist(mat)))Heatmap(mat, name = "mat", cluster_rows = dend, column_title = "reorder by dendsort", column_dend_gp = gpar(col = c("blue","red")), row_dend_gp = gpar(col = c("blue","red")))
写在最后
在R的学习中，我为什么更喜欢可视化呢？首先肯定是因为可视化可以画出很多赏心悦目的图片啦，当然除此之外，可视化最重要的是提供数据所隐藏的潜在规律，《R数据科学》一书中引用了John Tukey的一句话：
The simple graph has brought more information to the data analyst’s mind than any other device.
                                                                                               ———John Tukey
我想，这也是数据可视化所带来的迷人之处，尤其在大数据时代，我们如何直观地在上百万甚至上千万的数据中找到潜藏的某种规律呢？用的还是可视化（请注意我这里说的是“直观地”）。最后，同样地，喝水不忘挖井人，感谢ComplexHeatmap包的开发者顾祖光博士！这个R包一直是我最爱的可视化R包之一,希望你们学会了后也能喜欢它。

 好了，今天我们先到这里。ComplexHeatmap的故事还将继续，下一期，我们将继续对热图进行修饰，并且结合更多的聚类方法花样玩热图。我们下期再见吧！o(*￣▽￣*)ブ
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文