ComplexHeatmap复杂热图绘制学习——2.单个热图（二

2021-10-20 本文已影响0人 denghb001

2.7热图分割

ComplexHeatmap包的一个主要优点是它支持按行和列拆分热图以更好地对特征进行分组进而使其突出显示。

以下参数控制拆分：row_km, row_split, column_km, column_split。下面，我们将分裂产生的子簇称为"切片"。

2.7.1通过k-means聚类分割

row_km并column_km应用 k-means分割。

Heatmap(mat, name = "mat", row_km = 2)

image

Heatmap(mat, name = "mat", column_km = 3)

image

行拆分和列拆分可以同时进行。

Heatmap(mat, name = "mat", row_km = 2, column_km = 3)

image

注意此时的行和列树状图中有虚线，在后续解释。

Heatmap()内部调用kmeans()具有随机起点，这会导致在某些情况下从重复运行中生成不同的类群。为了解决这个问题，row_km_repeats和column_km_repeats可以设置为一个数字大于1运行kmeans()多次，最终的共有K-均值聚类被使用。请注意，形成共识 k-means 的最终集群数量可能小于row_km和中设置的数量column_km。

# of course, it will be a little bit slower
Heatmap(mat, name = "mat", 
    row_km = 2, row_km_repeats = 100,
    column_km = 3, column_km_repeats = 100)

2.7.2按分类变量拆分

一般来说，row_split或者column_split可以设置为分类向量或数据框，其中不同的级别组合将热图中的行/列分开。如何控制切片的顺序在2.7.4节介绍。

# split by a vector
Heatmap(mat, name = "mat", 
    row_split = rep(c("A", "B"), 9), column_split = rep(c("C", "D"), 12))

image

# split by a data frame
Heatmap(mat, name = "mat", 
    row_split = data.frame(rep(c("A", "B"), 9), rep(c("C", "D"), each = 9)))

image

# split on both dimensions
Heatmap(mat, name = "mat", row_split = factor(rep(c("A", "B"), 9)),
    column_split = factor(rep(c("C", "D"), 12)))

image

实际上，k-means 聚类只是生成一个簇类向量并附加到row_splitor column_split。row_km/column_km可以与row_split和混合使用column_split。

Heatmap(mat, name = "mat", row_split = rep(c("A", "B"), 9), row_km = 2)

image

# code only for demonstration
cl = kmeans(mat, centers = 2)$cluster
# classes from k-means are always put as the first column in `row_split`
Heatmap(mat, name = "mat", row_split = cbind(cl, rep(c("A", "B"), 9)))

如果您对默认的 k-means 分区不满意，只需将分区向量分配给row_split/ column_split即可轻松使用其他分区方法。

pa = cluster::pam(mat, k = 3)
Heatmap(mat, name = "mat", row_split = paste0("pam", pa$clustering))

image

如果设置了row_order或column_order，在每一行/列切片中，它仍然是有序的。

# remember when `row_order` is set, row clustering is turned off
Heatmap(mat, name = "mat", row_order = 18:1, row_km = 2)

image

字符矩阵只能被row_split/column_split参数分割。

# split by the first column in `discrete_mat`
Heatmap(discrete_mat, name = "mat", col = 1:4, row_split = discrete_mat[, 1])

image

如果row_km/column_km被设置或row_split/column_split被设置为向量或数据框，则首先将层次聚类应用于生成k树状图的每个切片，然后基于每个切片的平均值生成父树状图。通过添加所有子切片中树状图的最大高度来调整父树状图的高度，并将父树状图添加到子树状图的顶部以形成单个全局树状图。这就是您在之前的热图中的树状图中看到虚线的原因。它们用于区分父树状图和子树状图，并提醒用户它们的计算方式不同。在Heatmap()这些虚线可以通过设置除去show_parent_dend_line = FALSE，或将其设置为全局选项： ht_opt$show_parent_dend_line = FALSE.

Heatmap(mat, name = "mat", row_km = 2, column_km = 3, show_parent_dend_line = FALSE)

image

2.7.3按树状图分割

拆分的第二种情况是用户可能仍然希望保留从完整矩阵生成的全局树状图，而不是首先将其拆分。在这种情况下，row_split/column_split可以设置为单个数字将应用于cutree()行/列树状图。当cluster_rows/cluster_columns设置为TRUE 或分配有hclust/dendrogram对象时，这会起作用。

对于这种情况，树状图仍然和原来的一样，只是树状图叶子的位置被切片之间的间隙稍微调整了一下。（没有虚线，因为这里树状图被计算为一个完整的树状图，并且没有父树状图或子树状图。）

Heatmap(mat, name = "mat", row_split = 2, column_split = 3)

image

dend = as.dendrogram(hclust(dist(mat)))
dend = color_branches(dend, k = 2)
Heatmap(mat, name = "mat", cluster_rows = dend, row_split = 2)

image

如果你想结合拆分cutree()和其他分类变量，你需要首先生成类cutree()，添加到数据帧，例如 row_split作用数据帧，然后将其发送到 row_split参数。

# code only for demonstration
split = data.frame(cutree(hclust(dist(mat)), k = 2), rep(c("A", "B"), 9))
Heatmap(mat, name = "mat", row_split = split)

2.7.4切片顺序

当row_split/column_split设置为分类变量（向量或数据框）或row_km/column_km设置时，默认情况下，会对切片的均值应用额外的聚类以显示切片级别的层次结构。在这种情况下，您无法精确控制切片的顺序，因为它是由切片的聚类控制的。

