1.问题
应该是会经常遇到这样的热图,就是说它的聚类情况无法与分组信息吻合。其实这个并不是错了,是因为用于聚类的行(选出的这部分基因)的表达模式在两组之间没有明显区别。
只要改变基因数量,或者换一组基因,聚类树就有可能变得和分组吻合。
假如你尝试换了很多次还是没有办法让它们变得吻合,那就取消聚类。直接使用热图的参数cluster_cols = F即可。
图就会变成:
取消了列聚类,那么热图列的顺序就会按照矩阵列的顺序排布了。那当样本数量多起来的时候,组内的数据规律就无法清晰的展示。
所以会需要组内聚类这样的操作咯。这个只能用complexheatmap来实现,我还顺便找到了注释条加分组标签的画法,变得比常规热图好看一丢。
2.代码如下
rm(list = ls())
library(tinyarray)
a = geo_download("GSE56649",destdir = tempdir())
a$exp = log2(a$exp+1)
find_anno(a$gpl)
## [1] "`library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)` and `ids <- AnnoProbe::idmap('GPL570')` are both avaliable"
Group = factor(rep(c("control","treat"),times = c(13,9)))
Group = factor(Group,levels = c("control","treat"))
ids <- AnnoProbe::idmap('GPL570')
exp = trans_array(a$exp,ids)
g = names(tail(sort(apply(exp, 1, sd)),100))
draw_heatmap(exp[g,],Group)
这就是常规的热图,默认行列都会聚类
draw_heatmap(exp[g,],Group,cluster_cols = F)
cluster_cols=F控制了列不聚类。
接下来是分组聚类热图:
library(ComplexHeatmap)
library(circlize)
col_fun = colorRamp2(c(-2, 0, 2), c("#2fa1dd", "white", "#f87669"))
top_annotation = HeatmapAnnotation(
cluster = anno_block(gp = gpar(fill = c("#2fa1dd", "#f87669")),
labels = c("control","treat"),
labels_gp = gpar(col = "white", fontsize = 12)))
m = Heatmap(t(scale(t(exp[g,]))),name = " ",
col = col_fun,
top_annotation = top_annotation,
column_split = Group,
show_heatmap_legend = T,
border = F,
show_column_names = F,
show_row_names = F,
column_title = NULL)
m