客户需求根据meta的物种的NMDS结果和ARO的NMDS结果做Procrustes分析，我们之前的个性化流程中有针对meta的Procrustes分析，但是是根据PCoA结果进行分析的。
现在整理出一套R代码，可以使用任意两个NMDS结果（物种，ARO，MGE等）进行Procrustes分析，可得到普适分析图。
使用任意两个NMDS结果（物种，ARO，MGE等）进行Procrustes分析。
1. NR物种丰度文件准备：
meta分析流程中Unigenes.absolute.g.xls
2. ARO丰度文件准备：
meta分析流程中stat.ARO.absolute.xls
1. R代码：
setwd("./")
# Procrustes Anaylsis ----------
# 评估物种群落结构与环境因子间是否具存在显著的相关性（两个数据集）
# 加载R包
library(vegan)
# 加载数据
#data(varespec)
#head(varespec)
#data(varechem)
#head(varechem)
# 需要先对数据计算其样本间的距离
# 一般物种群落使用bray-curtis距离，而环境因子使用欧式距离
#spe.dist <- vegdist(varespec) # 默认Bray-Curtis
#env.dist <- vegdist(scale(varechem), "euclid")
#?procrustes
# 在进行普鲁克分析前，需要先对两个数据进行降为分析，这里使用的是NMDS
# 也可以使用PCA或者PCoA，然后进行普鲁克分析，计算显著性
#mds.s <- monoMDS(spe.dist)
#mds.e <- monoMDS(env.dist)
#mds.s
##流程NMDS方法-cjw
library(permute)
library(lattice)
library(ggplot2)
otu1=read.table("Unigenes.absolute.g.xls",header=T,row.names=1,sep="\t",comment.char="")
otu1t=t(otu1)
#dist=vegdist(otu,method="bray")
nmds1=metaMDS(otu1t)
capture.output(nmds1,file = "Stress1.txt" )
nmds_scores1=scores(nmds1, choices = c(1,2))
write.table(nmds_scores1,file="NMDS_scores1.txt")
otu2=read.table("stat.ARO.absolute.xls",header=T,row.names=1,sep="\t",comment.char="")
otu2t=t(otu2)
#dist=vegdist(otu,method="bray")
nmds2=metaMDS(otu2t)
capture.output(nmds2,file = "Stress2.txt" )
nmds_scores2=scores(nmds2, choices = c(1,2))
write.table(nmds_scores2,file="NMDS_scores2.txt")
# 以对称模式为例进行普氏分析（symmetric = TRUE）
#pro.s.e <- procrustes(mds.s,mds.e, symmetric = TRUE)
#summary(pro.s.e)
pro.s.e <- procrustes(nmds1,nmds2, symmetric = TRUE)
summary(pro.s.e)
# 一致性
plot(pro.s.e, kind = 2)
residuals(pro.s.e)
# 普氏分析中M2统计量的显著性检验
set.seed(1)
pro.s.e_t <- protest(nmds1,nmds2, permutations = 999)
pro.s.e_t
#偏差平方和(M2统计量)
M2<-round(pro.s.e_t$ss,3)
# 对应p值结果
PVAL<-pro.s.e_t$signif
# 获得x和y轴的坐标及旋转过的坐标
Pro_Y <- cbind(data.frame(pro.s.e$Yrot), data.frame(pro.s.e$X))
Pro_X <- data.frame(pro.s.e$rotation)
# 绘图
p1<-ggplot(data=Pro_Y) +
geom_segment(aes(x = X1 , y = X2, xend =NMDS1,yend =NMDS2),color = "grey", size = 1)+
#geom_segment(aes(x = X1, y = X2,
# xend = (X1 + NMDS1)/2, yend = (X2 + NMDS2)/2),
# geom_segment 绘制两点间的直线
#arrow = arrow(length = unit(0, 'cm')),
# color = "#9BBB59", size = 1) +
#geom_segment(aes(x = (X1 + NMDS1)/2, y = (X2 + NMDS2)/2,
# xend = NMDS1, yend = NMDS2),
#arrow = arrow(length = unit(0.2, 'cm')),
# color = "#957DB1", size = 1) +
geom_point(aes(X1, X2), color = "blue", size = 3, shape = 16) +
geom_point(aes(NMDS1, NMDS2), color = "red", size = 3, shape = 16) +
theme(panel.grid = element_blank(), # 绘制背景
panel.background = element_rect(color = 'black',
fill = 'transparent'),
legend.key = element_rect(fill = 'transparent'),
#axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color='black'),
axis.line = element_line(colour = "black"),
axis.title.x=element_text(colour='black', size=24),
axis.title.y=element_text(colour='black', size=24),
axis.text=element_text(colour='black',size=12)) +
labs(x = 'Dimension 1', y = 'Dimension 2', color = '') +
#labs(title="Correlation between community and environment") +
#geom_vline(xintercept = 0, color = 'gray', linetype = 2, size = 0.3) +
#geom_hline(yintercept = 0, color = 'gray', linetype = 2, size = 0.3) +
#geom_abline(intercept = 0, slope = Pro_X[1,2]/Pro_X[1,1], size = 0.3) +
#geom_abline(intercept = 0, slope = Pro_X[2,2]/Pro_X[2,1], size = 0.3) +
#annotate('text', label = paste("Procrustes analysis:M2 = ",M2,",p-value = ",PVAL,sep=""),
# x = -0.3, y = 0.3, size = 4,hjust = 0) +
theme(plot.title = element_text(size=14,colour = "black",
hjust = 0.5,face = "bold"))
p1
ggsave("PACo.pdf",height=12,width=15)
dev.off()
2. 运行方式：
/PUBLIC/software/public/System/R-2.15.3/bin/R -f procrustes.R
完整代码路径：
/TJPROJ5/META_ASS/16s/chenjiawei/script/META/gxh/Procrustes/procrustes.R