用 ComplexHeatmap 包绘制复杂热图

生信菜鸟团

发布于 2021-01-18 14:38:33

6.8K0

发布于 2021-01-18 14:38:33

最近碰到个画热图的需求，以前一直用的 pheatmap，但这次的图有些复杂，靠我的水平用 pheatmap 还是做不出来。早已听师兄推荐过 ComplexHeatmap 包的强大之处，今天有空就把他学了一下（虽然在前两天我已经用 ggplot 和 Y 叔的 aplot 搞定了）。

安装

if (!requireNamespace("BiocManager", quietly=TRUE))    install.packages("BiocManager")BiocManager::install("ComplexHeatmap")

热图组成

在 ComplexHeatmap 中单个热图由热图主体和热图组件组成。热图主体可按行或列进行拆分。热图组件包括标题，进化树，矩阵名称和热图注释，可分别放置于热图主体的四个侧面上，这些组件也可根据热图主体的顺序进行重新排序或拆分。

除了画单个热图之外，ComplexHeatmap 还支持组合多个热图，即称之为热图列表 heatmap list ，一系列热图和热图注释的集合。在热图列表周围，可设置全局级别的标题和图例。

当然除了横向排列的热图列表外，还可以纵向排列。

ComplexHeatmap 包以面向对象的方式实现。为了描述热图列表，主要有以下几类：

•Heatmap 类：单个热图，其中包含热图主体，行/列名称，标题，进化树和行/列注释。•HeatmapList 类：热图和热图注释的列表。•HeatmapAnnotation 类：定义行注释和列注释的列表。热图注释可以是热图的组成部分，也可以独立于热图。

还有一些内部类：

•SingleAnnotation 类：定义单个行注释或列注释。HeatmapAnnotation 对象其实就是 SingleAnnotation 对象的列表。•ColorMapping 类：从值到颜色的映射。主矩阵和注释的颜色映射由 ColorMapping 类控制。•AnnotationFunction 类：构建用户自定义的注释。

ComplexHeatmap 通过 grid 绘图系统实现，所以需要对基本的 grid 函数有一定了解才能充分利用这个包。

单个热图

先准备示例数据，生成一个随机矩阵，其中列分为三组，行分为三组：

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)] # random shuffle rows and columnsrownames(mat) = paste0("row", seq_len(nr))colnames(mat) = paste0("column", seq_len(nc))

使用基础函数 Heatmap() ，一行代码即可绘制默认参数下的热图：

Heatmap(mat)

图例的标题默认为热图的“名称”。若名称没有设定则默认为 matrix_%d 格式，如上图中的 matrix_1

颜色

可使用 circlize::colorRamp2() 函数在 Heatmap() 中生成颜色映射函数。 colorRamp2() 的两个参数分别是值的向量和对应颜色的向量。在下面的示例中，我们将根据 -2 到 2 的区域生成对应颜色，大于 2 的值都映射为红色，小于 -2 的值都映射为绿色（所以热图的颜色不会受异常值的影响）。

library(circlize)col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))col_fun(seq(-3, 3))

## [1] "#00FF00FF" "#00FF00FF" "#B1FF9AFF" "#FFFFFFFF" "#FF9E81FF" "#FF0000FF"## [7] "#FF0000FF"

# 这里用 name 为热图指定了名称，可以看到图中图例的名称改为了 nameHeatmap(mat, name = "mat", col = col_fun)

对于包含 NA 的矩阵，我们可用 na_col 设置对应颜色：

mat_with_na = matna_index = sample(c(TRUE, FALSE), nrow(mat)*ncol(mat), replace = TRUE, prob = c(1, 9))mat_with_na[na_index] = NAHeatmap(mat_with_na, name = "mat", na_col = "black",    column_title = "a matrix with NA values")

可以通过 border/border_gp 和 rect_gp 参数设置热图边框的颜色。

•border/border_gp 控制热图主体的全局边框。border 的值可以是逻辑值（TRUE 对应黑色）或颜色字符（例如 red）。border 参数的使用仅出于历史原因，在这里我们还可以设置 border_gp 参数，该参数应该是 gpar 对象。•rect_gp 控制热图中网格/单元的边框。rect_gp 是一个 gpar 对象。

Heatmap(mat, name = "mat", border_gp = gpar(col = "black", lty = 2),    column_title = "set heatmap borders")

Heatmap(mat, name = "mat", rect_gp = gpar(col = "white", lwd = 2),    column_title = "set cell borders")

标题

可使用 column_title 和 row_title 分别指定列名和行名：

Heatmap(mat, name = "mat", column_title = "I am a column title",     row_title = "I am a row title")

另外可用 column_title_side 参数指定标题位置：

Heatmap(mat, name = "mat", column_title = "I am a column title at the bottom",     column_title_side = "bottom")

