聚类树的合并展示

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