Heatmap与层次聚类

Tang Ming gave an implementation of improved heatmap using ComplexHeatmap.

Do not show too many cells on a single screen.

  • When you have too many cells (> 10,000), the use_raster option really helps. Also consider downsample the Seurat object to a smaller number of cells for plotting the heatmap. Your screen resolution is not as high as 300,000 pixels if you have 300,000 cells (columns).

https://divingintogeneticsandgenomics.rbind.io/post/enhancement-of-scrnaseq-heatmap-using-complexheatmap/

For hierarchical clustering, consult this hclust tutorial

https://uc-r.github.io/hc_clustering#algorithms

 

基于Seurat数据创建clustered heatmap

0. 导入需要的包

suppressPackageStartupMessages({
    library(ComplexHeatmap)
    library(circlize)
})

1. 准备数据

# 1.1 收集表达矩阵
mat  = GetAssayData(t.combined, slot = "data", assay = "RNA")

# 1.2 把细胞按某种想要的顺序排好
ordered.cells <- rownames(t.combined@meta.data)

# 1.3 选择一列感兴趣的基因用于绘图
features.sel <-  unique(c("CD3E","CD3D","CD4","CD8A","CD8B"))

# 1.4 重新组织原始矩阵按行列
mat = as.matrix( mat[features.sel, ordered.cells] )

2. 对数据作层次聚类(可选)

# 2.1 层次聚类
## 注意,mat是从Seurat来的feature-by-sample矩阵,应该转置成sample-by-feature矩阵;
## 此处选择一些基因做聚类,然后用另一些基因画zheng
distobj   <- dist(t(mat[c("CD3D","CD3E","CD4","CD8A","CD8B"),]), method = "euclidean")
hclustobj <- hclust(distobj, method = "ward.D2" )

# 2.2 切隔树形成聚类
sub_grp <- cutree(hclustobj, k = 10) #参数k直接指定类别数,也可用参数h指定切割深度
# 打印各类里的细胞数看看
table(sub_grp) 
# 打印聚类树看看
options(repr.plot.width=15, repr.plot.height=5)
plot(hclustobj, cex = 0.1, label=F) #画🌲,label=F隐藏样本名,不然过密不美观
rect.hclust(hclustobj, k = 10) #画方框
# 按聚类树顺序画heatmap看看
library(pheatmap)
options(repr.plot.width=15, repr.plot.height=10)
pheatmap(mat[,hclustobj$order], cluster_rows = F, cluster_cols = F, 
         border_color=NA,
         show_colnames = F, use_raster=TRUE)

# 看看要是还行,就把层次聚类的结果标签弄进metadata去
t.combined@meta.data["subgroup"]<-as.character(sub_grp)

3. 创建meta数据框,并构建annotations

meta = t.combined@meta.data[ordered.cells,]

# configure colors for cell type annotations
library(paletteer)  
cl_levels = unique(meta$subgroup)
blockcol = paletteer_d("ggsci::default_igv")[1:length(cl_levels)] %>% as.vector
names(blockcol) <-cl_levels


colann <- HeatmapAnnotation(
    cluster = meta$cluster,
    organ   = meta$organ,
    subgroup= meta$subgroup,
    col = list( subgroup = blockcol ),
    annotation_legend_param=list(
        cluster = list(nrow=5),
        organ = list(nrow=3),
        subgroup=list(nrow=1)
    )
)

4. 绘制Heatmap

hm<-Heatmap(mat, name = "Normalized expression", 
        cluster_rows = T, 
        cluster_columns = hclustobj, show_column_names=FALSE,
        column_dend_height = unit(4, "cm"),
        
        #column_split=meta$organ, cluster_column_slices=T,
        col= colorRamp2(c(0,1.5,3), c("#486E9E", "white", "#D84B59")),
        column_title_rot=90, column_gap=unit(2, "mm"),
        top_annotation=colann, heatmap_legend_param = list(direction = "horizontal"),
        use_raster = TRUE)

options(warn=-1)
options(repr.plot.width=20, repr.plot.height=20)
#cairo_pdf("XXXX.hm.pdf",width=15, height = 18)
draw(hm,
     padding = unit(c(10, 10, 30, 3), "mm"), #下左上右
     merge_legend = TRUE,
     heatmap_legend_side = "bottom", 
     annotation_legend_side = "bottom")
#dev.off()

Reply

您的电子邮箱地址不会被公开。 必填项已用 * 标注