不过，你可以设置cluster_row_slices或cluster_column_slices为 FALSE进而关闭聚类上切片，现在你可以精确地控制片的顺序。

当没有切片聚类时，每个切片的顺序可以由levels中row_split/column_split 的每个变量的控制（在这种情况下，每个变量都应该是一个因子）。如果所有变量都是字符，则默认顺序为unique(row_split)orunique(column_split) 。比较以下热图：

Heatmap(mat, name = "mat", 
  row_split = rep(LETTERS[1:3], 6),
    column_split = rep(letters[1:6], 4))

image

# clustering is similar as previous heatmap with branches in some nodes in the dendrogram flipped
Heatmap(mat, name = "mat", 
  row_split = factor(rep(LETTERS[1:3], 6), levels = LETTERS[3:1]),
    column_split = factor(rep(letters[1:6], 4), levels = letters[6:1]))

image

# now the order is exactly what we set
Heatmap(mat, name = "mat", 
  row_split = factor(rep(LETTERS[1:3], 6), levels = LETTERS[3:1]),
    column_split = factor(rep(letters[1:6], 4), levels = letters[6:1]),
    cluster_row_slices = FALSE, 
    cluster_column_slices = FALSE)

image

2.7.5拆分标题

当row_split/column_split被设置为单个数字，只有一个分类变量，而当row_km/column_km和 row_split/column_split同时被设置为分类变量，将存在多个分类变量。默认情况下，标题的形式 "level1,level2,..."对应于所有分类变量中的每个级别组合。拆分的标题可以通过“标题模板”来控制。

ComplexHeatmap支持三种类型的模板。第一个是被相应级别替换的sprintf()地方%s。在下面的例子中，由于所有的组合split都是A,C，A,D，B,C 并且B,D，如果row_title设置为%s|%s，四大行的标题将是 A|C，A|D，B|C，B|D。

split = data.frame(rep(c("A", "B"), 9), rep(c("C", "D"), each = 9))
Heatmap(mat, name = "mat", row_split = split, row_title = "%s|%s")

image

对于sprintf()模板，您只能放置A,B,C,D 的级别在标题中，并且C,D始终在A,B之后（即始终A,C和 A,D）。但是，在制作热图时，您可能希望放置更有意义的文本而不是内部级别。一旦知道如何将文本与级别对应起来，就可以通过以下两种模板方法来添加它。

在以下两个模板方法中，特殊标记用于标记可执行的R代码（称为变量插值，其中提取并执行代码并将返回值放回字符串）。有两种类型的模板标记@{}和{}. 第一个来自 GetoptLong包，在您安装ComplexHeatmap包时应该已经安装，第二个来自您需要先安装的glue包。

x当您使用后两个模板时，您应该使用一个内部变量。x只是一个包含当前类别级别的简单向量（例如 c("A", "C")）。

# We only run the code for the first heatmap
map = c("A" = "aaa", "B" = "bbb", "C" = "333", "D" = "444")
Heatmap(mat, name = "mat", row_split = split, row_title = "@{map[ x[1] ]}|@{map[ x[2] ]}")
Heatmap(mat, name = "mat", row_split = split, row_title = "{map[ x[1] ]}|{map[ x[2] ]}")

image

行标题默认旋转，您可以设置row_title_rot = 0为水平：

Heatmap(mat, name = "mat", row_split = split, row_title = "%s|%s", row_title_rot = 0)

image

当row_split/column_split设置为数字时，您还可以使用模板来调整切片的标题。

Heatmap(mat, name = "mat", row_split = 2, row_title = "cluster_%s")

image

如果你知道最后的行切片数，你可以直接设置一个 titles 向量为row_title。请注意行切片的数量并不总是与nlevel_1*nlevel_2*....

# we know there are four slices
Heatmap(mat, name = "mat", row_split = split, 
    row_title = c("top_slice", "middle_top_slice", "middle_bottom_slice", "bottom_slice"),
    row_title_rot = 0)

image

如果将row_title指定为单个字符串长度，则它就像所有切片的单个标题。

Heatmap(mat, name = "mat", row_split = split, row_title = "there are four slices")

image

如果您仍然想要每个切片的标题，但也需要全局标题，您可以执行以下操作。

ht = Heatmap(mat, name = "mat", row_split = split, row_title = "%s|%s")
# This row_title is actually a heatmap-list-level row title
draw(ht, row_title = "I am a row title")

image

实际上row_title是热图列表draw()函数的行标题，尽管在示例中只有一个热图。该draw() 功能和热图列表将后续介绍。

如果row_title设置为NULL，则不绘制行标题。

Heatmap(mat, name = "mat", row_split = split, row_title = NULL)

image

所有这些规则也适用于列标题切片。

2.7.6分割的图形参数

在行/列上应用拆分时，可以将行/列标题和行/列名称的图形参数指定为与切片数相同的长度。

# by defalt, there no space on the top of the title, here we add 4pt to the top.
# it can be reset by `ht_opt(RESET = TRUE)`
ht_opt$TITLE_PADDING = unit(c(4, 4), "points")
Heatmap(mat, name = "mat", 
    row_km = 2, row_title_gp = gpar(col = c("red", "blue"), font = 1:2),
    row_names_gp = gpar(col = c("green", "orange"), fontsize = c(10, 14)),
    column_km = 3, column_title_gp = gpar(fill = c("red", "blue", "green"), font = 1:3),
    column_names_gp = gpar(col = c("green", "orange", "purple"), fontsize = c(10, 14, 8)))

image

2.7.7切片之间的间隙

