社区首页 >专栏 >聚类树的合并展示

聚类树的合并展示

SYSU星空

发布于 2022-05-05 06:15:10

54900

代码可运行

文章被收录于专栏：微生态与微进化微生态与微进化

运行总次数：0

代码可运行

往期回顾

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

群落结构

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

#读取物种和群落信息
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)

最终作图结果如下所示：

环境数据

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

#读取物种和环境因子信息
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")