circlize 之 circos.heatmap()

感谢老俊俊的大力支持。我们会每日跟新，欢迎您关注老俊俊的生信笔记。

circlize 之 circos.heatmap()

在上期我们讲到使用 circos.rect()基本函数绘制热图，circlize 还有个 circos.heatmap()高级绘制热图函数，让绘制热图更加方便简单，接下来我们介绍 circos.heatmap()函数。

首先创建一个矩阵：

set.seed(123)
mat1 = rbind(cbind(matrix(rnorm(50*5, mean = 1), nr = 50),
                   matrix(rnorm(50*5, mean = -1), nr = 50)),
             cbind(matrix(rnorm(50*5, mean = -1), nr = 50),
                   matrix(rnorm(50*5, mean = 1), nr = 50))
            )
rownames(mat1) = paste0("R", 1:100)
colnames(mat1) = paste0("C", 1:10)
mat1 = mat1[sample(100, 100), ] # randomly permute rows
split = sample(letters[1:5], 100, replace = TRUE)
split = factor(split, levels = letters[1:5])

绘制标准热图：

library(ComplexHeatmap)
Heatmap(mat1, row_split = split)

1、输入数据

circos.heatmap()的输入应该是一个矩阵(或一个将被转换为单列矩阵的向量)。如果将矩阵划分为多个组，则必须使用 split 参数指定分类变量。注意，split 的值应该是字符向量或因子。如果是数字向量，则在内部将其转换为字符。

必须用定义的颜色模式指定 col 参数。如果矩阵是连续的数字，col 的 value 应该是 colorRamp2()生成的颜色映射，如果矩阵是字符，col 的 value 应该是一个命名的颜色向量。

library(circlize) # >= 0.4.10
# 设置颜色
col_fun1 = colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
# 绘图
circos.heatmap(mat1, split = split, col = col_fun1)
circos.clear()

如果未指定 split，则只有一个热图的大扇区：

circos.heatmap(mat1, col = col_fun1)
circos.clear()

2、环形布局

在绘制圆形图之前可以通过 circos.par()控制圆形布局。

可以在 circos.heatmap()函数中控制热图轨道，如轨迹高度和轨道背景等：

circos.par(start.degree = 90, gap.degree = 10)
circos.heatmap(mat1, split = split, col = col_fun1, track.height = 0.4,
    bg.border = "green", bg.lwd = 2, bg.lty = 2, show.sector.labels = TRUE)
circos.clear()

如果 split 参数的值是一个因子，那么因子级别的顺序将控制热图的顺序。如果 split 是一个简单向量，则 heatmap 的顺序是 unique(split)：

# note since circos.clear() was called in the previous plot,
# now the layout starts from theta = 0 (the first sector is 'e')
circos.heatmap(mat1, split = factor(split, levels = c("e", "d", "c", "b", "a")),
    col = col_fun1, show.sector.labels = TRUE)
circos.clear()

3、聚类树和行名

默认情况下，数值矩阵是按行聚类的，dend.side 参数控制树状图相对于热图轨道的位置。注意，树状图在一个独立的轨道上：

circos.heatmap(mat1, split = split, col = col_fun1, dend.side = "inside")
circos.clear()
circos.heatmap(mat1, split = split, col = col_fun1, dend.side = "outside")
circos.clear()

dend.track.height 参数控制树的高度。rownames.side 参数控制行名的方向：

circos.heatmap(mat1, split = split, col = col_fun1, rownames.side = "inside")
circos.clear()
text(0, 0, 'rownames.side = "inside"')
circos.heatmap(mat1, split = split, col = col_fun1, rownames.side = "outside")
circos.clear()
text(0, 0, 'rownames.side = "outside"')

矩阵的行名和树状图都可以画出来。当然，它们不可能在热图轨道的同一侧：

circos.heatmap(mat1, split = split, col = col_fun1, dend.side = "inside",
    rownames.side = "outside")
circos.clear()
circos.heatmap(mat1, split = split, col = col_fun1, dend.side = "outside",
    rownames.side = "inside")
circos.clear()

行名的图形参数可以设置为长度与矩阵中的行数相同的标量或向量：

circos.heatmap(mat1, split = split, col = col_fun1, rownames.side = "outside",
    rownames.col = 1:nrow(mat1) %% 10 + 1,
    rownames.cex = runif(nrow(mat1), min = 0.3, max = 2),
    rownames.font = 1:nrow(mat1) %% 4 + 1)
circos.clear()

4、聚类

默认情况下，数值矩阵按行聚类。cluster 参数可以设置为 FALSE 来关闭聚类。

当然，当关闭了聚类选项时，即使设置了 dend.side参数，也没有聚类树图：

