前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >聚类树的合并展示

聚类树的合并展示

作者头像
SYSU星空
发布2022-05-05 14:15:10
5050
发布2022-05-05 14:15:10
举报
文章被收录于专栏:微生态与微进化

往期回顾

层次聚类(hierarchical clustering)就是通过对数据集按照某种方法进行层次分解,直到满足某种条件为止,常用的方法有UPGMA、ward.D2等。聚类树是层次聚类最常用的可视化方法,我们可通过比较聚类来确定最佳分类,详见往期文章层次聚类与聚类树比较聚类

群落结构

通过层次聚类我们可以对微生物群落进行聚类并以聚类树的形式进行展示,但是要分析其生态学意义,我们需要结合更多的数据来对聚类簇进行解读。首先我们可以比较不同聚类簇中样品的群落结构的差异,分析不同微生物类群的变化规律,方法如下所示:

代码语言:javascript
复制
#读取物种和群落信息
data=read.table(file="otu_table.txt", header=TRUE, check.names=FALSE)
rownames(data)=data[, 1]
otu=as.matrix(data[,- 1])
community=read.table(file="phylum.txt",header=TRUE, check.names=FALSE)
rownames(community)=community[, 1]
com=as.matrix(community[, -1])
#计算物种相对丰度
library(vegan)
otu=decostand(otu, MARGIN=2, "total")*100
#物种丰度分组求均值
n=length(colnames(otu))
m=length(rownames(otu))
mean=rep(0, m*n/3)
otumean=matrix(mean, nrow=m, ncol=n/3)
rownames(otumean)=rownames(otu)
colnames(otumean)=LETTERS[1:22]
for (i in 1:m) {
  for (j in 1:(n/3)) {
    otumean[i, j]=mean(otu[i,(3*j-2):(3*j)])
  }
}
otumean=t(otumean)
#群落丰度数据求均值
p=ncol(com)
q=nrow(com)
commean=matrix(rep(0, p*q/3), nrow=q, ncol=p/3)
rownames(commean)=rownames(com)
colnames(commean)=LETTERS[1:22]
for (i in 1:q) {
  for (j in 1:(p/3)) {
    commean[i, j]=mean(com[i,(3*j-2):(3*j)])
  }
}
#计算距离矩阵
otu_dist=vegdist(otumean, method="bray", diag=TRUE, upper=TRUE, p=2)
#进行聚类分析
hclust=hclust(otu_dist, method="average")
#确定最佳聚类簇数目(这里省略,我们选聚类簇数目为3)
#聚类结果绘图
layout(matrix(c(1,2,3), 1, 3, byrow=TRUE), widths=c(1, 2, 1))
orhclust=reorder(hclust, otu_dist)
library(dendextend)
library(RColorBrewer)
tree=as.dendrogram(orhclust) %>% set("labels_cex", 2)
clusMember=cutree(tree, 3)
labelColors=brewer.pal(n=3, name="Set1")
colLab=function(n) {
  if (is.leaf(n)) {
    a=attributes(n)
    labCol=labelColors[clusMember[which(names(clusMember)==a$label)]]
    attr(n, "nodePar")=c(a$nodePar, lab.col=labCol)
  }
  n
}
#为聚类树设置颜色宽度等
tree=tree %>% set("labels_cex", 2) %>% set("branches_lwd", 2) %>% 
  set("branches_k_color", k=3) 
clusDendro=dendrapply(tree, colLab)
par(mar=c(5,2,3,2))
plot(clusDendro, type="rectangle", horiz=TRUE, xlab="Height")
#群落结构柱状图
#调整样品顺序与聚类树一致
commean=commean[,orhclust$order]
par(mar=c(5,1,3,1))
mycol=c(52,619,453,71,134,448,548,655,574,36,544,89,120,131,596,147,576,58,429,386,122,87,404,466,124,463,552,147,45,30,54,84,256,100,652,31,610,477,150,50,588,621,99,81,503,562,76,96,495)
mycol=colors()[mycol]
x=barplot(commean*100, xlab="Abundance (%)", ylab=NULL, col=mycol, border=NA, cex.axis=0.8, cex.names=0.8, xlim=c(0,100), yaxt="n", horiz=TRUE)
#添加样品名,用来检查样品名是否对应,最终作图可去除
#axis(2,pos=c(1.5,0), labels=colnames(commean), at=x, las=2, tick=FALSE, cex.axis=0.8 )
#添加图例
par(mar=c(1,0,1,1))
plot.new()
legend(x=0, y=0.98, legend=rownames(commean), ncol=1, fill=mycol, cex=1.5, border=NA, box.lwd=NA)