行/列切片之间的间隙空间可以通过row_gap/ column_gap控制。该值可以是单个单位或单位向量。

Heatmap(mat, name = "mat", row_km = 3, row_gap = unit(5, "mm"))

image

Heatmap(mat, name = "mat", row_km = 3, row_gap = unit(c(2, 4), "mm"))

image

Heatmap(mat, name = "mat", row_km = 3, row_gap = unit(c(2, 4), "mm"),
    column_km = 3, column_gap = unit(c(2, 4), "mm"))

image

通过设置添加热图边框时，会添加border = TRUE每个切片的边框。

Heatmap(mat, name = "mat", row_km = 2, column_km = 3, border = TRUE)

image

如果您将间隙大小设置为零，则热图看起来就像是由垂直线和水平线分隔的。

Heatmap(mat, name = "mat", row_km = 2, column_km = 3, 
    row_gap = unit(0, "mm"), column_gap = unit(0, "mm"), border = TRUE)

image

2.7.8拆分热图注释

当热图被分割时，所有的热图组件也被相应地分割。以下为您提供了一个简单的例子，热图注释将在第三章介绍。

Heatmap(mat, name = "mat", row_km = 2, column_km = 3,
    top_annotation = HeatmapAnnotation(foo1 = 1:24, bar1 = anno_points(runif(24))),
    right_annotation = rowAnnotation(foo2 = 18:1, bar2 = anno_barplot(runif(18)))
)

image

2.8热图作为光栅图像

当我们要生成所谓的“高品质的图形”，通常我们保存数据，以作为矢量图形的格式为：例如PDF或SVG。矢量图形基本上存储了单个图形元素的细节，因此，如果将一个由非常大的矩阵制成的热图保存为矢量图形，则最终的文件大小会非常大。另一方面，当在屏幕上可视化：例如pdf 文件时，由于屏幕尺寸有限，热图中的多个网格实际上仅映射到单个像素。因此，需要一些方法来有效地缩小原始图像，并且不需要为热图存储完整的矩阵。

光栅化是一种将矢量图形转换为颜色矩阵的方法。在这种情况下，图像表示为 RGB 值的矩阵，称为光栅图像。如果热图大于屏幕的大小或当前图形设备可以支持的像素，我们可以转换热图并缩小它，通过将其保存为与设备具有相同维度的颜色矩阵的形式。

让我们假设一个矩阵有 nr 行和 nc列。当它在某个图形设备上绘制时，例如屏幕设备，相应的热图主体具有pr 和 pc分别用于行和列的像素（或点）。当nr>pr或 nc>pc，矩阵中的多个值被映射到单个像素。这里我们需要减少nr或 nc，如果它们大于pr或 pc。

为了简单起见，我假设两者 nr>pr和 nc>pc. 对于矩阵只有一维比设备大的场景，原理基本相同。
在ComplexHeatmap 2.5.4 版本中，有以下图像光栅化实现。请注意，实现与早期版本略有不同（当然，比早期版本更好）。

首先是一个图像（以特定格式，例如png 或 jpeg）(pr⋅a)×(pc⋅a)分辨率保存在一个临时文件中 a是一个缩放因子，接下来raster通过e.g png::readPNG()或 jpeg::readJPEG()将其作为对象读回，然后通过将光栅对象填充到热图主体中grid::grid.raster()。所以我们可以说，光栅化是由光栅图像设备（png()或jpeg()）完成的。

当行数或列数超过 2000 时（您将看到一条消息。它不会静默地发生），这种类型的光栅化会自动打开（如果没有安装 magick 包）。也可以通过设置use_raster参数来手动控制：

Heatmap(..., use_raster = TRUE)

缩放因子由raster_quality参数控制。大于 1 的值会生成更大的文件。

Heatmap(..., use_raster = TRUE, raster_quality = 5)

只需将原始矩阵减少到pr×pc现在每个单个值都可以对应单个像素。在规定中，应用用户定义的函数来汇总子矩阵。

这可以通过raster_resize_mat参数设置：

# the default summary function is mean()
Heatmap(..., use_raster = TRUE, raster_resize_mat = TRUE)
# use max() as the summary function
Heatmap(..., use_raster = TRUE, raster_resize_mat = max)
# randomly pick one
Heatmap(..., use_raster = TRUE, raster_resize_mat = function(x) sample(x, 1))

具有分辨率的临时图像 nr×nc首先生成，这里 magick::image_resize()用来缩小图片到尺寸pr×pc，将缩小的图像作为raster对象读取并填充到热图主体中。magick提供了很多“resizing”/“scaling”图像的方法，在magick术语下称为“过滤方法” 。所有过滤方法都可以通过获得magick::filter_types()。

这种类型的光栅化可以通过设置raster_by_magick = TRUE 和选择合适的raster_magick_filter.

Heatmap(..., use_raster = TRUE, raster_by_magick = TRUE)
Heatmap(..., use_raster = TRUE, raster_by_magick = TRUE, raster_magick_filter = ...)

临时图像的类型由raster_device参数控制。所有支持“光栅设备”是：png，CairoPNG，agg_png，jpeg，tiff，CairoJPEG和 CairoTIFF。CairoPNG被视为默认光栅设备，因为以前的默认设备png偶尔会在光栅图像中产生白色的垂直和水平线。

在ComplexHeatmap 中，use_raster如果矩阵中的行数或列数超过 2000 ，则默认打开。

在这篇文章的以下部分，我比较了不同图像光栅化方法之间的视觉差异。

