聚类热图怎么画-ARGs热图怎么画

前言:

仍然是三代测序数据的分析,宏基因组的文章中经常出现聚类热图和物种丰度图,用来直观地识别与某些疾病或者表型相关的菌群构成。

1.读取数据

一共有11个样本,每一个样本的测序reads都经过Nanopore官方的Epi2Me程序鉴定了物种,下表中第一列是被鉴定的菌种,第二列是该样本中每个物种产生的reads数目。

首先导入到R语言中,合并所有样本到一个数据框:

2.绘制热图

经过上一步,我们得到了列名为样本,行名为菌种的reads数据框,然后就可以绘制热图,进行聚类分析了:

绘制结果:

3.绘制物种丰度图

丰度图,其实就是堆积图,把每个样本的reads数目转换为百分数,然后作图就可以了:

绘制结果:

先通过K均值聚类算法,把这些基因聚成30类来画热图。

绘制热图,最常用的是R语言,当然也能用GraphPad Prism做。但如果数据量不是特别大的话,这里推荐一个在线工具ClustVis,但数据量大时还是容易卡顿。不同基因在不同样本中的表达量,可能会跨越好几个数量级,这样就会把色键分得很细,所以通常会对表达量取对数,或取Z-score,这样色键的范围一般就在个位数内变化。

另一个不太卡的在线工具叫Morpheus,博大研究所开发的。还有桌面软件HemI,华中科技大学开发的。

本文是要把昨天画的聚类热图给“掰弯”,

gene <- (file = "clipboard",header = T,sep = "\t", = F)  #读入数据

(gene) <- letters    #添加上基因名,如果直接有基因名,则省去这一步

gg <- hclust(dist(gene))   #行聚类

zz <- hclust(dist(t(gene)))  #列聚类

gene <- gene[,zz$order]   #聚类后排序

gene <- gene[gg$order,]

gene <- cbind(name = (gene),gene)

data <- melt(gene, = "name")

data$num <- rep(c(1:26),6)

data$x <- rep(c(1:6),each = 26)

tex <- (x = c(13:18),y = rep(0,6),variable = colnames(gene)[2:7])   #建立一个数据框,用于给环形热图添加标签

ggplot(data)+

geom_tree (data = gg,continuous = c("size"),size = 1)+    #生成聚类树

geom_tile(data = data,aes(x = x+12 ,y = num,alpha = 0.7,fill = value),color = "gray")+  #在热图的左边空出12列,方便画出空心的圆

scale_fill_gradient2(low = "green",high = "darkred",mid = "white",midpoint = 6,name = "FPKM",

guide = guide_colourbar(barheight = unit(10,"cm"), = element_text(size = rel(15)) ))+

scale_y_discrete(position = c("right"),limits = factor(c(1:28)),label = c(gene$name," "," "))+   # 在热图的底部空出两行(26+2),设置环形的开口

scale_x_discrete(label = NULL)+

scale_size_continuous(range = c(4,9),guide = NULL)+

scale_alpha(guide = NULL)+

labs(x = NULL,y = NULL)+

coord_polar(theta = "y",start = 0)+    #环形热图的本质是采用极坐标系

geom_text(data = tex,aes(x,y,label = variable))+

theme( = element_blank(),

= element_blank(),

.y = element_blank(),

.x = element_text(size = rel(2),hjust = 0.5),

.x = element_text(size = rel(2),hjust = 0.5),

.y = element_text(hjust = 0.5,size = rel(2)),

= element_blank())