什么是熱圖狭莱？
熱圖是一個(gè)以顏色變化來(lái)顯示數(shù)據(jù)的可視化矩陣鹅髓，Toussaint Loua在1873年就曾使用過(guò)熱圖來(lái)繪制對(duì)巴黎各區(qū)的社會(huì)學(xué)統(tǒng)計(jì)筐带。我們就拿這張簡(jiǎn)單樸素的熱圖來(lái)講一下熱圖怎么看彪薛。
首先映入我們眼簾的是有的地方是黑的，有的地方是白的（顏色）坯汤，每一塊顏色都有對(duì)應(yīng)的XY軸虐唠。言下之意，對(duì)象X的屬性Y的值是用顏色表征的惰聂。顏色的聚集代表相應(yīng)對(duì)象X的屬性Y具有相似性（模式疆偿，pattern）。本質(zhì)上它是表現(xiàn)一個(gè)數(shù)值矩陣庶近，圖上每一個(gè)小方格都是一個(gè)數(shù)值翁脆，按一條預(yù)設(shè)好的色彩變化尺（稱為色鍵，Color Key）鼻种，給每個(gè)數(shù)值分配顏色。
有時(shí)候我們還能看到對(duì)象X或者屬性Y的聚類結(jié)果也繪制在熱圖的旁邊沙热，但是這就不屬于熱圖的部分了叉钥，因?yàn)樗呀?jīng)不熱了（熱，就是有的地方冷篙贸，有的地方熱）投队。
熱圖能說(shuō)明哪些問題
- 表達(dá)量
廣泛的應(yīng)用就是用熱圖來(lái)可視化表達(dá)量。我們想象一下一個(gè)9個(gè)樣本50個(gè)基因的表達(dá)譜爵川，人類一眼看過(guò)去就是一堆數(shù)字敷鸦，而表達(dá)量數(shù)值大小映射到顏色的深淺上，看起來(lái)就很清楚了寝贡。
很多時(shí)候扒披，為了同一個(gè)基因在不同樣本中的表達(dá)量有可比性，需要對(duì)表達(dá)量取對(duì)數(shù)圃泡，或取Z-score碟案，把數(shù)據(jù)標(biāo)準(zhǔn)化到一個(gè)水平上。
- 相關(guān)性
計(jì)算兩個(gè)矩陣的相關(guān)性颇蜡，可以得到兩兩的相關(guān)性价说，這時(shí)辆亏，用熱圖的顏色來(lái)表示相關(guān)性可以看出哪些配對(duì)相關(guān)性較高。
在單細(xì)胞中的應(yīng)用
- 表達(dá)量
這是一張典型的seurat做的熱圖鳖目，可以清楚地看出不同分群有著不同的表達(dá)模式扮叨。這里的每一個(gè)色塊是一個(gè)細(xì)胞某基因的表達(dá)量。cluster可以看做是細(xì)胞的聚類领迈，Y軸的基因彻磁，我們看到也是聚類了的（很可能是手動(dòng)的，每一類基因作者都給出了注釋）惦费。所以這張熱圖的關(guān)鍵是什么兵迅？ 細(xì)胞和基因及其順序。選擇合適的細(xì)胞和基因（一般是每個(gè)群的差異高表達(dá)基因）后薪贫，為什么我做的圖是一團(tuán)黑恍箭？很可能是因?yàn)椋?/p>
- 基因的順序沒排好。
- 沒有標(biāo)準(zhǔn)化
人們經(jīng)常需要根據(jù)差異分析的結(jié)果來(lái)探索基因列表的排序瞧省，如SC3的策略扯夭。差異基因的計(jì)算采用非參數(shù)Kruskal-Wallis檢驗(yàn)。SC3提供了調(diào)整p值< 0.01的所有差異表達(dá)基因的列表鞍匾，并繪制了p值最低的50個(gè)基因的基因表達(dá)譜交洗。值得注意的是，聚類后的差異表達(dá)計(jì)算可能會(huì)在p值的分布中引入偏差橡淑，因此我們建議僅使用p值對(duì)基因進(jìn)行排序构拳。
這類圖無(wú)疑反映了某geneList在某cluster的表達(dá)情況。如果巧了梁棠，這個(gè)geneList是某個(gè)細(xì)胞類型的marker基因置森，或者是某個(gè)功能的主要集合，熱圖有助于細(xì)胞群功能和類型的鑒定符糊。熱圖很好地將對(duì)象（X凫海，一般是我們的細(xì)胞）與它的屬性（Y，一般是我們的基因）聯(lián)系起來(lái)男娄。
在monocle2 中我們還看到一種熱圖將基因的表達(dá)情況與細(xì)胞發(fā)育軌跡結(jié)合到一起行贪。可視化所有明顯依賴于分支的基因的變化（如果愿意也可以自己定義geneList）模闲。這張熱圖同時(shí)顯示了兩種命運(yùn)的變化建瘫，它還要求選擇分支點(diǎn)（branch_point ）。列是偽時(shí)間中的點(diǎn)围橡，行是基因暖混，偽時(shí)間的開始在熱圖的中間。當(dāng)你從熱圖的中間讀到右邊的時(shí)候翁授，你正在跟隨一個(gè)偽時(shí)間譜系拣播。當(dāng)你讀到左邊時(shí)晾咪，另一個(gè)。這些基因是分層聚類的贮配，因此您可以可視化具有類似的依賴于序列的基因模塊.
- logFC
有時(shí)候?yàn)榱送怀稣故救号c群之間的差異基因,可以把差異分析的指標(biāo)反映在熱圖上谍倦，如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.
- 相關(guān)性
提到相關(guān)性，我們很容易注意到WGCNA(weighted correlation network analysis泪勒，加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析）, 用于提取與性狀或臨床特征相關(guān)的基因模塊昼蛀，解析與表達(dá)量相關(guān)生物學(xué)過(guò)程。這是除了富集分析之外另一個(gè)尋找好的geneList的方法圆存。這里的顏色不再是表達(dá)量的度量而是相似性的度量叼旋。
ComplexHeatmap在單細(xì)胞數(shù)據(jù)可視化中的應(yīng)用
人們針對(duì)單細(xì)胞發(fā)展了相應(yīng)的數(shù)據(jù)結(jié)構(gòu)如seurat的S4類，monocle的CDS沦辙，SingleCellExperiment的sce,scanpy的anndata等夫植，可見單細(xì)胞的故事遠(yuǎn)大于一張二維的表達(dá)譜。那么一張熱圖往往也不能完全的說(shuō)明問題油讯，于是我們希望能夠靈活地操縱熱圖來(lái)講更多的故事详民。于是，我們發(fā)現(xiàn)ComplexHeatmap這個(gè)R包真的是熱圖神器陌兑。
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
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
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
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
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()
數(shù)據(jù)可視化的過(guò)程就是一段探索意義的旅程沈跨，給每一種顏色、每一種形狀兔综、每一種聚集和離散找到一種生物學(xué)意義饿凛。這讓我想起海子的《面朝大海，春暖花開》：
給每一條河每一座山取一個(gè)溫暖的名字
陌生人软驰，我也為你祝福
愿你有一個(gè)燦爛的前程
愿你有情人終成眷屬
愿你在塵世獲得幸福
我只愿面朝大海笤喳，春暖花開
ComplexHeatmap
R數(shù)據(jù)可視化3：熱圖
如何畫熱圖
10秒鐘-完美掌握-熱圖（heatmap）繪制 - 所有人都可以！
