R绘制甲基化和表达谱联合分析热图

library(ComplexHeatmap)
library(circlize)
res_list = readRDS("meth.rds")
summary(res_list)
type = res_list$type
mat_meth = res_list$mat_meth
mat_expr = res_list$mat_expr
direction = res_list$direction
#甲基化与甲基化关联基因相关性p值
cor_pvalue = res_list$cor_pvalue
#基因类型(如蛋白编码基因或lincRNA等)
gene_type = res_list$gene_type
#DMR注释到基因的功能区间(如intergenic/intragenic或者TSS
anno_gene = res_list$anno_gene
#DMR到关联基因TSS的距离
dist = res_list$dist
#与增强子重叠的DMR的部分
anno_enhancer = res_list$anno_enhancer

##首先计算甲基化矩阵的列聚类，以便可以将表达矩阵中的列调整为具有与甲基化矩阵中相同的列顺序。
column_tree = hclust(dist(t(mat_meth)))
column_order = column_tree$order

#颜色定义参见上一篇(如何让你的图变得高大上之ComplexHeatmap())
library(RColorBrewer)
#定义甲基化表达水平颜色，从0/blue-0.5/white-1/red渐变
meth_col_fun = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))
#定义甲基化变化方向对应颜色
direction_col = c("hyper" = "red", "hypo" = "blue")
#定义表达水平颜色
expr_col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
#定义相关性p值颜色
pvalue_col_fun = colorRamp2(c(0, 2, 4), c("white", "white", "red"))
#定义基因类型颜色
gene_type_col = structure(brewer.pal(length(unique(gene_type)), "Set3"), 
                          names = unique(gene_type))
#定义注释model颜色
anno_gene_col = structure(brewer.pal(length(unique(anno_gene)), "Set1"), 
                          names = unique(anno_gene))
#定义距离颜色
dist_col_fun = colorRamp2(c(0, 10000), c("black", "white"))
#定义增强子相关颜色
enhancer_col_fun = colorRamp2(c(0, 1), c("white", "orange"))

##我们首先定义两个列注释，然后制作复杂的热图。
#ht_global_opt()是一个可选函数，它会全局控制一些参数。我们可以通过此全局函数同时为所有热图/注释设置一些参数。需要注意的是，一定将它放在热图代码(也就是Heatmap())之前，并在绘制热图后重置所有选项值以消除对下一个热图的影响。
#可以通过?ComplexHeatmap::ht_global_opt查看此函数的帮助
names(ht_global_opt())        #可查看该函数可定义的参数

ht_global_opt(
  heatmap_legend_title_gp = gpar(fontsize = 8, fontface = "bold"), 
  heatmap_legend_labels_gp = gpar(fontsize = 8), 
  heatmap_column_names_gp = gpar(fontsize = 8),
  heatmap_column_title_gp = gpar(fontsize = 10),
  heatmap_row_title_gp = gpar(fontsize = 8),ADD=TRUE
)
#利用HeatmapAnnotation()对行或列注释。HeatmapAnnotation()函数会生成一个注释用的list对象。该函数的主要格式是HeatmapAnnotation(df/数据框, name/注释名称, col/注释颜色列表, show_legend/是否显示数据框中每一列的图例)
#样本类型注释，Tumor样本为pink，Control样本为royalbule，名称在左边
ha = HeatmapAnnotation(type = type, 
                       col = list(type = c("Tumor" = "pink", "Control" = "royalblue")),
                       annotation_name_side = "left")
#不显示图例名称
ha2 = HeatmapAnnotation(type = type, 
                        col = list(type = c("Tumor" = "pink", "Control" = "royalblue")), 
                        show_legend = FALSE)

#Heatmap()实际上是单一热图的类构造函数。如果需要组合超过一个热图，用户可以通过+操作符添加热图。默认情况下，将两个热图通过+连接后，第二个热图的行聚类树会去掉，行的顺序会与是第一个热图的顺序保持一致。
ht_list = Heatmap(mat_meth, name = "methylation", col = meth_col_fun,
                  column_order= column_order,
                  top_annotation = ha, column_title = "Methylation")+
  Heatmap(direction, name = "direction", col = direction_col) +
  Heatmap(mat_expr[, column_tree$order], name = "expression", 
          col = expr_col_fun, 
          column_order = column_order, 
          top_annotation = ha2, column_title = "Expression") +
  Heatmap(cor_pvalue, name = "-log10(cor_p)", col = pvalue_col_fun)+ 
  Heatmap(gene_type, name = "gene type", col = gene_type_col) +
  Heatmap(anno_gene, name = "anno_gene", col = anno_gene_col) +
  Heatmap(dist, name = "dist_tss", col = dist_col_fun) +
  Heatmap(anno_enhancer, name = "anno_enhancer", col = enhancer_col_fun, 
          cluster_columns = FALSE, column_title = "Enhancer")
draw(ht_list, km = 2, split = direction,
     column_title = "Comprehensive correspondence between methylation,
expression and other genomic features", 
     column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
     merge_legends = TRUE, heatmap_legend_side = "bottom")
#重置全局参数消除影响  
ht_global_opt(RESET = TRUE)