circos.heatmap(mat1, split = split, cluster = FALSE, col = col_fun1)
circos.clear()

clustering.method 和 distance.method 控制聚类方法和计算距离方法。

注意，circos.heatmap()不支持对列聚类，若需要对列聚类，最好在绘图前对列聚好类：

column_od = hclust(dist(t(mat1)))$order
circos.heatmap(mat1[, column_od])
circos.clear()

5、树图回调

回调函数可以应用到每个在相应扇区生成的树状图上。回调函数编辑树状图，如 1：重新排列树状图，或 2：给系统树图设置颜色。

需要自定义 dend.callback 回调函数，至少包含 3 个参数：

dend：当前扇区的树。
m：当前扇区的矩阵。
si：扇区的所以或者名字。

默认的回调函数定义如下，它通过加权矩阵行平均值来重新排列树状图：

function(dend, m, si) reorder(dend, rowMeans(m))

下面的例子通过 dendsort::dendsort()对每个扇区的树图重新排序：

library(dendsort)
circos.heatmap(mat1, split = split, col = col_fun1, dend.side = "inside",
    dend.callback = function(dend, m, si) {
        dendsort(dend)
    }
)
circos.clear()

我们可以使用 denextend 包中的 color_branches()来渲染树状图的边缘。例如，为五个扇区的树状图分配不同的颜色。这里，树状图轨迹的高度通过 dend.track.height 参数控制：

library(dendextend)
dend_col = structure(1:5, names = letters[1:5])
circos.heatmap(mat1, split = split, col = col_fun1, dend.side = "inside",
    dend.track.height = 0.2,
    dend.callback = function(dend, m, si) {
        # when k = 1, it renders one same color for the whole dendrogram
        color_branches(dend, k = 1, col = dend_col[si])
    }
)
circos.clear()

或者如果矩阵没有被分割，我们可以用不同的颜色给子树状图赋值：

circos.heatmap(mat1, col = col_fun1, dend.side = "inside",
    dend.track.height = 0.2,
    dend.callback = function(dend, m, si) {
        color_branches(dend, k = 4, col = 2:5)
    }
)
circos.clear()

6、多个轨道的热图

如果绘制一个只包含一个热图轨道的环形图，那么使用 circos.heatmap()非常简单。如果你绘制一个包含多个轨道的更复杂的环形图，接下来了解更多的功能。

绘制两层热图：

mat2 = mat1[sample(100, 100), ] # randomly permute mat1 by rows
col_fun2 = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))

circos.heatmap(mat1, split = split, col = col_fun1, dend.side = "outside")
circos.heatmap(mat2, col = col_fun2)
circos.clear()

交换颜色，聚类还是在第一个轨道：

circos.heatmap(mat2, split = split, col = col_fun2, dend.side = "outside")
circos.heatmap(mat1, col = col_fun1)
circos.clear()

对内外两个热图聚类：

circos.heatmap.initialize(mat1, split = split)
circos.heatmap(mat2, col = col_fun2, dend.side = "outside")
circos.heatmap(mat1, col = col_fun1)
circos.clear()

两个热图具有不同列数：

circos.heatmap.initialize(mat1, split = split)
circos.heatmap(mat1[, 1:5], col = col_fun1)
circos.heatmap(mat1[, 6:10], col = col_fun1)
circos.clear()

7、添加其它图形

在完成热图布局之后，可以通过特殊变量 CELL_META 检索用于轨道/扇区/单元格的基本信息。扇区/单元格的数据为如下，它们对于正确对应的热图轨道非常重要：

CELL_META$dend：当前扇区的树信息。
CELL_META$order：当前扇区聚类后的矩阵顺序信息。
CELL_META$subset：完整矩阵的子集，升序排列。

下面是示例圆形热图中第一个扇区的信息输出：

CELL_META$row_dend
## 'dendrogram' with 2 branches and 14 members total, at height 10.51736

CELL_META$row_order
## [1]  2  6  4 12  8  1  5 10  7  9 13 11  3 14

CELL_META$subset
## [1]  8  9 14 18 20 37 55 62 66 72 78 85 93 97

在下面的例子中，我添加了一个轨迹，它可视化了 mat1 中前五列的行表示。我添加了 cell.padding = c(0.02, 0, 0.02, 0)，使最大和最小点不会与单元格的顶部和底部边界重叠：

circos.heatmap(mat1, split = split, col = col_fun1)
row_mean = rowMeans(mat1[, 1:5])
circos.track(ylim = range(row_mean), panel.fun = function(x, y) {
    y = row_mean[CELL_META$subset]
    y = y[CELL_META$row_order]
    circos.lines(CELL_META$cell.xlim, c(0, 0), lty = 2, col = "grey")
    circos.points(seq_along(y) - 0.5, y, col = ifelse(y > 0, "red", "blue"))
}, cell.padding = c(0.02, 0, 0.02, 0))
circos.clear()