下面的例子来自Guillaume Devailly 的模拟数据，但做了小改动。此示例显示了图顶部中心的富集模式。

在下面的例子中，我不会展示制作热图的代码，因为热图太多，具体设置已经写成每个热图的行标题。

我为所有热图设置了相同的颜色映射，以便您可以看到不同的光栅化如何改变原始模式。

为了进行比较，我生成了许多热图。它们可以分为三组，对应于我之前提到的三种光栅化方法。

通过png()/ CairoPNG()/ agg_png()/ jpeg()/ tiff()：光栅化方法1。第一排中的以下热图。
raster_resize_mat = *: 光栅化方法2，有不同的汇总方法。以下热图中的第二行。
filter = *：光栅化方法3，使用不同的filterring方法。字符串filter 应该是raster_magick_filter. 它被截断，以便行标题不会被绘图区域剪切。

image

更多关于图像光栅化的例子可以在这篇博文中找到：ComplexHeatmap 中的光栅化。

根据已显示的示例，我会说通过magick包进行光栅化效果更好，因此，默认情况下，在 ComplexHeatmap 中，光栅化由magick完成（"Lanczos" 作为默认过滤器方法），如果未安装magick，则使用 CairoPNG()并打印一条友好消息，建议用户安装 magick。

以下示例将 PDF 文件大小与不同光栅设备的光栅图像进行比较。

library(GetoptLong)
set.seed(123)
mat2 = matrix(rnorm(10000*100), ncol = 100)
pdf(qq("heatmap_no_rasteration.pdf"), width = 8, height = 8)
ht = Heatmap(mat2, cluster_rows = FALSE, cluster_columns = FALSE, use_raster = FALSE)
draw(ht)
dev.off()

# Here I removed CairoTIFF because it is not working on my laptop
all_devices = c("png", "CairoPNG", "agg_png", "jpeg", "CairoJPEG", "tiff")
for(device in all_devices) {
    pdf(qq("heatmap_@{device}.pdf"), width = 8, height = 8)
    ht = Heatmap(mat2, cluster_rows = FALSE, cluster_columns = FALSE, use_raster = TRUE, 
        raster_device = device)
    draw(ht)
    dev.off()
}

all_files = qq("heatmap_@{all_devices}.pdf", collapse = FALSE)
all_files = c("heatmap_no_rasteration.pdf", all_files)
fs = file.size(all_files)
names(fs) = all_files
sapply(fs, function(x) paste(round(x/1024), "KB"))

## heatmap_no_rasteration.pdf            heatmap_png.pdf 
##                  "6356 KB"                   "175 KB" 
##       heatmap_CairoPNG.pdf        heatmap_agg_png.pdf 
##                   "177 KB"                   "175 KB" 
##           heatmap_jpeg.pdf      heatmap_CairoJPEG.pdf 
##                   "667 KB"                   "677 KB" 
##           heatmap_tiff.pdf 
##                   "175 KB"

CairoPNG 生成较小的文件大小，因此它用作默认光栅设备。

2.9自定义热图主体

热图主体可以自定义添加更多类型的图形。默认情况下，热图主体由具有不同填充颜色的小矩形矩阵（在本文档的其他部分可能称为网格，但在此处称为“cells”）组成。但是，也可以在热图上添加更多图形或符号作为附加层。有两个参数cell_fun，layer_fun它们都应该是用户定义的函数。

2.9.1 cell_fun

cell_fun绘制在每个cells中重复，这是在两个嵌套内部执行for循环，而layer_fun是的向量化版本cell_fun。 cell_fun更容易理解，但layer_fun执行速度更快，更可定制。

cell_fun 期望一个具有 7 个参数的函数（参数名称可以与以下不同，但顺序必须相同），它们是：

j: 矩阵中的列索引。列索引对应于视图中的 x 方向，这j就是将其作为第一个参数的原因。
i: 矩阵中的行索引。
x：在热图主体的视图中测量的单元格中点的 x 坐标。
y：在热图主体的视图中测量的单元格中点的 y 坐标。
width: 单元格的宽度。值unit(1/ncol(sub_mat), "npc") ，其中sub_mat对应于由行分裂和分裂列的子矩阵。
height: 单元格的高度。该值为unit(1/nrow(sub_mat), "npc")。
fill: 单元格的颜色。

在每个单元格中执行时，七个参数的值会自动发送到函数。

最常见的用途是将矩阵中的值添加到热图中：

small_mat = mat[1:9, 1:9]
col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
Heatmap(small_mat, name = "mat", col = col_fun,
    cell_fun = function(j, i, x, y, width, height, fill) {
        grid.text(sprintf("%.1f", small_mat[i, j]), x, y, gp = gpar(fontsize = 10))
})

image

我们也可以选择只为具有正值的单元格添加文本：

Heatmap(small_mat, name = "mat",  col = col_fun,
    cell_fun = function(j, i, x, y, width, height, fill) {
        if(small_mat[i, j] > 0)
            grid.text(sprintf("%.1f", small_mat[i, j]), x, y, gp = gpar(fontsize = 10))
})

image

您可以拆分热图而无需执行任何额外操作cell_fun：

Heatmap(small_mat, name = "mat", col = col_fun,
    row_km = 2, column_km = 2,
    cell_fun = function(j, i, x, y, width, height, fill) {
        grid.text(sprintf("%.1f", small_mat[i, j]), x, y, gp = gpar(fontsize = 10))
})

image

在下面的例子中，我们制作了一个热图，它显示了与corrplot包类似的相关矩阵：

