scRNA-如何改造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)
但是这个图在后期调整时会遇到很多障碍，因此最好用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)
当然，这个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的结果更像
