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

Added PCA graphs

parent 02aa9b65
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 63 additions and 90 deletions
01_Mesobromion.R
+ 63
− 90
View file @ 80d08026
...@@ -603,22 +603,28 @@ CWM.wide <- species %>% ...@@ -603,22 +603,28 @@ CWM.wide <- species %>%
#### 4.1 PCA of Y (Bealls) matrix #### #### 4.1 PCA of Y (Bealls) matrix ####
library(vegan) library(vegan)
W.beals <- as.data.frame(beals(species, include=T, type=2)) W.beals <- as.data.frame(beals(species, include=T, type=2))
#save(W.beals, file="_derived/Mesobromion/MatrixY_Beals.RData")
write.table(W.beals, sep="\t", file="_derived/Mesobromion/MatrixY_Beals.csv") write.table(W.beals, sep="\t", file="_derived/Mesobromion/MatrixY_Beals.csv")
pca.out <- rda(W.beals) pca.out <- rda(W.beals)
varexpl <- round((pca.out$CA$eig)/sum(pca.out$CA$eig)*100,1) varexpl <- round((pca.out$CA$eig)/sum(pca.out$CA$eig)*100,1)
cwms.envfit <- envfit(pca.out, CWM.wide, na.rm = T, choices = 1:3) cwms.envfit <- envfit(pca.out, CWM.wide, na.rm = T, choices = 1:3)
env.envfit <- envfit(pca.out, env %>%
dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC), choices=1:3)
PCA1_2 <- ggplot() + 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]),
aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) + aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2), 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) + aes(x=0, xend=PC1, y=0, yend=PC2), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>% geom_segment(data=as.data.frame(env.envfit$vectors$arrows*.2),
rownames_to_column("Trait"), aes(x=0, xend=PC1, y=0, yend=PC2), col="Dark red", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
aes(x=PC1, y=PC2, label=Trait), col="Dark blue", size=2, 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) ) + position = position_dodge(1) ) +
geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
rownames_to_column("Trait"),
aes(x=PC1, y=PC2, label=Trait), col="Dark blue", size=2,
position = position_dodge(1) ) +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) + xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep="")) + ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
theme_bw() + theme_bw() +
...@@ -631,6 +637,12 @@ PCA1_2 <- ggplot() + ...@@ -631,6 +637,12 @@ PCA1_2 <- ggplot() +
aes(x=PC1, y=PC3), pch="+", size=2, alpha=0.8) + aes(x=PC1, y=PC3), pch="+", size=2, alpha=0.8) +
geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2), geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2),
aes(x=0, xend=PC1, y=0, yend=PC3), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) + aes(x=0, xend=PC1, y=0, yend=PC3), 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=PC3), 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=PC3, label=Env), col="Dark red", size=2,
position = position_dodge(1)) +
geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>% geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
rownames_to_column("Trait"), rownames_to_column("Trait"),
aes(x=PC1, y=PC3, label=Trait), col="Dark blue", size=2, aes(x=PC1, y=PC3, label=Trait), col="Dark blue", size=2,
...@@ -673,29 +685,55 @@ ggplot() + ...@@ -673,29 +685,55 @@ ggplot() +
theme(panel.grid = element_blank()) + theme(panel.grid = element_blank()) +
title("PCA of CWMs") title("PCA of CWMs")
#### 4.3 PCA of individual traits ##### #### 4.3 PCoA of individual traits #####
trait.pca <- vegan::rda(as.matrix(traits %>% # assign weights to traits [Growth forms and root structures are coded as binary, but should be downweighted]
filter(!is.na(rowSums(.)))), na.action=na.omit) myweights <- c(rep(1/9, 9), rep(1/17, 17), rep(1, 26))
varexpl.t <- round(trait.pca$CA$eig/sum(trait.pca$CA$eig)*100,1) traits.gowdis <- FD::gowdis(traits, w=myweights, ord="podani")
pcoa.species <- wcmdscale(traits.gowdis, eig=T)
ggplot() + varexpl <- round((pcoa.species$eig)/sum(pcoa.species$eig)*100,1)
geom_point(data=as.data.frame(trait.pca$CA$u[,1:2]), traits.envfit <- envfit(pcoa.species, traits %>%
aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) + dplyr::select(any_of(traits.sign.alone)),
geom_segment(data=as.data.frame(trait.pca$CA$v*2) %>% na.rm = T, choices = 1:3)
rownames_to_column("Trait") %>%
mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")), set.seed(14)
aes(x=0, xend=PC1, y=0, yend=PC2, col=in.try), arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) + PCoA.t1 <- ggplot() +
geom_text(data=as.data.frame(trait.pca$CA$v*2.1) %>% geom_point(data=as.data.frame(pcoa.species$points[,1:2]),
rownames_to_column("Trait") %>% aes(x=Dim1, y=Dim2), pch="+", size=2, alpha=0.8) +
mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")), geom_segment(data=as.data.frame(traits.envfit$vectors$arrows*0.2),
aes(x=PC1, y=PC2, label=Trait, col=in.try), size=3 ) + aes(x=0, xend=Dim1, y=0, yend=Dim2, col="Dark green"),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_label(data=as.data.frame(jitter(traits.envfit$vectors$arrows*0.22, factor = 300)) %>%
rownames_to_column("Trait"),
aes(x=Dim1, y=Dim2, label=Trait, col="Dark green"), size=2,
position = position_dodge(2)) +
scale_color_identity() + scale_color_identity() +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) + coord_equal(xlim = c(-.25, .25), ylim=c(-.25, .25)) +
ylab(paste("PC2 (", varexpl[2], "%)", sep="")) + xlab(paste("PCo1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PCo2 (", varexpl[2], "%)", sep="")) +
theme_bw() + theme_bw() +
theme(panel.grid = element_blank()) + theme(panel.grid = element_blank())# +
title("PCA of species-level traits") #ggtitle("PCoA of species-level traits")
PCoA.t2 <- ggplot() +
geom_point(data=as.data.frame(pcoa.species$points[,1:3]),
aes(x=Dim1, y=Dim3), pch="+", size=2, alpha=0.8) +
geom_segment(data=as.data.frame(traits.envfit$vectors$arrows*0.2),
aes(x=0, xend=Dim1, y=0, yend=Dim3, col="Dark green"),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_label(data=as.data.frame(traits.envfit$vectors$arrows*0.22) %>%
rownames_to_column("Trait"),
aes(x=Dim1, y=Dim3, label=Trait, col="Dark green"), size=2,
position = position_dodge(2)) +
scale_color_identity() +
coord_equal(xlim = c(-.25, .25), ylim=c(-.25, .25)) +
xlab(paste("PCo1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PCo3 (", varexpl[3], "%)", sep="")) +
theme_bw() +
theme(panel.grid = element_blank())
PC_traits <- cowplot::plot_grid(PCoA.t1, PCoA.t2, nrow=1)
ggsave("_pics/PCoA_Traits.png", width=10, height=5, dpi=300, PC_traits)
#### 4.4 PCoA of X (Fuzzy weighted) matrix #### #### 4.4 PCoA of X (Fuzzy weighted) matrix ####
...@@ -770,69 +808,4 @@ ggsave(filename="_pics/PlotDistribution.png", height=8, width=6, dpi=300, basema ...@@ -770,69 +808,4 @@ ggsave(filename="_pics/PlotDistribution.png", height=8, width=6, dpi=300, basema
### renaming variables
"LF_Macroph = LEB_F_Makrophanerophyt
LF_Nanoh = LEB_F_Nanophanerophyt
LF_Hemicr = LEB_F_Hemikryptophyt
LF_Geoph = LEB_F_Geophyt
LF_Hemiph = LEB_F_Hemiphanerophyt
LF_Theropf = LEB_F_Therophyt
LF_Hydroph = LEB_F_Hydrophyt
LF_Pseudoph = LEB_F_Pseudophanerophyt
LF_Chamaeph = LEB_F_Chamaephyt
LL_Pluri_poll = LEB_D_plurienn_pollakanth
LL_Pluri_hapa = LEB_D_plurienn_hapaxanth
LL_Annual = LEB_D_annuell
LL_Bienn = LEB_D_bienn
VS_Absent = V_V_VER_absent
VS_RootSprout = V_VER_Wurzelspross
VS_Runners = V_VER_Ausläufer
VS_Rhizon = V_VER_Rhizom
VS_Tuber = V_VER_Innovationsknopse.mit.Wurzelknolle
VS_StorageRoot = V_VER_Innovationsknospe.mit.Speicherwurzel
VS_RunningTuber = V_VER_Ausläuferknolle
VS_Brutsprösschen = V_VER_Brutsprösschen ## How to translate?
VS_Fragment = V_VER_Fragmentation
VS_Turio = V_VER_Turio
24 V_VER_Sprossknolle
25 V_VER_phyllogener_Spross
26 V_VER_Rhizompleiokorm
27 V_VER_Zwiebel
28 V_VER_Ausläuferrhizom
29 V_VER_Ausläuferzwiebel
30 V_VER_Bulbille
31 ROS_T_rosettenlose.Pflanzen
32 ROS_T_Halbrosettenpflanze
33 ROS_T_Ganzrosettenpflanzen
34 BL_AUSD_immergrün
35 BL_AUSD_sommergrün
36 BL_AUSD_vorsommergrün
37 BL_AUSD_überwinternd_grün
38 BL_ANAT_skleromorph
39 BL_ANAT_mesomorph
40 BL_ANAT_hygromorph
41 BL_ANAT_hydromorph
42 BL_ANAT_blattsukkulent
43 BL_ANAT_helomorph
44 BL_BEG
45 BL_END
46 BL_DAU
47 LeafArea
48 SLA
49 LeafC.perdrymass
50 LeafN
51 LeafP
52 PlantHeight
53 SeedMass
54 Seed.length
55 LDMC
56 LeafNperArea
57 LeafNPratio
58 Leaf.delta.15N
59 Seed.num.rep.unit
60 Leaffreshmass
61 Disp.unit.leng
62 BLU_KL
63 REPR_T
64 STRAT_T"
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment