R语言- ComplexHeatmap 绘制复杂热图示例

热图(Heatmap)作为展示基因表达模式最直观的工具。它经常用在分子生物学文章里(尤其是microarray, RNA-seq相关论文)直观地呈现多样本多个基因的全局表达量变化,和呈现多样本或多基因表达量的聚类关系。绘制复杂热图最好用的是 ComplexHeatmap包,,它提供了灵活、高效、易于定制的方法来绘制多种类型的热图,并支持多种数据类型和数据格式,可以处理大型数据集,并在短时间内生成高质量的热图。

本文重点介绍如何用ComplexHeatmap::Heatmap的方法在热图上展示全样本的基因表达量情况的同时,在样本组水平上展示簇间的聚类关系(如果样本太多,则展示样本间的聚类关系将变得非常不直观)。

簇级聚类表达热图

1、示例数据准备

future::plan("multiprocess", workers = 6);options(future.globals.maxSize = 100000 * 1024^5) #设置任务多线程
##Data process >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
scRNA <- readRDS("test.RDS")  #加载Seurat对象数据
### 计算seurat_cluster间的 差异基因
markers <- FindAllMarkers(object = scRNA, only.pos = FALSE, min.pct = 0.25, logfc.threshold = 0.25)
### 挑选每个cluster top5 的基因画表达热图
select.features <- markers %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC) 
### 随机为 细胞样本分配一个模拟标签
scRNA@meta.data$sample <- sample(c("species.1","species.2","species.3","species.4"),size = ncol(scRNA),replace = T) 
### 提取Seurat绘制热图的矩阵数据
data <- Seurat::DoHeatmap(scRNA, features = select.features$gene, group.by = "seurat_clusters", group.bar = T, size = 4)$data 
### 提取细胞标识
cell.meta <- scRNA@meta.data %>% tibble::rownames_to_column(var="Cell") %>% 
  select(Cell,seurat_clusters,sample) %>% 
  arrange(seurat_clusters) %>% mutate(Identity = seurat_clusters)

### 长数据转换宽数据
counts <- data %>% select(!Identity) %>% 
  tidyr::pivot_wider(names_from =  Cell, values_from = Expression) %>% 
  data.frame() %>% tibble::column_to_rownames(var = "Feature") %>% select(cell.meta$Cell)
mat <- as.matrix(counts) #格式转换为矩阵

预处理后用于热图展示的数据集

2、安装本示例需要使用的软件包

BiocManager::install("ComplexHeatmap"); library(ComplexHeatmap) #画热图的包
BiocManager::install("dendextend"); library(dendextend) #用于树聚类的包
BiocManager::install("magick") #用于图像元光栅化处理的包
ht_opt$message = FALSE #忽略ComplexHeatmap包的提示信息

3、创建本示例需要使用的调色板

colormaps <- c(RColorBrewer::brewer.pal(name="Dark2", n = 8),RColorBrewer::brewer.pal(name="Paired", n = 12),RColorBrewer::brewer.pal(name="Set1", n = 9))
scales::show_col(colormaps)  
29色调色板

4、创建热图注释块

  • 簇级聚类树
# Dg_tree
### 创建表达矩阵的样本组间聚类树【计算组内样本均值进行建树】
dend1 = cluster_between_groups(mat, cell.meta$Identity)
### 你希望这些样本被聚类成几簇 【按树枝颜色区分】
dend1 = color_branches(dend1, k = 5)  
### 树样式调整
dend1 = dend1 %>% set("branches_lwd", 2) # 聚类树树枝线条 厚度
### dend1 = dend1 %>% raise.dendrogram (3) #聚类树底端线条厚度
### dend1 = dend1 %>% highlight_branches_col(viridis::viridis(100)) #聚类树颜色调整
### dend1 = dend1 %>% highlight_branches_col(rev(viridis::magma(1000)))  #聚类树颜色调整
  • 顶部列注释- 1 - 插入一个空注释行
# 注释1 empty
ha_top_1 <- HeatmapAnnotation(
  empty = anno_empty(border = FALSE,height = unit(0.1, "cm")), #添加空的注释块
  annotation_name_side = "left",which = "column" 
)
  • 顶部列注释- 2 - 插入块注释显示样本的簇编号
# 注释2 Group
### 获取树聚类后的矩阵样本的排列顺序
HM <- Heatmap(mat,cluster_columns = dend1)
HM = draw(HM)
### 根据树聚类的样本排列顺序 来排列细胞信息表cell.meta
group.data <- cell.meta[column_order(HM),] 
### 提取按树聚类排布的样本簇标签顺序
group_order_label <- unique(group.data$Identity,fromLast = F) %>% as.vector()

### 创建样本簇标识 色板
color.cl <- colormaps[seq(length(unique(cell.meta$Identity)))] 
### 创建簇标识色块注释对象
ha_top_2 <- HeatmapAnnotation( 
  Group = anno_block(gp = gpar(fill = color.cl,col = 0),
                     labels = group_order_label, #块注释标签
                     labels_gp = gpar(col = "white", fontface = "bold") , #注释文本样式
                     show_name = TRUE , #显示注释对象名
                     height = unit(0.5,"cm") # 注释对象的整体高度
                     # weight = unit(10,"cm") # 注释对象的整体宽度
                     ),
  annotation_name_side = "left",#注释对象名显示方向
  which = "column"
  )
