PCA with vegan

PCA_with_vegan.R

source('prep_data.R')
library(vegan)
library(RColorBrewer)
 
# select env/process data
meta <- data[ , 1:( which( names(data) == "amplib_id" )) ]
env <- data[ , ( which( names(data) == "primary_sed" )):( which( names(data) == "otu_0" ) - 1) ]
otu <- as.matrix(data[, (which(names(data) == "otu_0" )):ncol(data)  ])
 
summary(env)
 
# PCA
hell.otu <- decostand(otu, method="hell")
pca.hell.otu <- rda(hell.otu)
 
screeplot(pca.hell.otu, type = 'line')
 
scores(pca.hell.otu, choices = 1:4, display = "species", scaling = 0)
pca.hell.otu$CA$v.eig  # eigen vectors
pca.hell.otu$CA$v   # unconstrained species scores
pca.hell.otu$CA$u  # unconstrained site scores
per.var <- round(pca.hell.otu$CA$eig / pca.hell.otu$CA$tot.chi * 100, 2)[1:5]
 
 
par(mfrow= c(1, 1))
 
# label point by metadata
plot(pca.hell.otu, type = "n", main= "Sites from PCA (hell corr)" )
points(pca.hell.otu, display = "sites", cex = 0.8, pch=1, col="black")
points(pca.hell.otu, display = "spec", pch=3, cex=0.7, col="blue")
 
with(meta, ordispider(pca.hell.otu, country, col="blue"))
with(meta, ordihull(pca.hell.otu, country, col="blue", label=TRUE))
with(meta, ordiellipse(pca.hell.otu, country, col="blue", label=TRUE))
 
 
# shows that countries and plants cluster significantly
plot(pca.hell.otu, type = "n", main= "Sites from PCA (hell corr)" )
points(pca.hell.otu, display = "sites", cex = 0.8, pch=1, col="black")
with(meta, ordispider(pca.hell.otu, country, col="blue"))
with(meta, ordihull(pca.hell.otu, country, col="blue", label=TRUE))
 
with(meta, ordispider(pca.hell.otu, plant, col="blue"))
with(meta, ordihull(pca.hell.otu, plant, col="blue", label=TRUE))
 
mf <- envfit(pca.hell.otu, meta, permutations= 999)
 
# label points with country code
plot(pca.hell.otu, type = "n", main= "Sites from PCA (hell corr)" )
text(pca.hell.otu, display = "sites", labels=meta$country, col="black")
 
# color points by country and add legend
numcols <- length(levels(meta$country))
cols <- brewer.pal(numcols, "Set1")
pal <- colorRampPalette(cols)
 
pl <- plot( pca.hell.otu, type = "n", )
points(pca.hell.otu, display = "sites",
       cex = 0.8, pch=20,
       col= cols[as.numeric(meta$country)] )
legend(-1, 1, levels(meta$country), pch=20, col = cols[1:length(meta$country)])
points(pca.hell.otu, display = "sp",
              cex = 0.8, pch=3, col= "blue")
identify(pl, "sp", labels=sp.names)
 
 
# plot sites against envfit
col.inf <- c('ebpr', 'Ca_add', 'glycerol_add', 'EtOH_add', 'MeOH_add')
env.red <- env[ , !(names(env) %in% col.inf) ]
ef <- envfit(pca.hell.otu, env.red, permutations=999)
ef
plot(pca.hell.otu, display="sites", type="p")
 
pl <- plot( pca.hell.otu, type = "n", )
points(pca.hell.otu, display = "sites",
       cex = 0.8, pch=20,
       col= cols[as.numeric(meta$country)] )
legend(1, 1, levels(meta$country), pch=20, col = cols[1:length(meta$country)])
plot(ef, p.max = 0.001, col="red")
 
with(meta, ordispider(pca.hell.otu, country, col="blue"))
with(meta, ordihull(pca.hell.otu, country, col="blue", label=TRUE))