5.3.2 Seurat V3 | 如何改造Seurat包的DoHeatmap函数？

刘小泽写于19.12.4

分析过单细胞数据的小伙伴应该都使用过Seurat包，其中有个函数叫DoHeatmap，具体操作可以看： 单细胞转录组学习笔记-17-用Seurat包分析文章数据

前言

走完Seurat流程，会得到分群结果FindClusters()，并找到marker基因FindAllMarkers()，然后想要对每群的前10个marker基因进行热图可视化

rm(list = ls()) 
options(warn=-1) 
options(stringsAsFactors = F)

install.packages('Seurat')
library(Seurat)
library(stringr)   
library(dplyr)  

load('sce_out_for_heatmap_all.Rdata')
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
DoHeatmap(sce,top10$gene,size=3)

plot_zoom_png

但是这个图在后期调整时会遇到很多障碍，因此最好用pheatmap重新画一下

如何用Pheatmap画这个结果？

其实，画一个热图最重要的是表达矩阵和分组信息

第1步：得到表达矩阵

那么现在有了sce对象，从中提取表达矩阵也不难

# 提取原始表达矩阵
cts <- GetAssayData(sce, slot = "counts")
> cts[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
               A1_01 A1_10 A1_11 A1_12
00R-AC107638.2      .      2      .      2
0610005C13Rik       .      .      .      .
0610007P14Rik       1     49      3    328
0610009B22Rik       .     12      .     78
# 然后对这个矩阵取log
cts <- log10(cts + 1)

第2步：得到小的top10表达矩阵

当然不能使用整个表达矩阵进行处理，可以直接使用Seurat得到的top10结果，它帮我们得到了每个cluster的marker基因，而我们只需要取出这个小表达矩阵即可

Seurat得到的top10计算结果是这样，那么矩阵的行就按基因名取：

因为这个热图结果是按照cluster进行排序展示的，因此我们也要将小表达矩阵的列按cluster从小到大排序：

# 原来的cluster分组信息存储在 sce$seurat_clusters 中
> head(sce$seurat_clusters)
A1_01 A1_10 A1_11 A1_12 A1_13 A1_14 
     4      4      4      4      4      4 
Levels: 0 1 2 3 4 5 6 7 8
# 可见并不是从第0个cluster开始的

# 排序之后
new_cluster <- sort(sce$seurat_clusters)
> head(new_cluster)
A10_09 A10_16 A10_18 A10_33 A10_36 A10_42 
      0       0       0       0       0       0 
Levels: 0 1 2 3 4 5 6 7 8

👌有了行和列的规定，我们就能很轻松地提取出整个小的表达矩阵：

cts <- as.matrix(cts[top10$gene, names(new_cluster)])

第3步：做一个列的注释

因为要做出来DoHeatmap的顶部0-8 cluster的展示，需要使用pheatmap的一个参数：annotation_col

这个参数接收一个数据框作为输入。因为是对列进行注释，所以这个数据框的行是矩阵的列名，而它的列在这里对应的就是cluster分群信息

ac=data.frame(cluster=new_cluster)
rownames(ac)=colnames(mat)
> head(ac)
        cluster
A10_09       0
A10_16       0
A10_18       0
A10_33       0
A10_36       0
A10_42       0

最后，就可以画图了：

library(pheatmap)
pheatmap(cts,show_colnames =F,show_rownames = T,
         cluster_rows = F,
         cluster_cols = F,
         annotation_col=ac)

plot_zoom_png-2

当然，这个pheatmap有很多参数可以调整，看它的参数就知道：

比如想加上每个cluster的分隔，需要用到参数gaps_col = ，不过需要提供每个cluster最后一个细胞的序号

【当然这个是后话了，最重要的是知道如何进行数据的转换】

除了pheatmap，当然还能用其他的包

画热图的包有很多，其中一个比较常用的就是ComplexHeatmap

例如：

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

# 列标题的颜色框
color <- rainbow(9)
names(color) <- levels(new_cluster)
top_color <- HeatmapAnnotation(
  cluster = anno_block(gp = gpar(fill = color), 
                       labels = levels(new_cluster), 
                       labels_gp = gpar(cex = 0.5, col = "white"))) 

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = TRUE,
        column_split = new_cluster,
        heatmap_legend_param = list(
          title = "log10(count+1)",
          title_position = "leftcenter-rot"
        ),
        top_annotation = top_anno,
        column_title = NULL)

也可以实现列的注释，并且将每个cluster的列都进行分隔，和DoHeatmap的结果更像

plot_zoom_png-3