library(permute)
library(lattice)
library(vegan)
otu <- read.delim('D:/Adata/First/otu.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))
distance <- vegdist(otu, method = 'bray')
pcoa <- cmdscale(distance, k = (nrow(otu) - 1), eig = TRUE)
write.table(as.matrix(distance), 'bray.csv', quote = F) #生成csv格式的文件,然后找到R语言文件列表下的bray.csv，对第一列进行分列，首行右移一个单元格（在第一个单元格插入，其余格右移）
beta = read.csv("bray.csv",header=T, row.names=1)
# k is dimension, 3 is recommended; eig is eigenvalues
pcoa = cmdscale(beta, k = 3, eig=T)
# get coordinate string, format to dataframme
points = as.data.frame(pcoa$points)
eig = pcoa$eig
design = read.table("D:/Adata/First/N48.txt", header=T, row.names=1, sep="\t") #导入OTU的分组信息表格
design$group1=design$group1
# rename group name
design$group2 = factor(design$group2, levels=c("RTB","MPB","NTB","RTR","MPR","NTR"))
points = cbind(points, design$group2)
colnames(points) = c("PC1", "PC2", "PC3", "group")
library(ggplot2)
p <- ggplot(points,aes(PC1, PC2,group=group,fill=group,shape = group))+
geom_point(data = points,aes(PC1, PC2,group=group,color=group,shape =group),size=8)+
geom_vline(xintercept = 0, color = 'black', size = 0.4,linetype="dashed") +
geom_hline(yintercept = 0, color = 'black', size = 0.4,linetype="dashed") +
labs(x=paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits=4), "%)", sep=""),
y=paste("PCoA 2 (", format(100 * eig[2] / sum(eig), digits=4), "%)", sep=""))+#x和y轴的坐标名称
scale_shape_manual(values = c(21,24,22, 21,24,22 )) + #可在这里修改点的形状
theme_bw()+theme(panel.grid=element_blank());p
p=p+
scale_color_manual(values =c("#00AAFF","#EE00FF","#AABB00","#0000DD","#FF0000","#88EE00"))+
scale_fill_manual(values =c("#00AAFF","#EE00FF","#AABB00","#0000DD","#FF0000","#88EE00"));p
p=p + stat_ellipse(aes(x=PC1, y=PC2, fill = group),
geom = "polygon", alpha = 0.2,level=0.68,linetype = 2,col="black");p