类似地，如果将点轨迹作为第一个轨迹放置，则布局应提前初始化：

circos.heatmap.initialize(mat1, split = split)
# This is the same as the previous example
circos.track(ylim = range(row_mean), panel.fun = function(x, y) {
    y = row_mean[CELL_META$subset]
    y = y[CELL_META$row_order]
    circos.lines(CELL_META$cell.xlim, c(0, 0), lty = 2, col = "grey")
    circos.points(seq_along(y) - 0.5, y, col = ifelse(y > 0, "red", "blue"))
}, cell.padding = c(0.02, 0, 0.02, 0))
circos.heatmap(mat1, col = col_fun1) # no need to specify 'split' here
circos.clear()

添加 boxplot：

circos.heatmap(mat1, split = split, col = col_fun1)
circos.track(ylim = range(mat1), panel.fun = function(x, y) {
    m = mat1[CELL_META$subset, 1:5, drop = FALSE]
    m = m[CELL_META$row_order, , drop = FALSE]
    n = nrow(m)
    # circos.boxplot is applied on matrix columns, so here we transpose it.
    circos.boxplot(t(m), pos = 1:n - 0.5, pch = 16, cex = 0.3)
    circos.lines(CELL_META$cell.xlim, c(0, 0), lty = 2, col = "grey")
}, cell.padding = c(0.02, 0, 0.02, 0))
circos.clear()

8、添加注释

扇区的标签可以通过设置 show.sector.labels = TRUE 来添加，但是，这并不提供对标签的任何细节调整。用户可以通过 panel.fun 函数来自定义自己的标签。如下。在这里，标签被添加到距离热图轨道 y 方向上的偏移量 2mm 的地方。

这里设置了 track.index = get.current.track.index()以确保标签总是添加到正确的轨道中：

circos.heatmap(mat1, split = split, col = col_fun1)
circos.track(track.index = get.current.track.index(), panel.fun = function(x, y) {
    circos.text(CELL_META$xcenter, CELL_META$cell.ylim[2] + convert_y(2, "mm"),
        paste0("this is group ", CELL_META$sector.index),
        facing = "bending.inside", cex = 0.8,
        adj = c(0.5, 0), niceFacing = TRUE)
}, bg.border = NA)
circos.clear()

circos.heatmap()不直接支持矩阵的列名，但也可以通过自定义 panel.fun 轻松添加列名。在下面的例子中，设置了更大的空间(10 度，用户通常需要尝试几个值来获得最佳空间)后的最后一个扇区(第五扇区)的间隙。在 circos.par()中的 parameter 之后，然后在 panel.fun 的最后一个扇区中绘制列名：

circos.par(gap.after = c(2, 2, 2, 2, 10))
circos.heatmap(mat1, split = split, col = col_fun1, track.height = 0.4)
circos.track(track.index = get.current.track.index(), panel.fun = function(x, y) {
    if(CELL_META$sector.numeric.index == 5) { # the last sector
        cn = colnames(mat1)
        n = length(cn)
        circos.text(rep(CELL_META$cell.xlim[2], n) + convert_x(1, "mm"),
            1:n - 0.5, cn,
            cex = 0.5, adj = c(0, 0.5), facing = "inside")
    }
}, bg.border = NA)
circos.clear()

给列添加分组：

circos.par(gap.after = c(2, 2, 2, 2, 10))
circos.heatmap(mat1, split = split, col = col_fun1, track.height = 0.4)
circos.track(track.index = get.current.track.index(), panel.fun = function(x, y) {
    if(CELL_META$sector.numeric.index == 5) { # the last sector
        circos.rect(CELL_META$cell.xlim[2] + convert_x(1, "mm"), 0,
                    CELL_META$cell.xlim[2] + convert_x(5, "mm"), 5,
                    col = "orange", border = NA)
        circos.text(CELL_META$cell.xlim[2] + convert_x(3, "mm"), 2.5,
                    "group 1", cex = 0.5, facing = "clockwise")

        circos.rect(CELL_META$cell.xlim[2] + convert_x(1, "mm"), 5,
                    CELL_META$cell.xlim[2] + convert_x(5, "mm"), 10,
                    col = "pink", border = NA)
        circos.text(CELL_META$cell.xlim[2] + convert_x(3, "mm"), 7.5,
                    "group 2", cex = 0.5, facing = "clockwise")
    }
}, bg.border = NA)
circos.clear()

添加图例：

circos.heatmap(mat1, split = split, col = col_fun1)
circos.clear()

library(ComplexHeatmap)
lgd = Legend(title = "mat1", col_fun = col_fun1)
grid.draw(lgd)