最终作图结果如下所示:

环境数据

前面的聚类均为对样品微生物群落数据进行分析,是一种非约束的聚类分析,我们可以根据聚类结果被动引入环境因子数据来进行比较,方法如下所示:

代码语言:javascript
复制
#读取物种和环境因子信息
data=read.table(file="otu_table.txt", header=TRUE, check.names=FALSE)
rownames(data)=data[, 1]
otu=as.matrix(data[,- 1])
environment=read.table(file="environment.txt",header=TRUE,check.names=FALSE)
rownames(environment)=environment[, 1]
env=data.frame(environment[, -1])
#计算物种相对丰度
library(vegan)
otu=decostand(otu, MARGIN=2, "total")*100
#物种丰度分组求均值
n=length(colnames(otu))
m=length(rownames(otu))
mean=rep(0, m*n/3)
otumean=matrix(mean, nrow=m, ncol=n/3)
rownames(otumean)=rownames(otu)
colnames(otumean)=LETTERS[1:22]
for (i in 1:m) {
  for (j in 1:(n/3)) {
    otumean[i, j]=mean(otu[i,(3*j-2):(3*j)])
  }
}
otumean=t(otumean)
#调整物种与环境因子样品ID顺序一致
otumean=otumean[rownames(env),]
#计算距离矩阵
otu_dist=vegdist(otumean, method="bray", diag=TRUE, upper=TRUE, p=2)
#进行聚类分析
hclust=hclust(otu_dist, method="average")
#确定最佳聚类簇数目(这里省略,我们选聚类簇数目为3)
#聚类结果绘图
layout(matrix(c(1,2,3,1,4,5), 2, 3, byrow=TRUE), widths=c(1, 1, 1))
orhclust=reorder(hclust, otu_dist)
library(dendextend)
library(RColorBrewer)
tree=as.dendrogram(orhclust) %>% set("labels_cex", 2)
clusMember=cutree(tree, 3)
labelColors=brewer.pal(n=3, name="Set1")
colLab=function(n) {
  if (is.leaf(n)) {
    a=attributes(n)
    labCol=labelColors[clusMember[which(names(clusMember)==a$label)]]
    attr(n, "nodePar")=c(a$nodePar, lab.col=labCol)
  }
  n
}
#为聚类树设置颜色宽度等
tree=tree %>% set("labels_cex", 2) %>% set("branches_lwd", 2) %>% 
  set("branches_k_color", k=3) 
clusDendro=dendrapply(tree, colLab)
par(mar=c(4,2,4,2))
plot(clusDendro, main ="UPGMA Tree", type="rectangle", horiz=TRUE)
#环境因子箱型图(环境因子筛选参照CCA等)
attach(env)
par(mar=c(6,5,8,2))
boxplot(Salinity~clusMember, pch="+", col=labelColors, range=0.8,  boxwex=0.5, ylab="Salinity")
boxplot(TN~clusMember, pch="+", col=labelColors, range=0.8,  boxwex=0.5, ylab="TN")
boxplot(TP~clusMember, pch="+", col=labelColors, range=0.8,  boxwex=0.5, ylab="TP")
boxplot(Temp~clusMember, pch="+", col=labelColors, range=0.8,  boxwex=0.5, ylab="Temperature")

作图结果如下所示:

示例数据下载链接:

https://pan.baidu.com/s/1HhLvdRY7rRrCIcGnCSLAYA

提取码:6rbl

END

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2020-11-25,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 微生态与微进化 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档