hello束昵，大家好，今天我們來實現(xiàn)一下下面的這張熱圖

圖片.png

我們最關(guān)注的就是圖的右邊只標記感興趣的基因醒颖，如何實現(xiàn)呢妻怎？？泞歉？

準備作圖文件

接下來逼侦，讓我們來實現(xiàn)一下上面的熱圖。

為了繪制這種熱圖腰耙，首先要準備兩類數(shù)據(jù)榛丢。

（1）基因表達矩陣，行是基因挺庞，列是樣本晰赞，表達值可以是FPKM或者log轉(zhuǎn)化后的表達值等都可以。

image

（2）基因名稱列表选侨，將待標識的重要基因名稱以一排的形式放在一個列表中掖鱼。

image

R包ComplexHeatmap的熱圖繪制

隨后，將上述兩個數(shù)據(jù)導入R中援制，繪制熱圖戏挡。

首先讀取基因表達值矩陣，繪制一個常規(guī)的無基因名稱的熱圖晨仑。

圖片.png

隨后褐墅，將重要的待展示名稱的基因名稱列表讀入到R中，并添加在熱圖右側(cè)顯示出洪己。

#讀取待展示的基因名稱妥凳，并添加到熱圖中
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自帶的熱圖函數(shù)DoHeatmap繪制的熱圖，其實是沒有這個效果答捕。于是我嘗試使用ComplexHeatmap這個R包來對結(jié)果進行展示逝钥。

個人覺得好的熱圖有三個要素

* 聚類: 能夠讓別人一眼就看到模式。如果本來就無法聚類拱镐，那圖也不好看晌缘。

* 注釋: 附加注釋能提供更多信息，比如說一些標記基因名

* 配色: 要符合直覺痢站，比如說大部分都會認為紅色是高表達磷箕，藍色是低表達

在正式開始之前，我們需要先獲取一下pbmc的數(shù)據(jù)阵难，Seurat提供了R包SeuratData專門用于獲取數(shù)據(jù)

devtools::install_github('satijalab/seurat-data')
library(SeuratData)
InstallData("pbmc3k")

加載數(shù)據(jù)并進行數(shù)據(jù)預處理岳枷，獲取繪制熱圖所需的數(shù)據(jù)

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包能在熱圖的基礎(chǔ)上整合很多信息。

先安裝并加載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繪制熱圖悟衩。對于單細胞這種數(shù)據(jù)，一定要設置如下4個參數(shù)

* cluster_rows= FALSE: 不作行聚類

* cluster_columns= FALSE: 不作列聚類

* show_column_names=FALSE: 不展示列名

* show_row_names=FALSE: 不展示行名栓拜，基因數(shù)目不多時候可以考慮設置為TRUE

Heatmap(mat,        
cluster_rows = FALSE,        
cluster_columns = FALSE,        
show_column_names = FALSE,        
show_row_names = TRUE)

圖片

從圖中座泳，我們可以發(fā)現(xiàn)以下幾個問題：

* 長寬比不合理，當然這和繪圖函數(shù)無關(guān)幕与，可以在保存時修改長寬比

* 基因名重疊挑势，考慮調(diào)整大小，或者不展示啦鸣，或者只展示重要的基因

* 顏色可以調(diào)整

* 缺少聚類信息

這些問題潮饱，我們可以通過在ComplexHeatmap Complete Reference查找對應信息來解決。

配色方案

在熱圖中會涉及到兩類配色赏陵，一種用來表示表達量的連續(xù)性變化饼齿，一種則是展示聚類饲漾。有一個神奇的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參數(shù)可以通過設置分類變量對熱圖進行分隔。更多對熱圖進行拆分僚楞，可以參考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函數(shù)膜赃，它對細胞的列進行注釋挺邀，而rowAnnotation函數(shù)可以對行進行注釋。這兩個函數(shù)能夠增加各種類型的注釋，包括條形圖端铛，點圖泣矛，折線圖，箱線圖禾蚕，密度圖等等您朽，這些函數(shù)的特征是anno_xxx，例如anno_block就用來繪制區(qū)塊圖换淆。

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參數(shù)用于設置各類圖形參數(shù)哗总，labels設置標簽，labels_gp設置和標簽相關(guān)的圖形參數(shù)产舞』臧拢可以用?gp來了解有哪些圖形參數(shù)。

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 ) # 不需要列標題

圖片

突出重要基因

由于基因很多直接展示出來易猫，根本看不清耻煤，我們可以強調(diào)幾個標記基因。用到兩個函數(shù)是rowAnnotation和anno_mark

已知不同類群的標記基因如下

我們需要給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)

圖片

關(guān)于如何增加標記注釋哈蝇，參考mark-annotation

調(diào)增圖例位置

目前的熱圖還有一個問題，也就是表示表達量范圍的圖例太占位置了攘已，有兩種解決方法

* 參數(shù)設置show_heatmap_legend=FALSE直接刪掉炮赦。

* 利用heatmap_legend_param參數(shù)更改樣式

我們根據(jù)legends這一節(jié)的內(nèi)容進行一些調(diào)整

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圖形系統(tǒng)，因此可以先繪制熱圖样勃，然后再用grid::draw繪制圖例吠勘，從而實現(xiàn)將條形圖的位置移動到圖中的任意位置。

先獲取繪制熱圖的對象

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        )

根據(jù)p@matrix_color_mapping獲取圖例的顏色的設置峡眶，然后用Legend構(gòu)建圖例

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繪制熱圖非常強大的工具剧防，大部分我想要的功能它都有，甚至我沒有想到的它也有辫樱，這個教程只是展示其中一小部分功能而已峭拘，還有很多功能要慢慢探索。

基礎(chǔ)知識狮暑，多多學習