circlize 绘制复杂热图
测试开头
点击上方关注老俊俊！
测试结尾
今日大寒
一个努力中的公众号
长的好看的人都关注了
1引言
跟着 circlize 官方文档绘制一个复杂的热图。
2数据介绍
加载测试数据
library(circlize)
# 加载测试数据
source("https://gist.githubusercontent.com/jokergoo/0ea5639ee25a7edae3871ed8252924a1/raw/57ca9426c2ed0cebcffd79db27a024033e5b8d52/random_matrices.R")
环境变量会多出这么些数据:
数据查看
差异甲基化数据:
mat_meth[1:5,1:5]
sample1 sample2 sample3 sample4 sample5
[1,] 0.09508993 0.10253955 0.1242478 0.2108422 0.1147976
[2,] 0.17894877 0.20436340 0.3663206 0.2851402 0.1415034
[3,] 0.19883432 0.09013638 0.1634693 0.1049659 0.1783066
[4,] 0.33309078 0.22795393 0.1966336 0.2611358 0.4468088
[5,] 0.37314951 0.37459873 0.4272505 0.3935484 0.4477683
差异甲基化数据对应的基因表达量数据:
mat_expr[1:5,1:5]
sample1 sample2 sample3 sample4 sample5
[1,] 1.2638068 0.41371194 0.1794433 -1.39548351 -0.4535142
[2,] -0.6937921 -0.01734942 1.2455341 -0.07400419 2.1798801
[3,] 0.5577426 0.62261785 1.5788161 0.18985637 -0.7171888
[4,] 0.4348747 1.90738289 0.1454953 1.59168129 -0.3831575
[5,] -0.6493806 -0.36898746 -0.9927603 -1.60419279 -0.9087783
甲基化数据的分类:
direction[1:5]
[1] "hypo" "hypo" "hypo" "hypo" "hyper"
table(direction)
direction
hyper hypo
595 405
甲基化和表达量-log10 的相关性数据:
cor_pvalue[1:5]
[1] 0.6560338 0.5782283 0.3494940 0.6809355 0.1278238
基因分类:
gene_type[1:5]
[1] "protein_coding" "psedo-gene" "psedo-gene" "lincRNA" "others"
table(gene_type)
gene_type
lincRNA microRNA others protein_coding psedo-gene
112 99 117 582 90
基因注释:
anno_gene[1:5]
[1] "intergenic" "intergenic" "intergenic" "intergenic" "intergenic"
table(anno_gene)
anno_gene
intergenic intragenic TSS
419 483 98
差异甲基化区域里 TSS 位点的距离:
dist[1:5]
[1] 388737 97280 36426 420011 355666
差异甲基化区域与 enhancer 重叠区域:
anno_enhancer[1:5,1:3]
enhancer_1 enhancer_2 enhancer_3
1 0.07196473 0.0000000 0.2879309
2 0.75837107 0.3554948 0.0000000
3 0.00000000 0.8477110 0.0000000
4 0.00000000 0.0000000 0.0000000
5 0.57110298 0.0000000 0.9173112
3热图
聚类提取类群编号:
set.seed(123)
km = kmeans(mat_meth, centers = 5)$cluster
绘制热图:
# 甲基化热图
col_meth = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))
circos.heatmap(mat_meth, split = km, col = col_meth, track.height = 0.12)
# 分类
col_direction = c("hyper" = "red", "hypo" = "blue")
circos.heatmap(direction, col = col_direction, track.height = 0.01)
# 基因表达热图
col_expr = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
circos.heatmap(mat_expr, col = col_expr, track.height = 0.12)
# 相关性p值
col_pvalue = colorRamp2(c(0, 2, 4), c("white", "white", "red"))
circos.heatmap(cor_pvalue, col = col_pvalue, track.height = 0.01)
# 基因分类
library(RColorBrewer)
col_gene_type = structure(brewer.pal(length(unique(gene_type)), "Set3"), names = unique(gene_type))
circos.heatmap(gene_type, col = col_gene_type, track.height = 0.01)
# 基因注释
col_anno_gene = structure(brewer.pal(length(unique(anno_gene)), "Set1"), names = unique(anno_gene))
circos.heatmap(anno_gene, col = col_anno_gene, track.height = 0.01)
# TSS距离
col_dist = colorRamp2(c(0, 10000), c("black", "white"))
circos.heatmap(dist, col = col_dist, track.height = 0.01)
# enhancer重叠数据
col_enhancer = colorRamp2(c(0, 1), c("white", "orange"))
circos.heatmap(anno_enhancer, col = col_enhancer, track.height = 0.03)
4添加连接线
# 制作连接数据
df_link = data.frame(
from_index = sample(nrow(mat_meth), 20),
to_index = sample(nrow(mat_meth), 20)
)
# 查看
head(df_link)
from_index to_index
1 938 519
2 818 426
3 118 649
4 299 766
5 229 211
6 244 932
# 绘图
for(i in seq_len(nrow(df_link))) {
circos.heatmap.link(df_link$from_index[i],
df_link$to_index[i],
col = rand_color(1))
}
以上绘图可以打包为一个函数:
# 打包为函数
circlize_plot = function() {
circos.heatmap(mat_meth, split = km, col = col_meth, track.height = 0.12)
circos.heatmap(direction, col = col_direction, track.height = 0.01)
circos.heatmap(mat_expr, col = col_expr, track.height = 0.12)
circos.heatmap(cor_pvalue, col = col_pvalue, track.height = 0.01)
circos.heatmap(gene_type, col = col_gene_type, track.height = 0.01)
circos.heatmap(anno_gene, col = col_anno_gene, track.height = 0.01)
circos.heatmap(dist, col = col_dist, track.height = 0.01)
circos.heatmap(anno_enhancer, col = col_enhancer, track.height = 0.03)
for(i in seq_len(nrow(df_link))) {
circos.heatmap.link(df_link$from_index[i],
df_link$to_index[i],
col = rand_color(1))
}
circos.clear()
}
5添加图例
# 绘制图例
library(ComplexHeatmap)
lgd_meth = Legend(title = "Methylation", col_fun = col_meth)
lgd_direction = Legend(title = "Direction", at = names(col_direction),
legend_gp = gpar(fill = col_direction))
lgd_expr = Legend(title = "Expression", col_fun = col_expr)
lgd_pvalue = Legend(title = "P-value", col_fun = col_pvalue, at = c(0, 2, 4),
labels = c(1, 0.01, 0.0001))
lgd_gene_type = Legend(title = "Gene type", at = names(col_gene_type),
legend_gp = gpar(fill = col_gene_type))
lgd_anno_gene = Legend(title = "Gene anno", at = names(col_anno_gene),
legend_gp = gpar(fill = col_anno_gene))
lgd_dist = Legend(title = "Dist to TSS", col_fun = col_dist,
at = c(0, 5000, 10000), labels = c("0kb", "5kb", "10kb"))
lgd_enhancer = Legend(title = "Enhancer overlap", col_fun = col_enhancer,
at = c(0, 0.25, 0.5, 0.75, 1), labels = c("0%", "25%", "50%", "75%", "100%"))
# 添加到图层
library(gridBase)
plot.new()
circle_size = unit(1, "snpc") # snpc unit gives you a square region
pushViewport(viewport(x = 0, y = 0.5, width = circle_size, height = circle_size,
just = c("left", "center")))
par(omi = gridOMI(), new = TRUE)
circlize_plot()
upViewport()
h = dev.size()[2]
lgd_list = packLegend(lgd_meth, lgd_direction, lgd_expr, lgd_pvalue, lgd_gene_type,
lgd_anno_gene, lgd_dist, lgd_enhancer, max_height = unit(0.9*h, "inch"))
draw(lgd_list, x = circle_size, just = "left")
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,数据代码已上传至QQ群,欢迎加入下载。
群二维码：
老俊俊微信：
知识星球：
所以今天你学习了吗？
欢迎小伙伴留言评论！
今天的分享就到这里了，敬请期待下一篇！
最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！
如果觉得对您帮助很大，赏杯快乐水喝喝吧！
往期回顾
◀组会文献分享 — miR455-3p 影响 HSF1 基因的 m6A 修饰参与结直肠癌的发生发展
◀基于 featureCounts 原理提取基因非冗余外显子长度
◀python 学习之 featureCounts 软件的基因长度是怎么算的?
◀...
请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！
