Skip to content
Snippets Groups Projects
Commit 8cbd3652 authored by Francesco Sabatini's avatar Francesco Sabatini
Browse files

Corrected problem with Fig 6 envfit. New version of Fig S6

parent 9c1db279
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 166 additions and 88 deletions
01_Mesobromion.R
+ 166
− 88
View file @ 8cbd3652
......@@ -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")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment