diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index 9f85ba9a0e38084e5c734a152d41650f946e2952..816b33edbcf1604dfb1732938142b3d27c589d05 100644
--- a/01_Mesobromion.R
+++ b/01_Mesobromion.R
@@ -145,9 +145,9 @@ releve08trait.samp <- sample(releve08trait, round(length(releve08trait)/10), rep
species <- species %>%
rownames_to_column("RELEVE_NR") %>%
filter(RELEVE_NR %in% releve08trait.samp) %>%
- column_to_rownames("RELEVE_NR") %>%
- as.tbl() %>%
- dplyr::select(one_of(traits$species0))
+ #column_to_rownames("RELEVE_NR") %>%
+ #as.tbl() %>%
+ dplyr::select(RELEVE_NR, one_of(traits$species0))
env <- env %>%
@@ -233,6 +233,8 @@ dim(traits) #783 53
dim(env) #558 8
+
+
######4. Extract Environmental Factors ######
### CHELSA
library(raster)
@@ -281,14 +283,18 @@ write_delim(traits, path="_data/Mesobromion/traits.out.10perc.txt", delim="\t")
write_delim(env, path="_data/Mesobromion/env.10perc.txt", delim="\t")
## version without missing species
- empty <- which(colSums(species)==0)
+ empty <- which(colSums(species[,-1])==0)
traits_nozero <- traits[-empty,]
- species_nozero <- species[,-empty]
+ species_nozero <- species[,-(empty+1)]
write_delim(species_nozero , path="_data/Mesobromion/species.out.10perc_nozero.txt", delim="\t")
write_delim(traits_nozero, path="_data/Mesobromion/traits.out.10perc_nozero.txt", delim="\t")
+write_delim(species %>%
+ dplyr::select(RELEVE_NR),
+ path="_derived/Mesobromion/ReleveList.txt", delim="\t")
+
@@ -300,6 +306,7 @@ write_delim(traits_nozero, path="_data/Mesobromion/traits.out.10perc_nozero.txt"
####1. Reimport data ################################
## calculate corr between species composition matrix and traits
species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
+species_nozero <- read_delim("_data/Mesobromion/species.out.10perc_nozero.txt", delim="\t")
traits <- read_delim("_data/Mesobromion/traits.out.10perc.txt", delim="\t")
trait.labs <- read_delim("_data/Mesobromion/TraitLabels_Long.csv", delim=",") %>%
rownames_to_column("Trait.comb")
@@ -598,9 +605,10 @@ break()
CWM.wide <- species %>%
- rownames_to_column("RELEVE_NR") %>%
- reshape2::melt(.id="RELEVE_NR") %>%
- rename(species0=variable, pres=value) %>%
+ #rownames_to_column("RELEVE_NR") %>%
+ # reshape2::melt(.id="RELEVE_NR") %>%
+ pivot_longer(-RELEVE_NR, names_to = "species0", values_to = "pres") %>%
+ #rename(species0=variable, pres=value) %>%
as.tbl() %>%
filter(pres>0) %>%
arrange(RELEVE_NR) %>%
@@ -616,9 +624,20 @@ CWM.wide <- species %>%
column_to_rownames("RELEVE_NR")
#### 4. PCA Graphs ####
-#### 4.1 PCA of Y (Bealls) matrix + CWM ####
library(vegan)
-W.beals <- as.data.frame(beals(species, include=T, type=2))
+library(ggrepel)
+##from https://stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2
+circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
+ r = diameter / 2
+ tt <- seq(0,2*pi,length.out = npoints)
+ xx <- center[1] + r * cos(tt)
+ yy <- center[2] + r * sin(tt)
+ return(data.frame(x = xx, y = yy))
+}
+
+#### 4.1 PCA of Y (Bealls) matrix + CWM ####
+W.beals <- as.data.frame(beals(species_nozero %>%
+ column_to_rownames("RELEVE_NR"), include=T, type=2))
write.table(W.beals, sep="\t", file="_derived/Mesobromion/MatrixY_Beals.csv")
pca.out <- rda(W.beals)
varexpl <- round((pca.out$CA$eig)/sum(pca.out$CA$eig)*100,1)
@@ -626,54 +645,80 @@ cwms.envfit <- envfit(pca.out, CWM.wide, na.rm = T, choices = 1:5)
env.envfit <- envfit(pca.out, env %>%
dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC), choices=1:5)
+### Transform to correlations and sink envfits
+### see https://www.davidzeleny.net/anadat-r/doku.php/en:suppl_vars_examples for procedure
+scores.pca <- pca.out$CA$u
+arrow_heads <- cwms.envfit$vectors$arrows # extracts matrix of coordinates of arrow heads from ef
+r2 <- cwms.envfit$vectors$r # extracts vector of r2 for each env. variable
+cwms.cor <- arrow_heads * sqrt (r2)
+cor(CWM.wide, scores.pca[,1:5]) #double check
+
+
+arrow_heads <- env.envfit$vectors$arrows # extracts matrix of coordinates of arrow heads from ef
+r2 <- env.envfit$vectors$r # extracts vector of r2 for each env. variable
+env.cor <- arrow_heads * sqrt (r2)
+cor(env %>%
+ dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC), scores.pca[,1:5]) #double check
+
+
+sink("_derived/Mesobromion/EnvFit_CWMs_env.txt")
+cwms.envfit$vectors
+env.envfit$vectors
+sink()
+
+
+dat <- circleFun(diameter = 2, npoints = 100)
+
+myvectors <- as.data.frame(env.cor) %>%
+ rownames_to_column("mylab") %>%
+ mutate(category="Env") %>%
+ bind_rows(as.data.frame(cwms.cor) %>%
+ rownames_to_column("mylab") %>%
+ mutate(category="Trait")) %>%
+ mutate(fontface0=ifelse(mylab %in% best.5traits, "bold", "plain")) %>%
+ mutate(category=as.factor(category)) %>%
+ mutate(mycol=ifelse(category=="Trait", "Dark blue", "Dark red"))
+
+
PCA1_2 <- ggplot() +
- geom_point(data=as.data.frame(pca.out$CA$u[,1:2]),
+ geom_point(data=as.data.frame(pca.out$CA$u[,1:2]*5),
aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2),
- aes(x=0, xend=PC1, y=0, yend=PC2), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_segment(data=as.data.frame(env.envfit$vectors$arrows*.2),
- aes(x=0, xend=PC1, y=0, yend=PC2), col="Dark red", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_label(data=as.data.frame(env.envfit$vectors$arrows*.23) %>%
- rownames_to_column("Env"),
- aes(x=PC1, y=PC2, label=Env), col="Dark red", size=2,
- position = position_dodge(1) ) +
- geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
- rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
- aes(x=PC1, y=PC2, label=Trait, fontface=fontface0), col="Dark blue", size=2,
- position = position_dodge(1) ) +
+ geom_segment(data=myvectors,
+ aes(x=0, xend=PC1, y=0, yend=PC2, col=mycol),
+ arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_label(data=myvectors %>%
+ mutate_if(~is.numeric(.), ~(.)*1.2),
+ aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
+ position = position_dodge(1)) +#, segment.alpha=0.5, segment.colour=gray(0.8)) +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
theme_bw() +
- scale_y_continuous(limits=c(-0.25, 0.25)) +
- scale_x_continuous(limits=c(-0.25, 0.25)) + coord_equal() +
+ scale_color_identity() +
+ scale_y_continuous(limits=c(-1, 1)) +
+ scale_x_continuous(limits=c(-1, 1)) + coord_equal() +
theme(panel.grid = element_blank())
-(PCA3_4 <- ggplot() +
- geom_point(data=as.data.frame(pca.out$CA$u[,1:4]),
+
+PCA3_4 <- ggplot() +
+ geom_point(data=as.data.frame(pca.out$CA$u[,3:4]*5),
aes(x=PC3, y=PC4), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2),
- aes(x=0, xend=PC3, y=0, yend=PC4), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_segment(data=as.data.frame(env.envfit$vectors$arrows*.2),
- aes(x=0, xend=PC3, y=0, yend=PC4), col="Dark red", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_label(data=as.data.frame(env.envfit$vectors$arrows*.23) %>%
- rownames_to_column("Env"),
- aes(x=PC3, y=PC4, label=Env), col="Dark red", size=2,
- position = position_dodge(1)) +
- geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
- rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.5traits, 2, 1)),
- aes(x=PC3, y=PC4, label=Trait, fontface=fontface0 ), col="Dark blue", size=2,
- position = position_dodge(1)
- ) +
+ geom_segment(data=myvectors,
+ aes(x=0, xend=PC3, y=0, yend=PC4, col=mycol),
+ arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_label_repel(data=myvectors,
+ aes(x=PC3, y=PC4, label=mylab, col=mycol, fontface=fontface0), size=2,
+ position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7), segment.size = 0.5) +
xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
ylab(paste("PC4 (", varexpl[4], "%)", sep="")) +
theme_bw() +
- scale_y_continuous(limits=c(-0.25, 0.25)) +
- scale_x_continuous(limits=c(-0.25, 0.25)) +
+ scale_color_identity() +
+ scale_y_continuous(limits=c(-1, 1)) +
+ scale_x_continuous(limits=c(-1, 1)) +
coord_equal() +
theme(panel.grid = element_blank())
-)
+
PC_beals <- cowplot::plot_grid(PCA1_2,PCA3_4, nrow=1)
ggsave("_pics/Fig6_PC_Beals_1-4.png", width=10, height=5, dpi=300, last_plot())
@@ -681,23 +726,6 @@ ggsave("_pics/Fig6_PC_Beals_1-4.png", width=10, height=5, dpi=300, last_plot())
-### Transform to correlations and sink envfits
-### see https://www.davidzeleny.net/anadat-r/doku.php/en:suppl_vars_examples for procedure
-arrow_heads <- cwms.envfit$vectors$arrows # extracts matrix of coordinates of arrow heads from ef
-r2 <- cwms.envfit$vectors$r # extracts vector of r2 for each env. variable
-cwms.cor <- arrow_heads * sqrt (r2)
-
-arrow_heads <- env.envfit$vectors$arrows # extracts matrix of coordinates of arrow heads from ef
-r2 <- env.envfit$vectors$r # extracts vector of r2 for each env. variable
-env.cor <- arrow_heads * sqrt (r2)
-
-
-sink("_derived/Mesobromion/EnvFit_toCor_CWMs_env.txt")
-rbind(cwms.cor, env.cor)
-sink()
-
-
-
#### 4.2 PCA of CWMs #####
@@ -727,7 +755,11 @@ ggplot() +
#myweights <- c(rep(1/9, 9), rep(1/17, 17), rep(1, 26))
#traits.gowdis <- FD::gowdis(traits, w=myweights, ord="podani")
data2 <- traits %>%
- dplyr::select(any_of(traits.sign.alone))
+ dplyr::select(any_of(traits.sign.alone)) %>%
+ rownames_to_column("species") %>%
+ filter(species %in% colnames(species_nozero)) %>%
+ column_to_rownames("species")
+
corr1 <- cor(data2, use="pairwise.complete.obs")
pca2 <- princomp(covmat=corr1)
#plot(pca2$loadings[,2]~pca2$loadings[,1])
@@ -740,50 +772,96 @@ pca.scores<- zdat %*% e1$vectors #scaled values x vectors
pca.scores <- pca.scores[,1:5]
colnames(pca.scores)<-c('pca1','pca2','pca3','pca4','pca5') #just quickly naming the columns
-PCA.t1 <- ggplot() +
- geom_point(data=as.data.frame(pca.scores[,1:2]),
- aes(x=pca1, y=pca2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(pca2$loadings[,1:5]*8),
+dat <- circleFun(diameter = 2, npoints = 100)
+
+baseplot <- ggplot(data=as.data.frame(pca.scores*.2)) +
+ scale_color_identity() +
+ coord_equal(xlim = c(-1, 1), ylim=c(-1, 1)) +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+ theme_bw() +
+ theme(panel.grid = element_blank())
+
+PCA.t1 <- baseplot +
+ geom_point(aes(x=pca1, y=pca2), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col="Dark green"),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_label(data=as.data.frame(jitter(pca2$loadings[,1:5]*8, factor = 300)) %>%
- #geom_label(data=as.data.frame(pca2$loadings[,1:5]*8) %>%
+ geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
rownames_to_column("Trait") %>%
mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
aes(x=Comp.1, y=Comp.2, label=Trait, col="Dark green", fontface=fontface0), size=2,
- position = position_dodge(2)) +
- scale_color_identity() +
- coord_equal(xlim = c(-5, 5), ylim=c(-5, 5)) +
- xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
- ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
- theme_bw() +
- theme(panel.grid = element_blank())# +
+ position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
#ggtitle("PCoA of species-level traits")
-PCA.t2 <- ggplot() +
- geom_point(data=as.data.frame(pca.scores[,1:4]),
- aes(x=pca3, y=pca4), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(pca2$loadings[,1:5]*6),
+PCA.t2 <- baseplot +
+ geom_point(aes(x=pca3, y=pca4), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col="Dark green"),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_label(data=as.data.frame(jitter(pca2$loadings[,1:5]*6, factor = 300)) %>%
+ geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
#geom_label(data=as.data.frame(pca2$loadings[,1:5]*8) %>%
rownames_to_column("Trait") %>%
mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
aes(x=Comp.3, y=Comp.4, label=Trait, col="Dark green", fontface=fontface0), size=2,
- position = position_dodge(2)) +
- scale_color_identity() +
- coord_equal(xlim = c(-5, 5), ylim=c(-5, 5)) +
- xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
- ylab(paste("PC4 (", varexpl[4], "%)", sep="")) +
- theme_bw() +
- theme(panel.grid = element_blank())
+ position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
PC_traits <- cowplot::plot_grid(PCA.t1, PCA.t2, nrow=1)
ggsave("_pics/FigS6_PCA_Traits_1-4_only11.png", width=10, height=5, dpi=300, PC_traits)
+## alternative version of figS6, showing the species
+tmp <- as.data.frame(pca.scores[,1:4]*.2) %>%
+ mutate(species0=rownames(pca.scores)) %>%
+ mutate(species=species0) %>%
+ separate(species0, sep="_", into=c("Gen", "Spe")) %>%
+ mutate(Gen=substr(Gen, 1, 3)) %>%
+ mutate(Spe=substr(Spe, 1, 3)) %>%
+ mutate(labels=paste(Gen, Spe, sep="_"))
+
+
+PCA.t1.sp <- baseplot %+% tmp +
+ geom_text(aes(x=pca1, y=pca2, label=labels), size=2, alpha=0.7) +
+ geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
+ aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col="Dark green"),
+ arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
+ rownames_to_column("Trait") %>%
+ mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
+ aes(x=Comp.1, y=Comp.2, label=Trait, col="Dark green", fontface=fontface0), size=2,
+ position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
+
+PCA.t2.sp <- baseplot %+% tmp +
+ geom_text(aes(x=pca3, y=pca4, label=labels), size=2, alpha=0.7) +
+ geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
+ aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col="Dark green"),
+ arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
+ #geom_label(data=as.data.frame(pca2$loadings[,1:5]*8) %>%
+ rownames_to_column("Trait") %>%
+ mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
+ aes(x=Comp.3, y=Comp.4, label=Trait, col="Dark green", fontface=fontface0), size=2,
+ position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
+#PC_traits_sp <- cowplot::plot_grid(PCA.t1.sp, PCA.t2.sp, nrow=1)
+
+ggsave("_pics/FigS6_PCA_Traits_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA.t1.sp)
+ggsave("_pics/FigS6_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA.t2.sp)
+
+## Create list of species labels
+write_delim(tmp %>%
+ dplyr::select(species, label=labels),
+ path = "_derived/Mesobromion/SpeciesList.txt")
+
+
+
+
+
+
#### 4.4 PCoA of X (Fuzzy weighted) matrix ####
library(vegan)
fuzzy <- read.table("_data/Mesobromion/X_Dry-Grassland_581plots_488spp_R.txt", header=T)#, delim="\t")