group.B.txt
| Sample | Group1 | Group2 |
|---|---|---|
| M3N2 | Control | Meth_7d |
| M3B2 | Control | Meth_7d |
| M3L4 | Control | Abst_14d |
| M3N1 | Control | Meth_0d |
1-B_bray_curtis.txt
| M1B1 | M1B2 | M1B3 | M1B4 | M1L1 | M1L2 | M1L3 | M1L4 | M1N1 | |
|---|---|---|---|---|---|---|---|---|---|
| M1B1 | 0 | 0.406891514 | 0.464991414 | 0.415630847 | 0.528281094 | 0.555857403 | 0.51824752 | 0.479445911 | 0.533112074 |
| M1B2 | 0.406891514 | 0 | 0.378623235 | 0.355250006 | 0.530049463 | 0.426138548 | 0.37125503 | 0.410992081 | 0.570619442 |
| M1B3 | 0.464991414 | 0.378623235 | 0 | 0.353481637 | 0.53799431 | 0.494976806 | 0.327378969 | 0.386427125 | 0.632358595 |
| M1B4 | 0.415630847 | 0.355250006 | 0.353481637 | 0 | 0.54327379 | 0.510866501 | 0.39643507 | 0.35325098 | 0.57880776 |
| M1L1 | 0.528281094 | 0.530049463 | 0.53799431 | 0.54327379 | 0 | 0.555716446 | 0.553281735 | 0.5228991 | 0.518093749 |
| M1L2 | 0.555857403 | 0.426138548 | 0.494976806 | 0.510866501 | 0.555716446 | 0 | 0.438799047 | 0.508508675 | 0.665893539 |
| M1L3 | 0.51824752 | 0.37125503 | 0.327378969 | 0.39643507 | 0.553281735 | 0.438799047 | 0 | 0.368858761 | 0.642802225 |
| M1L4 | 0.479445911 | 0.410992081 | 0.386427125 | 0.35325098 | 0.5228991 | 0.508508675 | 0.368858761 | 0 | 0.565058049 |
| M1N1 | 0.533112074 | 0.570619442 | 0.632358595 | 0.57880776 | 0.518093749 | 0.665893539 | 0.642802225 | 0.565058049 | 0 |
| M1N2 | 0.449844947 | 0.433673339 | 0.55933007 | 0.465196443 | 0.479958482 | 0.612099234 | 0.50192214 | 0.492516466 | 0.452817858 |
| M1N3 | 0.443463441 | 0.408582998 | 0.3720367 | 0.438606833 | 0.470065865 | 0.551782465 | 0.4882365 | 0.435185422 | 0.557190087 |
| M1R2 | 0.438594018 | 0.420602783 | 0.52598734 | 0.47296189 | 0.410453881 | 0.566941746 | 0.53967298 | 0.516286937 | 0.419206028 |
代碼
library(ggplot)
library(vegan)
design<-read.table("group.B.txt",header=T,row.names=1,sep="\t")
data<-read.table("1-B_bray_curtis.txt",sep="\t",header=T,check.names=F)
idx=rownames(design) %in% colnames(data)
sub_design=design[idx,]
data=data[rownames(sub_design),rownames(sub_design)]
pcoa=cmdscale(data,k=3,eig=T)
points=as.data.frame(pcoa\$points)
colnames(points)=c("x","y","z")
eig=pcoa$eig
points=cbind(points,sub_design[match(rownames(points),rownames(sub_design)),])
ggplot(points,aes(x=x,y=y,color=Group2)) +
geom_point(alpha=.7,size=2) +
stat_ellipse(level=0.95,show.legend=F) +
theme(panel.background=element_blank(),panel.border=element_rect(linetype="solid",fill=NA),
axis.text=element_text(size=10,color="black"),axis.title=element_text(size=12,face="bold",color="black")) +
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=""),
title="bray_curtis PCoA")