cor_mat = cor(small_mat)
od = hclust(dist(cor_mat))$order
cor_mat = cor_mat[od, od]
nm = rownames(cor_mat)
col_fun = circlize::colorRamp2(c(-1, 0, 1), c("green", "white", "red"))
# `col = col_fun` here is used to generate the legend
Heatmap(cor_mat, name = "correlation", col = col_fun, rect_gp = gpar(type = "none"), 
    cell_fun = function(j, i, x, y, width, height, fill) {
        grid.rect(x = x, y = y, width = width, height = height, 
            gp = gpar(col = "grey", fill = NA))
        if(i == j) {
            grid.text(nm[i], x = x, y = y)
        } else if(i > j) {
            grid.circle(x = x, y = y, r = abs(cor_mat[i, j])/2 * min(unit.c(width, height)), 
                gp = gpar(fill = col_fun(cor_mat[i, j]), col = NA))
        } else {
            grid.text(sprintf("%.1f", cor_mat[i, j]), x, y, gp = gpar(fontsize = 10))
        }
    }, cluster_rows = FALSE, cluster_columns = FALSE,
    show_row_names = FALSE, show_column_names = FALSE)

image

正如您在前面的图中所看到的，设置非标准参数时 rect_gp = gpar(type = "none")，会执行聚类，但在热图主体上没有绘制任何内容。

最后一个例子是可视化一个GO 游戏。输入数据记录游戏中的动作。

str = "B[cp];W[pq];B[dc];W[qd];B[eq];W[od];B[de];W[jc];B[qk];W[qn]
;B[qh];W[ck];B[ci];W[cn];B[hc];W[je];B[jq];W[df];B[ee];W[cf]
;B[ei];W[bc];B[ce];W[be];B[bd];W[cd];B[bf];W[ad];B[bg];W[cc]
;B[eb];W[db];B[ec];W[lq];B[nq];W[jp];B[iq];W[kq];B[pp];W[op]
;B[po];W[oq];B[rp];W[ql];B[oo];W[no];B[pl];W[pm];B[np];W[qq]
;B[om];W[ol];B[pk];W[qp];B[on];W[rm];B[mo];W[nr];B[rl];W[rk]
;B[qm];W[dp];B[dq];W[ql];B[or];W[mp];B[nn];W[mq];B[qm];W[bp]
;B[co];W[ql];B[no];W[pr];B[qm];W[dd];B[pn];W[ed];B[bo];W[eg]
;B[ef];W[dg];B[ge];W[gh];B[gf];W[gg];B[ek];W[ig];B[fd];W[en]
;B[bn];W[ip];B[dm];W[ff];B[cb];W[fe];B[hp];W[ho];B[hq];W[el]
;B[dl];W[fk];B[ej];W[fp];B[go];W[hn];B[fo];W[em];B[dn];W[eo]
;B[gp];W[ib];B[gc];W[pg];B[qg];W[ng];B[qc];W[re];B[pf];W[of]
;B[rc];W[ob];B[ph];W[qo];B[rn];W[mi];B[og];W[oe];B[qe];W[rd]
;B[rf];W[pd];B[gm];W[gl];B[fm];W[fl];B[lj];W[mj];B[lk];W[ro]
;B[hl];W[hk];B[ik];W[dk];B[bi];W[di];B[dj];W[dh];B[hj];W[gj]
;B[li];W[lh];B[kh];W[lg];B[jn];W[do];B[cl];W[ij];B[gk];W[bl]
;B[cm];W[hk];B[jk];W[lo];B[hi];W[hm];B[gk];W[bm];B[cn];W[hk]
;B[il];W[cq];B[bq];W[ii];B[sm];W[jo];B[kn];W[fq];B[ep];W[cj]
;B[bk];W[er];B[cr];W[gr];B[gk];W[fj];B[ko];W[kp];B[hr];W[jr]
;B[nh];W[mh];B[mk];W[bb];B[da];W[jh];B[ic];W[id];B[hb];W[jb]
;B[oj];W[fn];B[fs];W[fr];B[gs];W[es];B[hs];W[gn];B[kr];W[is]
;B[dr];W[fi];B[bj];W[hd];B[gd];W[ln];B[lm];W[oi];B[oh];W[ni]
;B[pi];W[ki];B[kj];W[ji];B[so];W[rq];B[if];W[jf];B[hh];W[hf]
;B[he];W[ie];B[hg];W[ba];B[ca];W[sp];B[im];W[sn];B[rm];W[pe]
;B[qf];W[if];B[hk];W[nj];B[nk];W[lr];B[mn];W[af];B[ag];W[ch]
;B[bh];W[lp];B[ia];W[ja];B[ha];W[sf];B[sg];W[se];B[eh];W[fh]
;B[in];W[ih];B[ae];W[so];B[af]"

我们将其转换为矩阵：

str = gsub("\\n", "", str)
step = strsplit(str, ";")[[1]]
type = gsub("(B|W).*", "\\1", step)
row = gsub("(B|W)\\[(.).\\]", "\\2", step)
column = gsub("(B|W)\\[.(.)\\]", "\\2", step)

go_mat = matrix(nrow = 19, ncol = 19)
rownames(go_mat) = letters[1:19]
colnames(go_mat) = letters[1:19]
for(i in seq_along(row)) {
    go_mat[row[i], column[i]] = type[i]
}
go_mat[1:4, 1:4]

##   a   b   c   d  
## a NA  NA  NA  "W"
## b "W" "W" "W" "B"
## c "B" "B" "W" "W"
## d "B" "W" "B" "W"

