「热图」ComplexHeatmap展示单细胞聚类

实用Seurat自带的热图函数DoHeatmap绘制的热图,感觉有点不上档次,于是我尝试使用ComplexHeatmap这个R包来对结果进行展示。

个人觉得好的热图有三个要素

  • 聚类: 能够让别人一眼就看到模式
  • 注释: 附加注释能提供更多信息
  • 配色: 要符合直觉,比如说大部分都会认为红色是高表达,蓝色是低表达

在正式开始之前,我们需要先获取一下pbmc的数据,Seurat提供了R包SeuratData专门用于获取数据

devtools::install_github('satijalab/seurat-data')
library(SeuratData)
InstallData("pbmc3k")

加载数据并进行数据预处理,获取绘制热图所需的数据

library(SeuratData)
library(Seurat)
data("pbmc3k")
pbmc <- pbmc3k
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

pbmc.markers <- FindAllMarkers(pbmc, 
                               only.pos = TRUE, 
                               min.pct = 0.25, 
                               logfc.threshold = 0.25)

先感受下Seurat自带热图

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
Seurat-heatmap

下面则是介绍如何用R包ComplexHeatmap进行组图,虽然这个R包名带着Complex,但是并不是说这个R包很复杂,这个Complex应该翻译成复合,也就是说这个R包能在热图的基础上整合很多信息。

先安装并加载R包。

BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)

为了手动绘制一个热图,要从Seurat对象中提取所需要的表达量矩阵。我提取的是原始的count值,然后用log2(count + 1)的方式进行标准化

mat <- GetAssayData(pbmc, slot = "counts")
mat <- log2(mat + 1)

获取基因和细胞聚类信息

gene_features <- top10
cluster_info <- sort(pbmc$seurat_annotations)

对表达量矩阵进行排序和筛选

mat <- as.matrix(mat[top10$gene, names(cluster_info)])

Heatmap绘制热图。对于单细胞这种数据,一定要设置如下4个参数

  • cluster_rows= FALSE: 不作行聚类
  • cluster_columns= FALSE: 不作列聚类
  • show_column_names=FALSE: 不展示列名
  • show_row_names=FALSE: 不展示行名,基因数目不多时候可以考虑设置为TRUE
Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = TRUE)
Heatmap-1

从图中,我们可以发现以下几个问题:

  • 长宽比不合理,当然这和绘图函数无关,可以在保存时修改长宽比
  • 基因名重叠,考虑调整大小,或者不展示,或者只展示重要的基因
  • 颜色可以调整
  • 缺少聚类信息

这些问题,我们可以通过在ComplexHeatmap Complete Reference查找对应信息来解决。

配色方案

在热图中会涉及到两类配色,一种用来表示表达量的连续性变化,一种则是展示聚类。有一个神奇的R包就是用于处理配色,他的Github地址为https://github.com/caleblareau/BuenColors

devtools::install_github("caleblareau/BuenColors")
library("BuenColors")

它提供了一些列预设的颜色,比方说jdb_color_maps

      HSC       MPP      LMPP       CMP       CLP       MEP       GMP 
"#00441B" "#46A040" "#00AF99" "#FFC179" "#98D9E9" "#F6313E" "#FFA300" 
      pDC      mono     GMP-A     GMP-B     GMP-C       Ery       CD4 
"#C390D4" "#FF5A00" "#AFAFAF" "#7D7D7D" "#4B4B4B" "#8F1336" "#0081C9" 
      CD8        NK         B 
"#001588" "#490C65" "#BA7FD0"

这些颜色就能用于命名单细胞的类群,比如说我选择了前9个

col <- jdb_color_maps[1:9]
names(col) <- levels(cluster_info)

增加列聚类信息

Heatmaprow_splitcolumn_split参数可以通过设置分类变量对热图进行分隔。更多对热图进行拆分,可以参考Heatmap split

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info)
Heatmap-2

只用文字描述可能不够好看,最好是带有颜色的分块图,其中里面的颜色和t-SNE或UMAP聚类颜色一致,才能更好的展示信息。