除此之外，可用 row_title_gp 或 column_title_gp 实现更多样的标题样式。

Heatmap(mat, name = "mat", column_title = "I am a big column title",     column_title_gp = gpar(fontsize = 20, fontface = "bold"))

Heatmap(mat, name = "mat", column_title = "I am a column title",     column_title_gp = gpar(fill = "red", col = "white", border = "blue"))

聚类

ComplexHeatmap 支持：

•预定义的距离算法（例如 "euclidean" 或 "pearson"），•其他距离算法•一个已包含聚类信息的对象（例如 hclust 或 dendrogram 对象）•其他聚类算法

同时还可以为不同的节点和分支渲染具有不同颜色和样式，以更好地展示进化树。

首先是一些常规设置

•关闭对列的聚类

Heatmap(mat, name = "mat", cluster_rows = FALSE) # turn off row clustering

•隐藏纵向的进化树

Heatmap(mat, name = "mat", show_column_dend = FALSE) # hide column dendrogram

•设置进化树位置：

Heatmap(mat, name = "mat", row_dend_side = "right", column_dend_side = "bottom")

•设置进化树高度

Heatmap(mat, name = "mat", column_dend_height = unit(4, "cm"),     row_dend_width = unit(4, "cm"))

•可用 clustering_distance_rows = "pearson" 指定聚类算法

除了用算法聚类外，我们仍然可以按手动设置行和列的顺序。如果例如设置row_order，这样默认情况下会关闭行聚类。

Heatmap(mat, name = "mat", row_order = order(as.numeric(gsub("row", "", rownames(mat)))),     column_order = order(as.numeric(gsub("column", "", colnames(mat)))),    column_title = "reorder matrix")

分割热图

ComplexHeatmap 支持按行和列拆分热图，以更好地对功能进行分组并突出显示模式。包括 row_km，row_split，column_km，column_split 这四个参数。

•使用 row_km 和 column_km 参数根据 k-means 聚类拆分行或列：

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

我们还可以手动设置分类并分割热图。

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

•使用 row_split，column_split 根据进化树拆分热图：

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

用 row_gap 或 column_gap 设置分割空隙大小

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

热图大小

heatmap_width 和 heatmap_height 分别控制整幅热图的宽和高（包括所有热图组件）而 width 和 height 参数仅控制热图主题区域的宽和高。

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

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

热图注释

ComplexHeatmap 包还具有强大的可扩展性，我们可以通过 top_annotation，bottom_annotation，left_annotation 和 right_annotation 参数为热图的四个侧面添加注释。

热图注释的简单用法如下：

set.seed(123)mat = matrix(rnorm(100), 10)rownames(mat) = paste0("R", 1:10)colnames(mat) = paste0("C", 1:10)column_ha = HeatmapAnnotation(foo1 = runif(10), bar1 = anno_barplot(runif(10)))row_ha = rowAnnotation(foo2 = runif(10), bar2 = anno_barplot(runif(10)))Heatmap(mat, name = "mat", top_annotation = column_ha, right_annotation = row_ha)

添加区块注释以区分不同分组：

Heatmap(matrix(rnorm(100), 10), name = "mat",    top_annotation = HeatmapAnnotation(foo = anno_block(gp = gpar(fill = 2:4))),    column_km = 3)

绘制更复杂的热图

为基因表达矩阵加入更多信息

热图常用于可视化基因表达矩阵，矩阵中的行与基因相对应，我们可以在表达热图后附加有关这些基因的更多信息。

在下面的示例中，大的热图展示了基因的表达量。右侧展示了基因的绝对表达量，基因长度和基因类型（即编码蛋白质或 lncRNA）。在热图的最左侧，是由 anno_block() 绘制的彩色矩形，用于区分根据 k-means 聚类识别出五个聚类。在“base mean”和“gene type”两个热图的顶部，汇总图（条形图和箱形图）显示了五个聚类的统计数据和分布。

library(ComplexHeatmap)library(circlize)expr = readRDS(system.file(package = "ComplexHeatmap", "extdata", "gene_expression.rds"))mat = as.matrix(expr[, grep("cell", colnames(expr))])base_mean = rowMeans(mat)mat_scaled = t(apply(mat, 1, scale))type = gsub("s\\d+_", "", colnames(mat))ha = HeatmapAnnotation(type = type, annotation_name_side = "left")ht_list = Heatmap(mat_scaled, name = "expression", row_km = 5,     col = colorRamp2(c(-2, 0, 2), c("green", "white", "red")),    top_annotation = ha,     show_column_names = FALSE, row_title = NULL, show_row_dend = FALSE) +Heatmap(base_mean, name = "base mean",     top_annotation = HeatmapAnnotation(summary = anno_summary(gp = gpar(fill = 2:6),         height = unit(2, "cm"))),    width = unit(15, "mm")) +rowAnnotation(length = anno_points(expr$length, pch = 16, size = unit(1, "mm"),     axis_param = list(at = c(0, 2e5, 4e5, 6e5),         labels = c("0kb", "200kb", "400kb", "600kb")),    width = unit(2, "cm"))) +Heatmap(expr$type, name = "gene type",     top_annotation = HeatmapAnnotation(summary = anno_summary(height = unit(2, "cm"))),    width = unit(15, "mm"))ht_list = rowAnnotation(block = anno_block(gp = gpar(fill = 2:6, col = NA)),     width = unit(2, "mm")) + ht_listdraw(ht_list)