根据矩阵中的值放置黑色和白色的石头：

Heatmap(go_mat, name = "go", rect_gp = gpar(type = "none"),
    cell_fun = function(j, i, x, y, w, h, col) {
        grid.rect(x, y, w, h, gp = gpar(fill = "#dcb35c", col = NA))
        if(i == 1) {
            grid.segments(x, y-h*0.5, x, y)
        } else if(i == nrow(go_mat)) {
            grid.segments(x, y, x, y+h*0.5)
        } else {
            grid.segments(x, y-h*0.5, x, y+h*0.5)
        }
        if(j == 1) {
            grid.segments(x, y, x+w*0.5, y)        
        } else if(j == ncol(go_mat)) {
            grid.segments(x-w*0.5, y, x, y)
        } else {
            grid.segments(x-w*0.5, y, x+w*0.5, y)
        }

        if(i %in% c(4, 10, 16) & j %in% c(4, 10, 16)) {
            grid.points(x, y, pch = 16, size = unit(2, "mm"))
        }

        r = min(unit.c(w, h))*0.45
        if(is.na(go_mat[i, j])) {
        } else if(go_mat[i, j] == "W") {
            grid.circle(x, y, r, gp = gpar(fill = "white", col = "white"))
        } else if(go_mat[i, j] == "B") {
            grid.circle(x, y, r, gp = gpar(fill = "black", col = "black"))
        }
    },
    col = c("B" = "black", "W" = "white"),
    show_row_names = FALSE, show_column_names = FALSE,
    column_title = "One famous GO game",
    heatmap_legend_param = list(title = "Player", at = c("B", "W"), 
        labels = c("player1", "player2"), border = "black")
)

image

2.9.2 layer_fun

cell_fun逐个单元添加图形，而layer_fun以块方式添加图形。与cell_fun,类似，layer_fun也需要七个参数，但它们都是向量形式（layer_fun也可以有第八个和第九个参数，这将在本节后面介绍）：

# code only for demonstration
Heatmap(..., layer_fun = function(j, i, x, y, w, h, fill) {...})
# or you can capitalize the arguments to mark they are vectors,
# the names of the argumetn do not matter
Heatmap(..., layer_fun = function(J, I, X, Y, W, H, F) {...})

j并且i仍然包含与原始矩阵对应的列和行索引，但因为现在layer_fun适用于一个单元格块（如果热图被分割，则为一个热图块），j并且i是当前热图切片中所有单元格的向量。类似y，x，w，h和fill是对应于在当前热图片中的所有小区中的所有载体。

由于j和i现在是向量，为了获得矩阵中的相应值，我们不能使用该形式，mat[j, i]因为它为您提供了一个包含length(i)行和length(j)列的子矩阵。相反，我们可以使用ComplexHeatmap中的pindex() 函数，这类似于矩阵的成对索引。请参阅以下示例：

mfoo = matrix(1:9, nr = 3)
mfoo[1:2, c(1, 3)]

##      [,1] [,2]
## [1,]    1    7
## [2,]    2    8

# but we actually want mfoo[1, 1] and mfoo[2, 3]
pindex(mfoo, 1:2, c(1, 3))

## [1] 1 8

下一个示例显示了layer_fun在热图上添加文本的版本。与cell_fun版本基本一致。

col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
Heatmap(small_mat, name = "mat", col = col_fun,
    layer_fun = function(j, i, x, y, width, height, fill) {
        # since grid.text can also be vectorized
        grid.text(sprintf("%.1f", pindex(small_mat, i, j)), x, y, gp = gpar(fontsize = 10))
})

image

并且只向具有正值的单元格添加文本：

Heatmap(small_mat, name = "mat", col = col_fun, 
    layer_fun = function(j, i, x, y, width, height, fill) {
        v = pindex(small_mat, i, j)
        l = v > 0
        grid.text(sprintf("%.1f", v[l]), x[l], y[l], gp = gpar(fontsize = 10))
})

image

当热图被分割时，layer_fun应用于每个切片。

Heatmap(small_mat, name = "mat", col = col_fun,
    row_km = 2, column_km = 2,
    layer_fun = function(j, i, x, y, width, height, fill) {
        v = pindex(small_mat, i, j)
        grid.text(sprintf("%.1f", v), x, y, gp = gpar(fontsize = 10))
        if(sum(v > 0)/length(v) > 0.75) {
            grid.rect(gp = gpar(lwd = 2, fill = "transparent"))
        }
})

image

layer_fun还可以有另外两个参数，它们是当前行切片和列切片的索引。例如，我们想为右上角和左下角的切片添加边框。

Heatmap(small_mat, name = "mat", col = col_fun,
    row_km = 2, column_km = 2,
    layer_fun = function(j, i, x, y, width, height, fill, slice_r, slice_c) {
        v = pindex(small_mat, i, j)
        grid.text(sprintf("%.1f", v), x, y, gp = gpar(fontsize = 10))
        if(slice_r != slice_c) {
            grid.rect(gp = gpar(lwd = 2, fill = "transparent"))
        }
})

image

这样做的好处layer_fun是不仅可以快速添加图形，还可以提供更多自定义热图的可能性。考虑以下可视化：对于热图中的每一行，如果相邻两列中的值具有相同的符号，我们根据两个值的符号添加一条红线或一条绿线。由于现在单元格中的图形依赖于其他单元格，因此只能通过layer_fun. （不要被下面的代码吓到。代码后面会解释。）

