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

2020-08-06 14:19:39 浏览数 (3)

随着时代的发展,单一研究转录组、蛋白代谢、甲基化等已经难以满足研究者越来越高的研究期望,大家更多地期望联合多种数据进行多组学联合分析。那么这时候,一种好的展示结果的方式无疑会为发表高分文章增光添彩。

本次,我们将展示一个甲基化与表达谱联合分析的热图。本着先学习再创造的态度,小编做了一下知识的搬运工,本篇所有代码均引用自:Zuguang Gu, Roland Eils and Matthias Schlesner, Complex heatmaps reveal patterns and correlations in multidimensional genomic data, Bioinformatics, 2016。

我们先上效果图:

其实代码并不长,关键在于如何准备作图所用数据以及对代码的理解上,所以下面将着重对这两点进行解释说明。这里的meth.rds从https://jokergoo.github.io/ComplexHeatmap-reference/data/meth.rds获取。

代码语言:javascript复制
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)

复杂的热图显示高度甲基化的DMR富含基因间和基因内区域,很少与增强子重叠。相反,低甲基化的DMR富含转录起始位点(TSS)和增强子。

知识点总结

1.ComplexHeatmap可实现单个热图的相加以实现数据之间的联合。

2.ht_global_opt()函数可实现整个热图的全局控制,但要注意使用结束后进行重置。

3.draw()函数在返回值是HeatmapList对象可以实现更多的控制。

0 人点赞