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

> 分析过单细胞数据的小伙伴应该都使用过Seurat包，其中有个函数叫`DoHeatmap`，具体操作可以看： [单细胞转录组学习笔记-17-用Seurat包分析文章数据](https://www.jianshu.com/p/f6f54ce92e24)

## 前言

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

```r
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](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-12-04-071009.png)

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

## 如何用Pheatmap画这个结果？

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

### **第1步：得到表达矩阵**

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

```r
# 提取原始表达矩阵
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计算结果是这样，那么**矩阵的行就按基因名取：**

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-12-04-072247.png)

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

```r
# 原来的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
```

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

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

### **第3步：做一个列的注释**

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

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

```r
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
```

最后，就可以画图了：

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

![plot\_zoom\_png-2](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-12-04-073859.png)

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

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

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

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-12-04-073958.png)

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

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

例如：

```r
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](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-12-04-075333.png)