Heatmap(small_mat, name = "mat", col = col_fun,
    row_km = 2, column_km = 2,
    layer_fun = function(j, i, x, y, w, h, fill) {
        # restore_matrix() is explained after this chunk of code
        ind_mat = restore_matrix(j, i, x, y)
        for(ir in seq_len(nrow(ind_mat))) {
            # start from the second column
            for(ic in seq_len(ncol(ind_mat))[-1]) {
                ind1 = ind_mat[ir, ic-1] # previous column
                ind2 = ind_mat[ir, ic]   # current column
                v1 = small_mat[i[ind1], j[ind1]]
                v2 = small_mat[i[ind2], j[ind2]]
                if(v1 * v2 > 0) { # if they have the same sign
                    col = ifelse(v1 > 0, "darkred", "darkgreen")
                    grid.segments(x[ind1], y[ind1], x[ind2], y[ind2],
                        gp = gpar(col = col, lwd = 2))
                    grid.points(x[c(ind1, ind2)], y[c(ind1, ind2)], 
                        pch = 16, gp = gpar(col = col), size = unit(4, "mm"))
                }
            }
        }
    }
)

image

发送到的值layer_fun都是向量（用于grid图形函数的向量化），但是，layer_fun应用到的热图切片仍然由矩阵表示，因此，如果所有参数都layer_fun可以，这将非常方便转换为当前切片的子矩阵。在这里，如上例所示， restore_matrix()完成这项工作。restore_matrix()直接接受in的前四个参数layer_fun并返回一个索引矩阵，其中行和列对应当前切片中的行和列，从上到下，从左到右。矩阵中的值是j当前切片中向量的自然顺序。

如果您运行以下代码：

Heatmap(small_mat, name = "mat", col = col_fun,
    row_km = 2, column_km = 2,
    layer_fun = function(j, i, x, y, w, h, fill) {
        ind_mat = restore_matrix(j, i, x, y)
        print(ind_mat)
    }
)

左上角切片的第一个输出：

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    4    7   10   13
[2,]    2    5    8   11   14
[3,]    3    6    9   12   15

如您所见，这是一个三行五列的索引矩阵，其中第一行对应于切片中的顶行。矩阵中的值对应于 in j, i, x, y, ... in的自然索引（即 1, 2, ...）layer_fun。现在，如果我们想在左上角切片的第二列上添加值，放在里面的代码layer_fun看起来像：

for(ind in ind_mat[, 2]) {
    grid.text(small_mat[i[ind], j[ind]], x[ind], y[ind], ...)
}

现在更容易理解第二个例子：我们想在每个切片的第二行和第三列添加点：

Heatmap(small_mat, name = "mat", col = col_fun,
    row_km = 2, column_km = 2,
    layer_fun = function(j, i, x, y, w, h, fill) {
        ind_mat = restore_matrix(j, i, x, y)
        ind = unique(c(ind_mat[2, ], ind_mat[, 3]))
        grid.points(x[ind], y[ind], pch = 16, size = unit(4, "mm"))
    }
)

image

2.10热图的大小

width, heatmap_width,height并heatmap_height控制热图的大小。默认情况下，所有热图组件都有固定的宽度或高度，例如行树状图的宽度为1cm。热图主体的宽度或高度填充最终绘图区域的其余区域，这意味着，如果您在交互式图形窗口中绘制它并通过拖动它来更改窗口的大小，则热图主体的大小为自动调整。

heatmap_width和heatmap_height控制整个热图包括所有热图组件（不包括图例）的宽度/高度，同时width 与height仅控制heamtap主体的宽度/高度。所有这四个参数都可以设置为绝对单位。

Heatmap(mat, name = "mat", width = unit(8, "cm"), height = unit(8, "cm"))

image

Heatmap(mat, name = "mat", heatmap_width = unit(8, "cm"), heatmap_height = unit(8, "cm"))

image

在调整热图列表中的大小时，这四个参数更为重要（参见第4.2节）。

当热图的大小设置为绝对单位时，图的大小可能大于热图的大小，这会在热图周围产生空白区域。热图的大小可以通过ComplexHeatmap:::width() 和ComplexHeatmap:::height()函数检索。

ht = Heatmap(mat, name = "mat", width = unit(8, "cm"), height = unit(8, "cm"))
ht = draw(ht) # you must call draw() to reassign the heatmap variable
ComplexHeatmap:::width(ht)

## [1] 118.851591171994mm

ComplexHeatmap:::height(ht)

## [1] 114.717557838661mm

ht = Heatmap(mat, name = "mat", heatmap_width = unit(8, "cm"), 
    heatmap_height = unit(8, "cm"))
ht = draw(ht)
ComplexHeatmap:::width(ht)

## [1] 94.8877245053272mm

ComplexHeatmap:::height(ht)

## [1] 83.8660578386606mm

如果要将热图保存在与热图大小完全相同的图形文件中，请查看此博客文章（在热图中设置单元格宽度/高度）。

2.11绘制热图

Heatmap()函数实际上只是一个构造函数，也就是说它只将所有的数据和配置放入Heatmap类中的对象中。只有在draw()调用该方法时才会执行聚类。在交互模式下（例如，您可以在交互式 R 终端中逐行键入 R 代码），直接调用Heatmap()而不返回任何对象会打印对象和类对象内部调用的打印方法（或 S4show()方法）。因此，如果您在 R 终端中输入，它看起来像是一个绘图函数，例如，您需要注意它实际上不是真的，并且在以下情况下您可能看不到任何绘图。Heatmap``draw()``Heatmap(...)``plot()