### anno_block 块注释的图例构建
lgd_Group <- Legend(title = "Group", labels = group_order_label,legend_gp = gpar(fill = color.cl)) 
  • 顶部列注释- 3 - 添加样本的sample标签注释
# 注释3 Batch
ha_top_3 <- HeatmapAnnotation(Batch = cell.meta$sample, 
                              annotation_legend_param = list(Batch = list(title = "Batch",ncol=1)), #注释图例参数调整
                              annotation_name_side = "left",which = "column")
  • 右侧行注释- 4 - 统计每个基因的表达量(基于seurat标准化后的data矩阵)
### 统计每个基因的表达量
sum_Normexpr <- scRNA@assays$RNA@data[rownames(mat),] %>% Matrix::rowSums()
ha_rig = rowAnnotation(sum_Normexpr = anno_barplot(sum_Normexpr, bar_width = 1,gp = gpar(fill = "yellow",col="red"),
                                                   border=F, #行注释对象外侧边框
                                                   width = unit(2,"cm"), # 行注释的宽度
                                                   axis_param =(list(side = "top",gp=gpar(fontsize=5,col="red"))) # 坐标轴参数
                                                   ),
                       show_annotation_name = FALSE, #不显示注释对象标题
                       annotation_name_side = "bottom",# 注释标题旋转位置
                       annotation_name_gp= gpar(fontsize = 8), #注释标题大小
                       annotation_name_rot = 0 #注释标题旋转
                       )
### anno_block注释图例对象创建
lgd_sumExpr <- Legend(title = "sum_Norm_expr",at = "",legend_gp = gpar(fill = "yellow")) 

5、合并模块创建热图

Heatmap(mat,
        cluster_columns = dend1, #列方向添加 簇级 树聚类
        column_split = length(unique(cell.meta$Identity)), #热图列方向按簇拆分
        #热图主体
        column_dend_height = unit(2, "cm"), #树的高度
        clustering_method_columns = "spearson", #树的聚类方法
        column_title = "_OH_MY_Doheatmap_", #列方向大标题
        column_title_side = "bottom",
        column_title_gp = gpar(fontsize = 15, fontface = "bold"), #列方向大标题样式
        
        name = "Expr", #热图名称,表达量图例名
        cluster_rows = FALSE, #关闭行方向聚类
        show_column_names = FALSE, #关闭显示列名
        show_row_names = TRUE, #打开显示行名
        col = viridis::viridis(200), #表达量梯度颜色设置
        na_col = "black", #空值单元格的颜色
        row_title = "cluster_between_groups", #行方向大标题
        row_title_gp = grid::gpar(fontsize = 20,fontface="bold"), #行方向大标题样式
        row_names_side = "left", #行名显示方向
        row_names_gp = grid::gpar(fontsize = 6,fontface="bold"), #行名大小调整
        border = TRUE, #热图图像外边框显示
        
        # 表达量图例 样式设置
        heatmap_legend_param = list(
          title = "Exp",
          border = "red",
          direction = "vertical",
          title_position = "topleft"
          # legend_height = unit(12, "cm") # 热图表达量图例大小
        ),
        # 顶部注释
        top_annotation = c(ha_top_1,ha_top_2,ha_top_3), # 合并多个注释对象
        # 右注释
        right_annotation = ha_rig,
        
        ##图像 光栅化转换
        use_raster = TRUE, raster_quality = 5
) %>% draw(merge_legend = TRUE,padding = unit(c(1, 1, 2, 1), "cm"), # panding:图像编剧下-左-上-右
           annotation_legend_list = list(lgd_Group,lgd_sumExpr) # 添加 自己创建的 legend 对象
           )
decorate_column_dend("Expr", {grid.yaxis()}) # 树聚类 修饰

6、转换 Heatmapggplot 并保存出PDF

p <- grid.grabExpr(draw(ht.p))
ggsave(filename = "</FILE NAME/>.pdf",plot = p,width = *,height = *)

ok,以上就是首图所示热图样式的绘制的所有具体代码了,如有不懂欢迎留言一起讨论...

热图注释 - 简书
树聚类简介Introduction to dendextend (r-project.org)
viridis color maps 调色板简介
Chapter 5 Legends | ComplexHeatmap Complete Reference (jokergoo.github.io)
Chapter 14 More Examples | ComplexHeatmap Complete Reference (jokergoo.github.io)
Cluster groups in ComplexHeatmap - A Bioinformagician (jokergoo.github.io)

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 219,490评论 6 508
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,581评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 165,830评论 0 356
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,957评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,974评论 6 393
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,754评论 1 307
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,464评论 3 420
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,357评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,847评论 1 317
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,995评论 3 338
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,137评论 1 351
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,819评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,482评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,023评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,149评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,409评论 3 373
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,086评论 2 355

推荐阅读更多精彩内容