甲基化，表达与其他基因组特征之间的相关性

# 下载 meth.rds：https://jokergoo.github.io/ComplexHeatmap-reference/data/meth.rdsres_list = readRDS("data/meth.rds")type = res_list$typemat_meth = res_list$mat_methmat_expr = res_list$mat_exprdirection = res_list$directioncor_pvalue = res_list$cor_pvaluegene_type = res_list$gene_typeanno_gene = res_list$anno_genedist = res_list$distanno_enhancer = res_list$anno_enhancer

这些变量的含义分别为：

•type：标记样本是肿瘤还是正常组织•mat_meth：矩阵，行对应于差异甲基化区域（DMR），矩阵中的值是每个样品中 DMR 的平均甲基化水平•mat_expr：矩阵，行对应于与 DMR 相关的基因（即与 DMR 最接近的基因），矩阵中的值是每个样品中对应基因的表达水平•direction：甲基化变化的方向•cor_pvalue：甲基化和基因表达之间的相关性 p 值•gene_type：基因的类型•anno_gene：基因的注释（基因间，基因内或 TSS）•dist：对应基因从 DMR 到 TSS 的距离•anno_enhancer：DMR中与增强子的重叠部分。

这些数据仅包括与其相关基因表达呈负相关的 DMR。

首先计算甲基化矩阵的列聚类，并据此调整表达矩阵中的列，使其顺序与甲基化矩阵保持一致。

column_tree = hclust(dist(t(mat_meth)))column_order = column_tree$order

library(RColorBrewer)meth_col_fun = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))direction_col = c("hyper" = "red", "hypo" = "blue")expr_col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))pvalue_col_fun = colorRamp2(c(0, 2, 4), c("white", "white", "red"))gene_type_col = structure(brewer.pal(length(unique(gene_type)), "Set3"),     names = unique(gene_type))anno_gene_col = structure(brewer.pal(length(unique(anno_gene)), "Set1"),     names = unique(anno_gene))dist_col_fun = colorRamp2(c(0, 10000), c("black", "white"))enhancer_col_fun = colorRamp2(c(0, 1), c("white", "orange"))

我们首先定义两列注释，然后进一步绘制更复杂的热图。

ht_opt(    legend_title_gp = gpar(fontsize = 8, fontface = "bold"),     legend_labels_gp = gpar(fontsize = 8),     heatmap_column_names_gp = gpar(fontsize = 8),    heatmap_column_title_gp = gpar(fontsize = 10),    heatmap_row_title_gp = gpar(fontsize = 8))ha = HeatmapAnnotation(type = type,     col = list(type = c("Tumor" = "pink", "Control" = "royalblue")),    annotation_name_side = "left")ha2 = HeatmapAnnotation(type = type,     col = list(type = c("Tumor" = "pink", "Control" = "royalblue")),     show_legend = FALSE)ht_list = Heatmap(mat_meth, name = "methylation", col = meth_col_fun,    column_order= column_order,    top_annotation = ha, column_title = "Methylation") +    Heatmap(direction, name = "direction", col = direction_col) +    Heatmap(mat_expr[, column_tree$order], name = "expression",         col = expr_col_fun,         column_order = column_order,         top_annotation = ha2, column_title = "Expression") +    Heatmap(cor_pvalue, name = "-log10(cor_p)", col = pvalue_col_fun) +    Heatmap(gene_type, name = "gene type", col = gene_type_col) +    Heatmap(anno_gene, name = "anno_gene", col = anno_gene_col) +    Heatmap(dist, name = "dist_tss", col = dist_col_fun) +    Heatmap(anno_enhancer, name = "anno_enhancer", col = enhancer_col_fun,         cluster_columns = FALSE, column_title = "Enhancer")draw(ht_list, row_km = 2, row_split = direction,    column_title = "Comprehensive correspondence between methylation, expression and other genomic features",     column_title_gp = gpar(fontsize = 12, fontface = "bold"),     merge_legends = TRUE, heatmap_legend_side = "bottom")