热图在单细胞数据分析中的应用

2020-04-08 10:55:34 浏览数 (1)

分享是一种态度

作者 | 周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

什么是热图？

热图是一个以颜色变化来显示数据的可视化矩阵，Toussaint Loua在1873年就曾使用过热图来绘制对巴黎各区的社会学统计。我们就拿这张简单朴素的热图来讲一下热图怎么看。

首先映入我们眼帘的是有的地方是黑的，有的地方是白的（颜色），每一块颜色都有对应的XY轴。言下之意，对象X的属性Y的值是用颜色表征的。颜色的聚集代表相应对象X的属性Y具有相似性（模式，pattern）。本质上它是表现一个数值矩阵，图上每一个小方格都是一个数值，按一条预设好的色彩变化尺（称为色键，Color Key），给每个数值分配颜色。

有时候我们还能看到对象X或者属性Y的聚类结果也绘制在热图的旁边，但是这就不属于热图的部分了，因为他已经不热了（热，就是有的地方冷，有的地方热）。

热图能说明哪些问题

表达量

广泛的应用就是用热图来可视化表达量。我们想象一下一个9个样本50个基因的表达谱，人类一眼看过去就是一堆数字，而表达量数值大小映射到颜色的深浅上，看起来就很清楚了。

很多时候，为了同一个基因在不同样本中的表达量有可比性，需要对表达量取对数，或取Z-score，把数据标准化到一个水平上。

相关性

计算两个矩阵的相关性，可以得到两两的相关性，这时，用热图的颜色来表示相关性可以看出哪些配对相关性较高。

在单细胞中的应用

表达量

seurat主题

这是一张典型的seurat做的热图，可以清楚地看出不同分群有着不同的表达模式。这里的每一个色块是一个细胞某基因的表达量。cluster可以看做是细胞的聚类，Y轴的基因，我们看到也是聚类了的（很可能是手动的，每一类基因作者都给出了注释）。所以这张热图的关键是什么？细胞和基因及其顺序。选择合适的细胞和基因（一般是每个群的差异高表达基因）后，为什么我做的图是一团黑？很可能是因为：

  1. 基因的顺序没排好。
  2. 没有标准化

人们经常需要根据差异分析的结果来探索基因列表的排序，如SC3的策略。差异基因的计算采用非参数Kruskal-Wallis检验。SC3提供了调整p值< 0.01的所有差异表达基因的列表，并绘制了p值最低的50个基因的基因表达谱。值得注意的是，聚类后的差异表达计算可能会在p值的分布中引入偏差，因此我们建议仅使用p值对基因进行排序。

SC3主题

这类图无疑反映了某geneList在某cluster的表达情况。如果巧了，这个geneList是某个细胞类型的marker基因，或者是某个功能的主要集合，热图有助于细胞群功能和类型的鉴定。热图很好地将对象（X，一般是我们的细胞）与它的属性（Y，一般是我们的基因）联系起来。

scanpy主题

在monocle2 中我们还看到一种热图将基因的表达情况与细胞发育轨迹结合到一起。可视化所有明显依赖于分支的基因的变化（如果愿意也可以自己定义geneList）。这张热图同时显示了两种命运的变化，它还要求选择分支点（branch_point ）。列是伪时间中的点，行是基因，伪时间的开始在热图的中间。当你从热图的中间读到右边的时候，你正在跟随一个伪时间谱系。当你读到左边时，另一个。这些基因是分层聚类的，因此您可以可视化具有类似的依赖于序列的基因模块。

monocle2主题

logFC

有时候为了突出展示群与群之间的差异基因,可以把差异分析的指标反映在热图上，如10X的loup软件那样：

Question 2: Which genes are displayed in the Loupe Cell Browser heatmap? How can I export the list of genes displayed?

Answer: In Loupe Cell Browser (version 2.0.0), the heatmap is a compact display of a subset of differentially expressed genes per cluster. Specifically, the gene list is the union of the top 120/N upregulated genes for each cluster ranked by log2 fold-change (N=total number of clusters).The gene names are on the plot when you export the heatmap. However, as of version 2.0.0, there is no 1-click function to export the associated information for the subset of heatmap genes.

