Script introduction
Find in the following lines the code corresponding to the analysis of the article: Heterogeneity in the isolation of patches may be essential for the action of metacommunity mechanisms.
All the suplementary material images, original scripts and used functions can be found in:
Article DOI:
Direct link to repositpory: https://github.com/Metacommunity-Lab/FEE_Coal_Rand_Pondscape
To properly follow the following steps we recommend to read the entire work to better contextualise the use and limitations of the current codes and functions.
The scripts can be found in the following repository
Link to the function: https://github.com/Metacommunity-Lab/FEE_Coal_Rand_Pondscape
Link to supplementary material:
# Load basic functions to run all the simulations
library(dplyr); library(tibble); library(ggplot2); library(gridExtra)Main developedfunctoins
To develpe and carry all the analysis that readers can found in the previously cited paper we developed two main functions that define our two main analysis:
1. Pondscape diversity along dispersal and centrality gradients (mainly providing Figure 3A and B columns)
2. Pondscape randomization and diveristy comparisons (mainly providing Figure 3C and D colums)
For each one of these analysis we provide in the repository the scripts developing the corresponding functions: 1) FEE_IterationModel_Ponderful.R and 2) FEE_Random_landscape_Ponderful.R.
In these scripts there is a detailed explanations of what each line is conducting and which is the specific objective deboted to each function. We begun the analysis by charging these two functions as well as the .RData archive where all the necessary information is contained (xy coordinates of each water body from each one of the studied pondscapes).
# Charging functoins
source("FEE_Random_landscape_Ponderful.R")
source("FEE_IterationModel_Ponderful.R")
# We load all the necessary information for each studied pondscape
load("Espacio_Retromed.RData")
CHARCOS_ROCHA <- read.delim("CHARCOS_ROCHA.txt")
Rocha<-CHARCOS_ROCHA
Rocha <- cbind(Rocha,0,1,100)Lines of analysis for each Pondscape
For each pondscape we run the same set of code lines where we first define the set of requirements for the two functions that we use. Then, we run the two functions that can take several hours to finish as the coalescent process is based on several iterations were we fill the communities with individuals of species until all the community spaces are full (see methods section from the above manuscript for a detailed explanation). The following lines can also be found in the script named Lines_Analysis.R contained in the main github repository of this project (see above).
Albera
# Seting elements for each function
Albera<-cbind(ALBERA,100)
Albera<-Albera[,-4]
colnames(Albera)[8]<-"J"
Albera<-Albera[,c(1:3, 5,7,6,4,8)]
M.distance.Albera<-as.matrix(dist(Albera[,2:3]))
Meta.comm.Albera<-Meta.comm.Albera
J.Albera<-Albera[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Albera$grado
Cc<-Albera$cc
Bt<-Albera$bc
N<-nrow(M.distance.Albera)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Albera<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Albera, m.pool=0.01, Js=J.Albera, id.module=Albera[,7],
M.dist=M.distance.Albera, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Albera,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Albera))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Albera<-rbind(out.d50.obs.landscape.Albera, out.t2)
}
out.d50.obs.landscape.Albera
save.image("~/Dropbox/H2020/RETROMED/Frontiers_Diciembre_2022/Frontiers_Respuestas/Retromed_H2020_febrero_2023.RData")
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Albera_Febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Albera, m.pool=0.01,
Js=Albera[,8],
id.module=Albera[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Albera[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1,
comm.fixed=Meta.comm.Albera,Lottery=F, it=0,
prop.dead.by.it=0,id.obs=1:nrow(Albera))Giara
# Seting elements for each function
Giara<-cbind(GIARA,100)
Giara<-Giara[,-4]
colnames(Giara)[8]<-"J"
Giara<-Giara[,c(1:3, 5,7,6,4,8)]
M.distance.Giara<-as.matrix(dist(Giara[,2:3]))
Meta.comm.Giara<-Meta.comm
J.Giara<-Giara[,8]
#D50s<-seq(50,3000,,50)
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Giara$grado
Cc<-Giara$cc
Bt<-Giara$bc
N<-nrow(M.distance.Giara)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Giara<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Giara, m.pool=0.01, Js=J.Giara, id.module=Giara[,7],
M.dist=M.distance.Giara, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Giara,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Giara))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Giara<-rbind(out.d50.obs.landscape.Giara, out.t2)
}
out.d50.obs.landscape.Giara
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Giara_febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Giara, m.pool=0.01, Js=Giara[,8],
id.module=Giara[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Giara[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1, comm.fixed=Meta.comm.Giara,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Giara))Guils
# Seting elements for each function
Guils<-cbind(GUILS,100)Figure_3_Albera.png
Guils<-Guils[,-4]
colnames(Guils)[8]<-"J"
Guils<-Guils[,c(1:3, 5,7,6,4,8)]
M.distance.Guils<-as.matrix(dist(Guils[,2:3]))
Meta.comm.Guils<-Meta.comm.Albera
J.Guils<-Guils[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Guils$grado
Cc<-Guils$cc
Bt<-Guils$bc
N<-nrow(M.distance.Guils)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Guils<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Guils, m.pool=0.01, Js=J.Guils, id.module=Guils[,7],
M.dist=M.distance.Guils, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Guils,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Guils))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Guils<-rbind(out.d50.obs.landscape.Guils, out.t2)
}
out.d50.obs.landscape.Guils
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Guils_febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 1000,Meta.pool=Meta.comm.Guils, m.pool=0.01,
Js=Guils[,8],
id.module=Guils[,7],
D50.min=1, D50.max=2000,m.max=1,M.X.Y=Guils[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1,
comm.fixed=Meta.comm.Guils,
Lottery=F, it=0, prop.dead.by.it=0,
id.obs=1:nrow(Guils))PNAE
# Seting elements for each function
Pnae<-cbind(PNAE,100)
Pnae<-Pnae[,-4]
colnames(Pnae)[8]<-"J"
Pnae<-Pnae[,c(1:3, 5,7,6,4,8)]
M.distance.Pnae<-as.matrix(dist(Pnae[,2:3]))
Meta.comm.Pnae<-Meta.comm
J.Pnae<-Pnae[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Pnae$grado
Cc<-Pnae$cc
Bt<-Pnae$bc
N<-nrow(M.distance.Pnae)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Pnae<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Pnae, m.pool=0.01, Js=J.Pnae, id.module=Pnae[,7],
M.dist=M.distance.Pnae, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Pnae,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Pnae))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Pnae<-rbind(out.d50.obs.landscape.Pnae, out.t2)
}
out.d50.obs.landscape.Pnae
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Pnae_Febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Pnae, m.pool=0.01, Js=Pnae[,8],
id.module=Pnae[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Pnae[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1, comm.fixed=Meta.comm.Pnae,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Pnae))Vila nova
# Seting elements for each function
Vilavona<-cbind(VILAVONA,100)
Vilavona<-Vilavona[,-4]
colnames(Vilavona)[8]<-"J"
Vilavona<-Vilavona[,c(1:3, 5,7,6,4,8)]
M.distance.Vilavona<-as.matrix(dist(Vilavona[,2:3]))
Meta.comm.Vilavona<-Meta.comm.Albera
J.Vilavona<-Vilavona[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Vilavona$grado
Cc<-Vilavona$cc
Bt<-Vilavona$bc
N<-nrow(M.distance.Vilavona)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Vilavona<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Vilavona, m.pool=0.01, Js=J.Vilavona, id.module=Vilavona[,7],
M.dist=M.distance.Vilavona, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Vilavona,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Vilavona))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Vilavona<-rbind(out.d50.obs.landscape.Vilavona, out.t2)
}
out.d50.obs.landscape.Vilavona
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Vilavona_Febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Vilavona,
m.pool=0.01, Js=Vilavona[,8],
id.module=Vilavona[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Vilavona[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1,
comm.fixed=Meta.comm.Vilavona,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Vilavona))Rocha
# Seting elements for each function
colnames(Rocha)[6:8]<-c("bc","modules", "J")
Rocha[1:10,7]<-2
M.distance.Rocha<-as.matrix(dist(Rocha[,2:3]))
Meta.comm.Rocha<-Meta.comm.Albera
J.Rocha<-Rocha[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Rocha$grado
Cc<-Rocha$cc
Bt<-Rocha$bc
N<-nrow(M.distance.Rocha)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Rocha<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Rocha, m.pool=0.01, Js=J.Rocha, id.module=Rocha[,7],
M.dist=M.distance.Rocha, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Rocha,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Rocha))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Rocha<-rbind(out.d50.obs.landscape.Rocha, out.t2)
}
out.d50.obs.landscape.Rocha
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Rocha_febrero2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Rocha, m.pool=0.01, Js=Rocha[,8],
id.module=Rocha[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Rocha[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1, comm.fixed=Meta.comm.Rocha,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Rocha))Lines of ploting and results representation for each Pondscape
Once all the results have been produced and all the simulations carried we extract the specific outputs of alpha diversity and betadiversity in order to plot them accrodingly and to represent the article figures. The following lines can also be found in the script named Lines_PLots.R contained in the main github repository of this project (see above).
Figure 1
Here we represent the pondscapes analysed in this study. Five European pondscapes (Albera, Clots de Guils, Empordà Wetlands (PNAE) and Vila Nova de Milfontes) and one Uruguayan (Rocha).
xy.ALB <- read.table(file = "Pondscapes_loc/CHARCOS_ALBERA.txt",sep = "\t", dec =",",header = T)
xy.GIA <- read.table(file = "Pondscapes_loc/CHARCOS_GIARA.txt",sep = "\t", dec =",",header = T)
xy.GUI <- read.table(file = "Pondscapes_loc/CHARCOS_GUILS.txt",sep = "\t", dec =",",header = T)
xy.PNAE <- read.table(file = "Pondscapes_loc/CHARCOS_PNAE.txt",sep = "\t", dec =",",header = T)
xy.VMIL <- read.table(file = "Pondscapes_loc/CHARCOS_VILANOVAMILFONTES.txt",sep = "\t", dec =",",header = T)
xy.ROC <- read.table(file = "Pondscapes_loc/CHARCOS_ROCHA.txt",sep = "\t", dec =".",header = T)
# Check structure and correct PNAE
str(xy.ALB)
str(xy.GIA)
str(xy.GUI)# Remember that they are repeated
xy.GUI<- xy.GUI[-which(duplicated(xy.GUI[,2:3])==T),]
str(xy.PNAE)# Remember that they are repeated
xy.PNAE<- xy.PNAE[-which(duplicated(xy.PNAE[,2:3])==T),]
str(xy.VMIL)
str(xy.ROC)
c(nrow(xy.ALB),nrow(xy.GIA),nrow(xy.GUI),nrow(xy.PNAE),nrow(xy.VMIL),nrow(xy.ROC))
c(542,160,86,255,146)
# Creation of minimum spanning tree (one link per pond)
GG_ALB <- mst(as.matrix(dist(xy.ALB)))
GG_GIA <- mst(as.matrix(dist(xy.GIA)))
GG_GUI <- mst(as.matrix(dist(xy.GUI)))
GG_PNAE <- mst(as.matrix(dist(xy.PNAE)))
GG_VMIL <- mst(as.matrix(dist(xy.VMIL)))
GG_ROC <- mst(as.matrix(dist(cbind(xy.ROC[,2],xy.ROC[,3]))))
# List creation
xy.values <- list(xy.ALB[,2:3],xy.PNAE[,2:3],xy.GIA[,2:3],xy.VMIL[,1:2], xy.GUI[,2:3], xy.ROC[,2:3])
GG.values <- list(GG_ALB,GG_PNAE, GG_GIA, GG_VMIL, GG_GUI, GG_ROC)
Pondscapes_names <- c("Albera","PNAE","Giara", "Vila Nova","Guils", "Rocha")
# Nice colours? CUNILLERA_palette is what you need
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R") # Also found in https://github.com/Cunillera-Montcusi/Cunillera_Pallete
library(sna)
gg_plot_network <- list()
for (u in 1:length(xy.values)) {
detach("package:sna", unload = TRUE)
library(igraph)
# In order to plot the network in a ggplot way do the follwing
weight_degree <- Pondscapes_outputs[[u]][which(Pondscapes_outputs[[u]][,1]==max_diff[u]),7]
detach("package:igraph", unload = TRUE)
library("sna")
n<- network(GG.values[[u]], directed=F, diag=F)
n %v% "family" <- weight_degree # Family is an standard name...
# ALL THESE LINES MUST BE RUNNED!!!
#from here_____________________________________________________________________
gg_plot_network[[u]] <- ggplot(n, layout=as.matrix(xy.values[[u]]),
aes(x = x, y = y, xend = xend, yend = yend))+
geom_edges(color = "grey30", size=0.5) +
geom_nodes(aes(fill=family, size=family) ,color ="black" ,shape=21, alpha=.55)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(x="",y="",title = Pondscapes_names[[u]])+
theme(axis.line = element_blank(),
plot.title = element_text(size=35, colour = "black",hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = "none",
panel.background = element_blank())
legend_plots<- get_legend(ggplot(n, layout=as.matrix(xy.values[[u]]),
aes(x = x, y = y, xend = xend, yend = yend))+
geom_edges(color = "black", size=0.5) +
geom_nodes(aes(fill=family) ,color ="grey50" ,shape=21, alpha=.55)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Weighted degree",
guide = guide_colorbar(ticks = T,
barheight = 5,
barwidth = 1))+
theme_classic()+
theme(legend.direction = "vertical",
legend.box="vertical"))
}
Figure 3
This is figure representing the main results obtained in the presented work were the interplay between landscape configuration and species dispersal ability as determinants of metacommunity biodiversity is plotted. Rows correspond to metacommunities of water bodies. A, B: Alpha and beta diversities along a dispersal ability gradient. Dispersal ability in x-axes is the distance at which dispersal falls to half its maximum rate (d_50). The colours of points represent the isolation-centrality gradient (weighted degree centrality) of communities. The black curves represent the set of more isolated (dashed) and central (continuous) communities in the pondscape. Inset figure is the ratio for alpha or beta diversity in the most peripheral vs. the most central community in the gradient of dispersal ability—i.e., 〖Ratio〗_alpha and 〖Ratio〗_beta. C, D: Importance of pondscape configuration for alpha and beta diversity along the gradient of dispersal abilities. Z-values correspond to standardised deviations of diversity predicted for the real pondscape structure (S ̅_obs) and the average diversity expected by the null model prediction when the location of communities was randomised (S ̅_null): Z-value=((S ̅_obs-S ̅_null))⁄(SD(S ̅_null)), for each dispersal ability (d_50).
# Nice colours? CUNILLERA_palette is what you need
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R")
Pondscapes_outputs <- list(out.d50.obs.landscape.Albera,
out.d50.obs.landscape.Pnae,
out.d50.obs.landscape.Giara,
out.d50.obs.landscape.Vilavona,
out.d50.obs.landscape.Guils,
out.d50.obs.landscape.Rocha)
Pondscapes__rand_outputs <- list(random_vs_real_D50.Albera_Febrero_2023,
random_vs_real_D50.Pnae_Febrero_2023,
random_vs_real_D50.Giara_febrero_2023,
random_vs_real_D50.Vilavona_Febrero_2023,
random_vs_real_D50.Guils_Febrero_2023,
random_vs_real_D50.Rocha_febrero2023)
Pondscapes_names <- c("Albera", "PNAE", "Giara", "Vila Nova", "Guils","Rocha")
final_plot_v1 <- list()
final_plot_v2 <- list()
supp_CV_plot <- list()Figure 3 Albera (first row)
#Albera________________________________________________________________________________________________________####
p=1
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=4, Smax=30, q=3, d50=150))
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=32, q=3, d50=100))
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>% ggplot(aes(d,S))+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(title="A)",subtitle = Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(40), xmax = log10(250), ymin = 20, ymax = 37)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=500))
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-3, d50=150))
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>% arrange(desc(d)) %>% ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(title="B)",subtitle = " ", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(40), xmax = log10(250), ymin = .5, ymax = .8)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = "C)", subtitle = " ")+
geom_hline(yintercept=0, linetype="dashed")+scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title="D)", subtitle = " ")+
theme_bw()+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- data.frame()
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Albera",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)
Figure 3 PNAE (second row)
#PNAE________________________________________________________________________________________________________####
p=2
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.01))
ii.max.g<-which(a[,6]>quantile(a[,6],.99))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=0, Smax=30, q=3, d50=500))
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=0, Smax=20, q=3, d50=100))
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle=Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(520), xmax = log10(2900), ymin = 1, ymax = 22)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=250) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-1, d50=200) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(510), xmax = log10(3100), ymin = .7, ymax = .99)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+#scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+#scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "PNAE",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)
Figure 3 Giara (third row)
#Giara________________________________________________________________________________________________________####
p=3
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=0, Smax=25, q=3, d50=500) )
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=35, q=1, d50=100) )
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle = Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(900), xmax = log10(3200), ymin = 0, ymax = 22)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=500) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=10, Bmax=1, q=-3, d50=50) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(40), xmax = log10(250), ymin = .42, ymax = .74)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+#scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+#scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Giara",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)
Figure 3 Vila nova (foruth row)
#Vilanova________________________________________________________________________________________________________####
p=4
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=4, Smax=30, q=3, d50=200) )
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=32, q=3, d50=100) )
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle=Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(400), xmax = log10(3000), ymin = 1, ymax = 27)
#### BETA
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=250) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-3, d50=150) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for Beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(500), xmax = log10(3000), ymin = .53, ymax = .99)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Vila Nova",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)