你将Heatmap(...)放入一个函数中，
你将Heatmap(...)放入一个代码块，如for或if-else
你将Heatmap(...)放入一个 Rscript 并在命令行下运行它。

原因是在以上三种情况下，show()方法不会被调用，因此draw()方法也不会被执行。因此，要绘制该图，您需要draw()显式调用：draw(Heatmap(...))或：

# code only for demonstration
ht = Heatmap(...)
draw(ht)

该draw()函数实际上应用于HeatmapList类中的热图列表。draw()单个Heatmap类的方法构造一个HeatmapList只有一个热图和类的调用draw()方法HeatmapList。该draw()函数接受更多参数，例如控制图例。它将在第4章中讨论。

draw(ht, heatmap_legend_side, padding, ...)

2.12获取顺序和树状图

热图的行/列顺序可以通过row_order()/column_order()函数获得。您可以直接应用于由返回的热图对象Heatmap()或由返回的对象draw()。下面，我们以row_order()为例。

small_mat = mat[1:9, 1:9]
ht1 = Heatmap(small_mat)
row_order(ht1)

## [1] 5 7 2 4 9 6 8 1 3

ht2 = draw(ht1)
row_order(ht2)

## [1] 5 7 2 4 9 6 8 1 3

如上一节所述，Heatmap()function 不执行聚类，因此，当直接 apply row_order()on 时ht1，将首先执行聚类。稍后在通过制作热图时draw(ht1)，将再次应用聚类。如果您在热图中设置 k-means 聚类，这可能是一个问题。由于聚类应用了两次，k-means 可能会给你不同的聚类，这意味着你可能有不同的结果row_order()，你可能有不同的热图。

在以下代码块中o1，o2和o3可能会有所不同，因为每次都执行 k 均值聚类。

# code only for demonstration
ht1 = Heatmap(small_mat, row_km = 2)
o1 = row_order(ht1)
o2 = row_order(ht1)
ht2 = draw(ht1)
o3 = row_order(ht2)
o4 = row_order(ht2)

draw()函数返回已经重新排序的热图（或更准确地说，热图列表），并且应用row_order()只是从对象中提取行顺序，这确保行顺序与热图中显示的行顺序完全相同。在上面的代码中，o3始终与o4.

因此，获取行/列顺序的优选方法如下。

# code only for demonstration
ht = Heatmap(small_mat)
ht = draw(ht)
row_order(ht)
column_order(ht)

如果行/列被拆分，行顺序或列顺序将是一个列表。

ht = Heatmap(small_mat, row_km = 2, column_km = 3)
ht = draw(ht)
row_order(ht)

## $`2`
## [1] 5 7 2 4
## 
## $`1`
## [1] 9 6 8 1 3

column_order(ht)

## $`1`
## [1] 6 2 7
## 
## $`3`
## [1] 1 3 4 5
## 
## $`2`
## [1] 8 9

同样，row_dend()/column_dend()函数返回树状图。它返回单个树状图或树状图列表，具体取决于热图是否被分割。

ht = Heatmap(small_mat, row_km = 2)
ht = draw(ht)
row_dend(ht)

## $`2`
## 'dendrogram' with 2 branches and 4 members total, at height 2.946428 
## 
## $`1`
## 'dendrogram' with 2 branches and 5 members total, at height 2.681351

column_dend(ht)

## 'dendrogram' with 2 branches and 9 members total, at height 4.574114

row_order(), column_order(),row_dend()和column_dend()也适用于热图列表，它将在第4.12节中介绍。

2.13热图子集

由于热图是矩阵的表示，因此Heatmap该类还有一个子集方法。

ht = Heatmap(mat, name = "mat")
dim(ht)

## [1] 18 24

ht[1:10, 1:10]

image

注释也相应地进行了子集化。

ht = Heatmap(mat, name = "mat", row_km = 2, column_km = 3,
    col = colorRamp2(c(-2, 0, 2), c("green", "white", "red")),
    top_annotation = HeatmapAnnotation(foo1 = 1:24, bar1 = anno_points(runif(24))),
    right_annotation = rowAnnotation(foo2 = 18:1, bar2 = anno_barplot(runif(18)))
)
ht[1:9*2 - 1, 1:12*2] # odd rows, even columns

image

如果热图组件类似于向量，则它们被子集化。热图中的一些配置在子化时保持不变，例如如果row_km 在原始热图中设置，则保留k-means的配置并在子热图中执行。所以在下面的例子中，k-means 聚类只在为制作热图时执行ht2。

ht = Heatmap(mat, name = "mat", row_km = 2)
ht2 = ht[1:10, 1:10]
ht2

image

子集热图的实现是非常实验性的。它并不总是有效，例如如果cell_fun已定义并使用外部矩阵，或者将聚类对象分配给cluster_rows或cluster_columns。

HeatmapAnnotation类（第3.19节）和HeatmapList类（第4.10节）也有子集方法，但两者都非常具有实验性。

ComplexHeatmap复杂热图绘制学习——2.单个热图（二

2.7热图分割

2.7.1通过k-means聚类分割

2.7.2按分类变量拆分

2.7.3按树状图分割

2.7.4切片顺序

2.7.5拆分标题

2.7.6分割的图形参数

2.7.7切片之间的间隙

2.7.8拆分热图注释

2.8热图作为光栅图像

2.9自定义热图主体

2.9.1 cell_fun

2.9.2 layer_fun

2.10热图的大小

2.11绘制热图

2.12获取顺序和树状图

2.13热图子集

猜你喜欢

热点阅读