cloup主题

相关性

提到相关性，我们很容易注意到WGCNA(weighted correlation network analysis，加权基因共表达网络分析）, 用于提取与性状或临床特征相关的基因模块，解析与表达量相关生物学过程。这是除了富集分析之外另一个寻找好的geneList的方法。这里的颜色不再是表达量的度量而是相似性的度量。

WGCNA主题

ComplexHeatmap在单细胞数据可视化中的应用

人们针对单细胞发展了相应的数据结构如seurat的S4类，monocle的CDS，SingleCellExperiment的sce,scanpy的anndata等，可见单细胞的故事远大于一张二维的表达谱。那么一张热图往往也不能完全的说明问题，于是我们希望能够灵活地操纵热图来讲更多的故事。于是，我们发现ComplexHeatmap这个R包真的是热图神器。

ComplexHeatmap 主题

代码语言：javascript复制
library(gdata)
library(ComplexHeatmap)
library(GetoptLong)
library(circlize)

getwd()
#?read.xls
expr <-readxl::read_excel("ComplexHeatmap//data//41587_2015_BFnbt3102_MOESM18_ESM.xlsx",sheet=2)
#View(expr)
#expr  <- as.matrix(expr)
expr = expr[!duplicated(expr[[1]]), ]
rownames(expr) = expr[[1]]
#expr = expr[-1]
colnames(expr<-as.data.frame(expr))
expr = expr[,-1]
expr = t(expr)
#expr = expr[-1,]
expr[1:3,1:3];dim(expr)
       Cell 1  Cell 2 Cell 3
Gnai3 3.23220 1.98320 2.2482
Cdc45 3.19810 1.17300 3.1705
Narf  0.29411 0.49389 1.6279
[1] 7073   81
代码语言：javascript复制
expr = expr[apply(expr, 1, function(x) sum(x > 0)/length(x) > 0.5), , drop = FALSE]
get_correlated_variable_genes = function(mat, n = nrow(mat), cor_cutoff = 0, n_cutoff = 0) {
  ind = order(apply(mat, 1, function(x) {
    q = quantile(x, c(0.1, 0.9))
    x = x[x < q[1] & x > q[2]]
    var(x)/mean(x)
  }), decreasing = TRUE)[1:n]
  mat2 = mat[ind, , drop = FALSE]
  dt = cor(t(mat2), method = "spearman")
  diag(dt) = 0
  dt[abs(dt) < cor_cutoff] = 0
  dt[dt < 0] = -1
  dt[dt > 0] = 1

  i = colSums(abs(dt)) > n_cutoff

  mat3 = mat2[i, ,drop = FALSE]
  return(mat3)
}

mat = get_correlated_variable_genes(expr, cor_cutoff = 0.5, n_cutoff = 20)
mat[1:4,1:4];dim(mat)
       Cell 1 Cell 2 Cell 3 Cell 4
Sdhd   3.7634 2.7851 3.0983 3.9585
Ccnd2  3.5098 4.2493 3.7505 4.5186
Dazap2 2.9348 2.0173 3.0520 4.7106
Lck    3.8987 1.8177 3.6390 4.2300
[1] 722  81
代码语言：javascript复制
mat2 = t(apply(mat, 1, function(x) {
  q10 = quantile(x, 0.1)
  q90 = quantile(x, 0.9)
  x[x < q10] = q10
  x[x > q90] = q90
  scale(x)
}))
colnames(mat2) = colnames(mat)
mat2[1:4,1:4]

           Cell 1     Cell 2      Cell 3   Cell 4
Sdhd    1.4222656 -0.1033475  0.38507327 1.458445
Ccnd2  -1.1291769  0.8112796 -0.49757737 1.517926
Dazap2 -0.2861832 -1.5546649 -0.05732535 1.585296
Lck     0.5384592 -1.5226640 -0.13645403 1.399448
代码语言：javascript复制
base_mean = rowMeans(mat)
cc = readRDS("ComplexHeatmap//data//mouse_cell_cycle_gene.rds")

head(cc)
ENSMUSG00000019256 ENSMUSG00000021866 ENSMUSG00000017716 ENSMUSG00000007815 ENSMUSG00000034218 ENSMUSG00000042750 
             "Ahr"           "Anxa11"            "Birc5"             "Rhoa"              "Atm"             "Bex2" 
ccl = rownames(mat) %in% cc
table(cc1)
cc_gene = rownames(mat)[ccl]
rp = readRDS("ComplexHeatmap//data//mouse_ribonucleoprotein.rds")

head(rp)
ENSMUSG00000029580 ENSMUSG00000028427 ENSMUSG00000067274 ENSMUSG00000000568 ENSMUSG00000022557 ENSMUSG00000017146 
            "Actb"             "Aqp7"            "Rplp0"           "Hnrnpd"             "Bop1"            "Brca1" 

rpl = rownames(mat) %in% rp
代码语言：javascript复制
ht_list = Heatmap(mat2, col = colorRamp2(c(-1.5, 0, 1.5), c("blue", "white", "red")), 
                  name = "scaled_expr", column_title = qq("relative expression for @{nrow(mat)} genes"),
                  show_column_names = FALSE, width = unit(8, "cm"),
                  heatmap_legend_param = list(title = "Scaled expr"))  
  Heatmap(base_mean, name = "base_expr", width = unit(5, "mm"),
          heatmap_legend_param = list(title = "Base expr"))  
  Heatmap(rpl   0, name = "ribonucleoprotein", col = c("0" = "white", "1" = "purple"), 
          show_heatmap_legend = FALSE, width = unit(5, "mm"))  
  Heatmap(ccl   0, name = "cell_cycle", col = c("0" = "white", "1" = "red"), 
          show_heatmap_legend = FALSE, width = unit(5, "mm"))  
  rowAnnotation(link = anno_mark(at = which(ccl & base_mean > quantile(base_mean, 0.25)), 
                                 labels = rownames(mat)[ccl & base_mean > quantile(base_mean, 0.25)], 
                                 labels_gp = gpar(fontsize = 10), padding = unit(1, "mm")))  
  Heatmap(cor(t(mat2)), name = "cor", 
          col = colorRamp2(c(-1, 0, 1), c("green", "white", "red")), 
          show_row_names = FALSE, show_column_names = FALSE, row_dend_side = "right", 
          show_column_dend = FALSE, column_title = "pairwise correlation between genes",
          heatmap_legend_param = list(title = "Correlation"))

ht_list = draw(ht_list, main_heatmap = "cor")

decorate_column_dend("scaled_expr", {
  tree = column_dend(ht_list)$scaled_expr
  ind = cutree(as.hclust(tree), k = 2)[order.dendrogram(tree)]

  first_index = function(l) which(l)[1]
  last_index = function(l) { x = which(l); x[length(x)] }
  x1 = c(first_index(ind == 1), first_index(ind == 2)) - 1
  x2 = c(last_index(ind == 1), last_index(ind == 2))
  grid.rect(x = x1/length(ind), width = (x2 - x1)/length(ind), just = "left",
            default.units = "npc", gp = gpar(fill = c("#FF000040", "#00FF0040"), col = NA))
})

dev.off()

数据可视化的过程就是一段探索意义的旅程，给每一种颜色、每一种形状、每一种聚集和离散找到一种生物学意义。这让我想起海子的《面朝大海，春暖花开》：

代码语言：javascript复制
给每一条河每一座山取一个温暖的名字
陌生人，我也为你祝福
愿你有一个灿烂的前程
愿你有情人终成眷属
愿你在尘世获得幸福
我只愿面朝大海，春暖花开

References

[1] ComplexHeatmap: https://links.jianshu.com/go?to=https://jokergoo.github.io/ComplexHeatmap-reference/book/ [2] R数据可视化3：热图: https://links.jianshu.com/go?to=https://db.cngb.org/science/ARTtoo54/ [3] 如何画热图: https://links.jianshu.com/go?to=https://zhuanlan.zhihu.com/p/21258579 [4] 10秒钟-完美掌握-热图（heatmap）绘制 - 所有人都可以！: https://www.jianshu.com/p/fb1170d03b97

数据安全

0 人点赞