为了增加聚类注释,我们需要用到HeatmapAnnotation函数,它对细胞的列进行注释,而rowAnnotation函数可以对行进行注释。这两个函数能够增加各种类型的注释,包括条形图,点图,折线图,箱线图,密度图等等,这些函数的特征是anno_xxx,例如anno_block就用来绘制区块图。

top_anno <- HeatmapAnnotation(
  cluster = anno_block(gp = gpar(fill = col), # 设置填充色
                       labels = levels(cluster_info), 
                       labels_gp = gpar(cex = 0.5, col = "white"))) # 设置字体

其中anno_block中的gp参数用于设置各类图形参数labels设置标签,labels_gp设置和标签相关的图形参数。可以用?gp来了解有哪些图形参数

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno, # 在热图上边增加注释
        column_title = NULL ) # 不需要列标题
Heatmap-3

突出重要基因

由于基因很多直接展示出来,根本看不清,我们可以强调几个标记基因。用到两个函数是rowAnnotationanno_mark

已知不同类群的标记基因如下

Cluster ID Markers Cell Type
0 IL7R, CCR7 Naive CD4+ T
1 IL7R, S100A4 Memory CD4+
2 CD14, LYZ CD14+ Mono
3 MS4A1 B
4 CD8A CD8+ T
5 FCGR3A, MS4A7 FCGR3A+ Mono
6 GNLY, NKG7 NK
7 FCER1A, CST3 DC
8 PPBP Platelet

我们需要给anno_mark提供基因所在行即可。

mark_gene <- c("IL7R","CCR7","IL7R","S100A4","CD14","LYZ","MS4A1","CD8A","FCGR3A","MS4A7","GNLY","NKG7","FCER1A", "CST3","PPBP")
gene_pos <- which(rownames(mat) %in% mark_gene)

row_anno <-  rowAnnotation(mark_gene = anno_mark(at = gene_pos, 
                                     labels = mark_gene))

接着绘制热图

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno,
        right_annotation = row_anno,
        column_title = NULL)
Heatmap-4

关于如何增加标记注释,参考mark-annotation

调增图例位置

目前的热图还有一个问题,也就是表示表达量范围的图例太占位置了,有两种解决方法

  • 参数设置show_heatmap_legend=FALSE直接删掉。
  • 利用heatmap_legend_param参数更改样式

我们根据legends这一节的内容进行一些调整

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno,
        right_annotation = row_anno,
        column_title = NULL,
        heatmap_legend_param = list(
          title = "log2(count+1)",
          title_position = "leftcenter-rot"
        ))
heatmap-5

因为ComplextHeatmap是基于Grid图形系统,因此可以先绘制热图,然后再用grid::draw绘制图例,从而实现将条形图的位置移动到图中的任意位置。

先获取绘制热图的对象

p <- Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno,
        right_annotation = row_anno,
        column_title = NULL,
        show_heatmap_legend = FALSE
        )

根据p@matrix_color_mapping获取图例的颜色的设置,然后用Legend构建图例

col_fun  <- circlize::colorRamp2(c(0, 1, 2 ,3, 4),
                                 c("#0000FFFF", "#9A70FBFF", "#D8C6F3FF", "#FFC8B9FF", "#FF7D5DFF"))
lgd <-  Legend(col_fun = col_fun, 
               title = "log2(count+1)", 
               title_gp = gpar(col="white", cex = 0.75),
               title_position = "leftcenter-rot",
               #direction = "horizontal"
               at = c(0, 1, 4), 
               labels = c("low", "median", "high"),
               labels_gp = gpar(col="white")
               )

绘制图形

grid.newpage() #新建画布
draw(p) # 绘制热图
draw(lgd, x = unit(0.05, "npc"), 
     y = unit(0.05, "npc"), 
     just = c("left", "bottom")) # 绘制图形
heatmap-6

ComplexHeatmap绘制热图非常强大的工具,大部分我想要的功能它都有,甚至我没有想到的它也有,这个教程只是展示其中一小部分功能而已,还有很多功能要慢慢探索。


版权声明:本博客所有文章除特别声明外,均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

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

推荐阅读更多精彩内容