文献中的热图代码实现（热图标记感兴趣的基因,基础知识）
2021-06-15 本文已影响0人
单细胞空间交响乐
hello，大家好，今天我们来实现一下下面的这张热图
图片.png
我们最关注的就是图的右边只标记感兴趣的基因，如何实现呢？？？
准备作图文件
接下来，让我们来实现一下上面的热图。
为了绘制这种热图，首先要准备两类数据。
（1）基因表达矩阵，行是基因，列是样本，表达值可以是FPKM或者log转化后的表达值等都可以。
image
（2）基因名称列表，将待标识的重要基因名称以一排的形式放在一个列表中。
image
R包ComplexHeatmap的热图绘制
随后，将上述两个数据导入R中，绘制热图。
首先读取基因表达值矩阵，绘制一个常规的无基因名称的热图。
图片.png
随后，将重要的待展示名称的基因名称列表读入到R中，并添加在热图右侧显示出。
#读取待展示的基因名称，并添加到热图中
name <- read.delim('name.txt', header = FALSE, check.names = FALSE)
heat + rowAnnotation(link = anno_mark(at = which(rownames(mat) %in% name$V1),
labels = name$V1, labels_gp = gpar(fontsize = 10)))
图片.png
接下来我们来一个详细版的对比
图片
但是用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 <- pbmc3kpbmc[["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()
图片
作为探索性分析，这张图是可用的，但是可能无法直接用于最后文章的展示。
下面则是介绍如何用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)
图片
从图中，我们可以发现以下几个问题：
* 长宽比不合理，当然这和绘图函数无关，可以在保存时修改长宽比
* 基因名重叠，考虑调整大小，或者不展示，或者只展示重要的基因
* 颜色可以调整
* 缺少聚类信息
这些问题，我们可以通过在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)
增加列聚类信息
Heatmap的row_split和column_split参数可以通过设置分类变量对热图进行分隔。更多对热图进行拆分，可以参考Heatmap split
Heatmap(mat,
cluster_rows = FALSE,
cluster_columns = FALSE,
show_column_names = FALSE,
show_row_names = FALSE,
column_split = cluster_info)
图片
只用文字描述可能不够好看，最好是带有颜色的分块图，其中里面的颜色和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 ) # 不需要列标题
图片
突出重要基因
由于基因很多直接展示出来，根本看不清，我们可以强调几个标记基因。用到两个函数是rowAnnotation和anno_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)
图片
关于如何增加标记注释，参考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" ))
图片
因为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")) # 绘制图形
图片
ComplexHeatmap绘制热图非常强大的工具，大部分我想要的功能它都有，甚至我没有想到的它也有，这个教程只是展示其中一小部分功能而已，还有很多功能要慢慢探索。
基础知识，多多学习
