跟着Nature学数据分析：R语言vegan包做主坐标分析（PCoA）并使用ggplot2可视化结果

inMBss<-read.table("20220420/nature/inMBss.tsv",
                   row.names = 1,
                   header=T)
dim(inMBss)

library(vegan)
distMatrix <- vegdist(inMBsss,method = "bray")
pCoa <- cmdscale(distMatrix, eig = T,k = 2 )

varExp <- (eigenvals(pCoa)/sum(eigenvals(pCoa)))[1:2]
xVar <- as.numeric(varExp[1]*100)
yVar <- as.numeric(varExp[2]*100)

pCoaVecs <- as.data.frame(pCoa$points)
colnames(pCoaVecs) <- paste0("PCo",c(1:2))
pCoaVecs$ID <- row.names(pCoaVecs)

inMBssID <- inMBsss
inMBssID$ID <- row.names(inMBssID)
inMBssPrev <- inMBssID[,c("ID","k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Prevotellaceae.g__Prevotella.s__Prevotella_copri")]
colnames(inMBssPrev) <- c("ID","P.copri")
pCoaVecsTaxa <- merge(pCoaVecs,inMBssPrev,by="ID")

library(ggplot2)
g <- ggplot(pCoaVecsTaxa,aes(x=PCo1,y=PCo2,color=P.copri)) + 
  geom_point(size=1.25) + 
  theme_classic() + 
  scale_color_gradientn(colours = rainbow(5)) + 
  xlab(paste0('PCo1 ',round(xVar,2),' %')) + 
  ylab(paste0('PCo2 ',round(yVar,2),' %')) + 
  theme(text = element_text(size = 16))
print(g)

library(patchwork)
g+g+plot_layout(guide="collect")