1）排序图添加环境梯度等值线
主要用到vegan包的ordisurf函数，主要思路是：
- 将样点携带的环境因子信息构建样条曲线（splines）
- 样点的坐标信息和样条曲线的坐标信息统一到一个坐标系
> ?ordisurf
Function ordisurf fits a smooth surface for given variable and plots the result on ordination diagram.
ordisurf(x, y, choices = c(1, 2), knots = 10,
family = "gaussian", col = "red", isotropic = TRUE,
thinplate = TRUE, bs = "tp", fx = FALSE, add = FALSE,
display = "sites", w = weights(x, display), main, nlevels = 10,
levels, npoints = 31, labcex = 0.6, bubble = FALSE,
cex = 1, select = TRUE, method = "REML", gamma = 1,
plot = TRUE, lwd.cl = par("lwd"), ...)
# 注意此处的x， R document 的解释如下：
# For `ordisurf` an ordination configuration, either a matrix or a result known by `scores`
> LatNMDS <- metaMDS(LatOTU.R) # NMDS分析
... ...
... Procrustes: rmse 1.962234e-06 max resid 5.019163e-06
... Similar to previous best
*** Solution reached
# 提取作图数据，此处笔者主要想展示纬度梯度
> LatNMDS_data <- data.frame(
+ NMDS1 = LatNMDS$points[,1],
+ NMDS2 = LatNMDS$points[,2],
+ MAT= metadata1$MAT,
+ Lat = metadata1$Latitude,
+ Type = metadata1$Type,
+ nam = rownames(LatOTU.R)
+ )
> library(metR)
> ordi <- ordisurf(LatNMDS_data[,1:2], as.numeric(LatNMDS_data$Lat), plot = FALSE, bs = 'ds')
> ordi.grid <- ordi$grid
> ordi.sp <- expand.grid(x = ordi.grid$x, y = ordi.grid$y)
> ordi.sp$z <- as.vector(ordi.grid$z)
> ordi.sp.na <- data.frame(na.omit(ordi.sp))
> ggplot(LatNMDS_data,aes(NMDS1,NMDS2))+
+ geom_point(size=4, aes(shape = Type))+
+ annotate("text",
+ x=min(LatNMDS_data$NMDS1),
+ y=min(LatNMDS_data$NMDS2),
+ hjust=0,
+ vjust=0,
+ label=paste("Stress:",round(LatNMDS$stress,4)))+
+ geom_text(label=LatNMDS_data$nam,
+ size=2,
+ hjust=-0.2,
+ vjust=-0.5)+
+ stat_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..))+
+ geom_text_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..), size = 5)+
+ scale_color_gradient(low = "firebrick3",high = "royalblue3")
结果如下图所示：
NMDS+LatGradient
既然函数ordisurf()的输入是排序后样点的坐标值，那么其他排序方法也可以添加环境因子梯度，如PCA分析
# vegan自带数据，意将pH梯度展示
# scores()函数提取坐标值
> data("varespec")
> data("varechem")
> vare.pca <- prcomp(varespec)
> pdata = scores(vare.pca, choices=c(1,2)) %>% as.data.frame()
>
> ordi <- ordisurf(pdata[,1:2], as.numeric(varechem$pH), plot = FALSE, bs = 'ds')
> ordi.grid <- ordi$grid
> ordi.sp <- expand.grid(x = ordi.grid$x, y = ordi.grid$y)
> ordi.sp$z <- as.vector(ordi.grid$z)
> ordi.sp.na <- data.frame(na.omit(ordi.sp))
> ggplot(pdata,aes(PC1,PC2))+
+ geom_point(size=4)+
+ stat_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..))+
+ geom_text_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..), size = 5)+
+ scale_color_gradient(low = "firebrick3",high = "royalblue3")
结果如图：
PCA+pHGradient
NMDS分析原理
image
image
image
image
image
image
image
image
image
image
image
image
