# Libraries ----
library(here)
library(arm)
library(DHARMa)
library(dplyr)
library(tidyr)
library(cxr)
library(igraph)
library(dplyr)
library(tidyr)
library(scales)
library(ggplot2)
###########################################
#### Prepare data: HWA(CONTROL) ###########
###########################################
#Prepare data ----
datacompcontrol<-read.table(here("MCT_analysis", "MCT_data_control.txt"),header = T, stringsAsFactors = T)
datacomp3control<-datacompcontrol
datacomp3control$ID<-interaction(datacomp3control$PLOT,datacomp3control$COLUMN,datacomp3control$ROW, datacomp3control$PLOIDY)
datacomp3control$ID<-make.unique(as.character(datacomp3control$ID))
datacomp4control<-datacomp3control[,c("Biomasa","ID", "PLOIDY", "Combinacion")]
datacomp5control<-datacomp4control %>%
separate(Combinacion, into = c("Combinacion1", "Combinacion2"), sep = "_")
datacomp6control <- datacomp5control %>%
mutate(PLOIDY_C = if_else(PLOIDY == Combinacion1, Combinacion2, Combinacion1))
datacomp7control<-datacomp6control[,c("ID", "Biomasa","PLOIDY", "PLOIDY_C")]
colnames(datacomp7control)[2]<-"fitness"
datacomp7control$PLOIDY<-paste0("P",datacomp7control$PLOIDY)
datacomp7control$PLOIDY_C<-paste0("P",datacomp7control$PLOIDY_C)
datacomp7control$fitness<-round(datacomp7control$fitness*1000,0)
datacomp8control<-datacomp7control %>%
mutate(`P2x` = as.integer(grepl('\\bP2x\\b', PLOIDY_C)),
`P4x` = as.integer(grepl('4x', PLOIDY_C)),
`P6x` = as.integer(grepl('6x', PLOIDY_C)),
`P12x` = as.integer(grepl('12x', PLOIDY_C)))
datacomp9control<-datacomp8control[,-c(3,4)]
datacomp9control
datalistcontrol <- split(datacomp9control, f = datacomp8control$PLOIDY)
datalist_cleancontrol<-lapply(datalistcontrol, na.omit)
#Interaction coefficients and niche differences
datacontrol <- datalist_cleancontrol
for(i in 1:length(datacontrol)){
datacontrol[[i]] <- datacontrol[[i]][,2:length(datacontrol[[i]])]
}
focal_column <- names(datacontrol)
model_family <- "RK"
optimization_method <- "bobyqa" # we use a bounded method
alpha_form <- "pairwise"
lambda_cov_form <- "none"
alpha_cov_form <- "none"
initial_values = list(lambda = 10,
alpha_intra = 0.1,
alpha_inter = 0.1)
lower_bounds = list(lambda = 0,
alpha_intra = 0.01,
alpha_inter = 0.01)
upper_bounds = list(lambda = 2000,
alpha_intra = 1,
alpha_inter = 1)
fixed_terms <- NULL
bootstrap_samples <- 0
all.sp.fitcontrol <- cxr_pm_multifit(data = datacontrol,
model_family = model_family,
focal_column = focal_column,
optimization_method = optimization_method,
alpha_form = alpha_form,
lambda_cov_form = lambda_cov_form,
alpha_cov_form = alpha_cov_form,
initial_values = initial_values,
lower_bounds = lower_bounds,
upper_bounds = upper_bounds,
fixed_terms = fixed_terms,
bootstrap_samples = bootstrap_samples)
all.sp.fitcontrol #Intrinsic growth rate and Interaction coefficients
niche_overlap_all_pairscontrol <- niche_overlap(cxr_multifit = all.sp.fitcontrol)
niche_overlap_all_pairscontrol #Niche overlap
niche_overlap_all_pairscontrol$MCT_niche_diff <- 1 - niche_overlap_all_pairscontrol$niche_overlap_MCT
niche_overlap_all_pairscontrol$MCT_niche_diff #Niche differences
# Average Fitness differences
avg_fitness_diff_all_pairscontrol <- avg_fitness_diff(cxr_multifit = all.sp.fitcontrol)
avg_fitness_diff_all_pairscontrol <- avg_fitness_diff_all_pairscontrol[complete.cases(
avg_fitness_diff_all_pairscontrol),] # as with niche overlap, remove NA values
avg_fitness_diff_all_pairscontrol
# Competitive ability
competitive_ability_all_pairscontrol <- competitive_ability(cxr_multifit = all.sp.fitcontrol)
competitive_ability_all_pairscontrol <- competitive_ability_all_pairscontrol[complete.cases(
competitive_ability_all_pairscontrol),]
competitive_ability_all_pairscontrol
boxplot(competitive_ability_all_pairscontrol$competitive_ability_sp1~competitive_ability_all_pairscontrol$sp1)
###########################################
#### Prepare data: LWA(DROUGHT) ###########
###########################################
# Prepare data ----
datacompdrought<-read.table(here("MCT_analysis", "MCT_data_drought.txt"),header = T, stringsAsFactors = T)
datacomp3drought<-datacompdrought
datacomp3drought$ID<-interaction(datacomp3drought$PLOT,datacomp3drought$COLUMN,datacomp3drought$ROW, datacomp3drought$PLOIDY)
datacomp3drought$ID<-make.unique(as.character(datacomp3drought$ID))
datacomp4drought<-datacomp3drought[,c("Biomasa","ID", "PLOIDY", "Combinacion")]
datacomp5drought<-datacomp4drought %>%
separate(Combinacion, into = c("Combinacion1", "Combinacion2"), sep = "_")
datacomp6drought <- datacomp5drought %>%
mutate(PLOIDY_C = if_else(PLOIDY == Combinacion1, Combinacion2, Combinacion1))
datacomp7drought<-datacomp6drought[,c("ID", "Biomasa","PLOIDY", "PLOIDY_C")]
colnames(datacomp7drought)[2]<-"fitness"
datacomp7drought$PLOIDY<-paste0("P",datacomp7drought$PLOIDY)
datacomp7drought$PLOIDY_C<-paste0("P",datacomp7drought$PLOIDY_C)
datacomp7drought$fitness<-round(datacomp7drought$fitness*1000,0)
datacomp8drought<-datacomp7drought %>%
mutate(`P2x` = as.integer(grepl('\\bP2x\\b', PLOIDY_C)),
`P4x` = as.integer(grepl('4x', PLOIDY_C)),
`P6x` = as.integer(grepl('6x', PLOIDY_C)),
`P12x` = as.integer(grepl('12x', PLOIDY_C)))
datacomp9drought<-datacomp8drought[,-c(3,4)]
datacomp9drought
datalistdrought <- split(datacomp9drought, f = datacomp8drought$PLOIDY)
datalist_cleandrought<-lapply(datalistdrought, na.omit)
#Interaction coefficients and niche differences
datadrought <- datalist_cleandrought
for(i in 1:length(datadrought)){
datadrought[[i]] <- datadrought[[i]][,2:length(datadrought[[i]])]
}
focal_column <- names(datadrought)
model_family <- "RK"
optimization_method <- "bobyqa" # we use a bounded method
alpha_form <- "pairwise"
lambda_cov_form <- "none"
alpha_cov_form <- "none"
initial_values = list(lambda = 10,
alpha_intra = 0.1,
alpha_inter = 0.1)
lower_bounds = list(lambda = 0,
alpha_intra = 0.01,
alpha_inter = 0.01)
upper_bounds = list(lambda = 2000,
alpha_intra = 1,
alpha_inter = 1)
fixed_terms <- NULL
bootstrap_samples <- 0
all.sp.fitdrought <- cxr_pm_multifit(data = datadrought,
model_family = model_family,
focal_column = focal_column,
optimization_method = optimization_method,
alpha_form = alpha_form,
lambda_cov_form = lambda_cov_form,
alpha_cov_form = alpha_cov_form,
initial_values = initial_values,
lower_bounds = lower_bounds,
upper_bounds = upper_bounds,
fixed_terms = fixed_terms,
bootstrap_samples = bootstrap_samples)
all.sp.fitdrought #Intrinsic growth rate and Interaction coefficients
niche_overlap_all_pairsdrought <- niche_overlap(cxr_multifit = all.sp.fitdrought)
niche_overlap_all_pairsdrought #Niche overlap
niche_overlap_all_pairsdrought$MCT_niche_diff <- 1 - niche_overlap_all_pairsdrought$niche_overlap_MCT
niche_overlap_all_pairsdrought$MCT_niche_diff #Niche difference
niche_overlap_all_pairsdrought$SA_niche_diff <- 1 - niche_overlap_all_pairsdrought$niche_overlap_SA
# Average Fitness differences
avg_fitness_diff_all_pairsdrought <- avg_fitness_diff(cxr_multifit = all.sp.fitdrought)
avg_fitness_diff_all_pairsdrought <- avg_fitness_diff_all_pairsdrought[complete.cases(
avg_fitness_diff_all_pairsdrought),] # as with niche overlap, remove NA values
avg_fitness_diff_all_pairsdrought
# Competitive ability
competitive_ability_all_pairsdrought <- competitive_ability(cxr_multifit = all.sp.fitdrought)
competitive_ability_all_pairsdrought <- competitive_ability_all_pairsdrought[complete.cases(
competitive_ability_all_pairsdrought),]
competitive_ability_all_pairsdrought
boxplot(competitive_ability_all_pairsdrought$competitive_ability_sp1~competitive_ability_all_pairsdrought$sp1)
##############################################
############### MCT Graphs ###################
##############################################
############ Coexistence Graphs ############
#HWA (CONTROL). Extract Average Fitness Differences and Niche Differences
AFDcontrol<- avg_fitness_diff_all_pairscontrol[,c(1,2,5)]
AFDcontrol
AFD2control<-subset(AFDcontrol, AFDcontrol$average_fitness_ratio>1)
AFD2control
ND2control<-niche_overlap_all_pairscontrol[,c(1,2,5)]
ND2control
#LWA (DROUGHT). Extract Average Fitness Differences and Niche Differences
AFDdrought<- avg_fitness_diff_all_pairsdrought[,c(1,2,5)]
AFD2drought<-subset(AFDdrought, AFDdrought$average_fitness_ratio>1)
AFD2drought
NDdrought<-niche_overlap_all_pairsdrought[,c(1,2,5)]
NDdrought
NDdrought<-NDdrought[c(3,1,5,2,4,6),] #Change if its necessary
NDdrought
FD<-seq(0, 20, length.out = 501)
rho<-seq(0,10,by = 0.02)
ND<-1-rho
plot(ND2control$MCT_niche_diff,AFD2control$average_fitness_ratio,type="n",
ylim=c(-0.5,16.6),
xlim=c(-4.7,1.5),
cex.lab=1.25, cex.axis=1.3,
xlab="Niche Differences", ylab="Average Fitness Differences")
abline(h=1, lty=2)
abline(v=0, lty=2)
lines(ND,1/rho)
lines(ND,1/-rho+2)
points(ND2control$MCT_niche_diff,AFD2control$average_fitness_ratio,pch=21,bg="lightblue",cex=2)
etiquetas_modificadas <- paste0(gsub("\\D", "", AFD2control$sp1), "x-", gsub("\\D", "", AFD2control$sp2), "x")
text(ND2control$MCT_niche_diff,AFD2control$average_fitness_ratio,
font=2,adj = c(0.5,-1), labels=etiquetas_modificadas,col="darkblue",cex=1.3)
points(NDdrought$MCT_niche_diff,AFD2drought$average_fitness_ratio,pch=21,bg="lightpink",cex=2)
etiquetas_modificadas2 <- paste0(gsub("\\D", "", AFD2drought$sp1), "x-", gsub("\\D", "", AFD2drought$sp2), "x")
text(NDdrought$MCT_niche_diff,AFD2drought$average_fitness_ratio,
font=2,adj = c(0.5,-1),labels=etiquetas_modificadas2,col="#800000",cex=1.3)
##############################################
######## Competitive ability Graph ###########
##############################################
#HWA (CONTROL)
# Assuming your data frame is named 'df'
df <- read.table(header=TRUE, text="
sp1 sp2 competitive_ability_sp1
2 P2x P12x 119.758280
3 P4x P12x 183.201314
4 P6x P12x 11.820549
5 P12x P2x 8.932419
7 P4x P2x 23.469274
8 P6x P2x 100.967627
9 P12x P4x 53.925714
10 P2x P4x 11.975828
12 P6x P4x 13.584266
13 P12x P6x 5.762743
14 P2x P6x 12.371999
15 P4x P6x 42.458210
")
df<-competitive_ability_all_pairscontrol
# Standardize pair order
df <- df %>%
rowwise() %>%
mutate(pair = paste(sort(c(sp1, sp2)), collapse = "-")) %>%
ungroup()
# Find the maximum value for each pair
agg_df <- df %>%
group_by(pair) %>%
summarize(competitive_ability_sp1 = max(competitive_ability_sp1, na.rm = TRUE))
agg_df2c<-df[df$competitive_ability_sp1%in%agg_df$competitive_ability_sp1,-4]
agg_df2c$sp1<-factor(substr(agg_df2c$sp1,2,4), levels=c("2x","4x","6x","12x"))
agg_df2c$sp2<-factor(substr(agg_df2c$sp2,2,4), levels=c("2x","4x","6x","12x"))
# Create a graph object from the data frame
graph <- graph_from_data_frame(agg_df2c, directed=TRUE)
# Calculate arrow widths based on competitive_ability_sp1
arrow_widths <- scales::rescale(agg_df2c$competitive_ability_sp1, to=c(0.8,6))
# Plot the graph using a circular layout
par(mar = c(0, 0, 0, 0))
plot(graph,
vertex.size=30,
vertex.color="#B0C4DE",
# vertex.color=col_p,
vertex.label.cex=2.1,
vertex.label.color="grey25",
edge.width=arrow_widths,
edge.label=round(E(graph)$competitive_ability_sp1, 2), # Add edge labels
edge.label.cex=1.6,
edge.label.color='black',
edge.arrow.size=1.8,
edge.color="black",
edge.label.y = c(0.35, -0.60, -0.60, 0.5,0.5,-0.1), # modifiy label position in y-axis
edge.label.x = c(0.20, -0.70, 0.7,-0.7,0.75,-0.45), # modifiy label position in x-axis
layout=layout_in_circle(graph,order=c("6x","2x","4x","12x")))
#LWA (DROUGHT)
# Assuming your data frame is named 'df'
df <- read.table(header=TRUE, text="
sp1 sp2 competitive_ability_sp1
2 P2x P12x 23.79534
3 P4x P12x 11.34881
4 P6x P12x 69.38080
5 P12x P2x 28.87310
7 P4x P2x 105.52710
8 P6x P2x 268.62062
9 P12x P4x 182.82808
10 P2x P4x 148.94970
12 P6x P4x 149.50079
13 P12x P6x 25.40063
14 P2x P6x 20.27258
15 P4x P6x 15.32260
")
df<-competitive_ability_all_pairsdrought
# Standardize pair order
df <- df %>%
rowwise() %>%
mutate(pair = paste(sort(c(sp1, sp2)), collapse = "-")) %>%
ungroup()
# Find the maximum value for each pair
agg_df <- df %>%
group_by(pair) %>%
summarize(competitive_ability_sp1 = max(competitive_ability_sp1, na.rm = TRUE))
agg_df2s<-df[df$competitive_ability_sp1%in%agg_df$competitive_ability_sp1,-4]
agg_df2s$sp1<-factor(substr(agg_df2s$sp1,2,4), levels=c("2x","4x","6x","12x"))
agg_df2s$sp2<-factor(substr(agg_df2s$sp2,2,4), levels=c("2x","4x","6x","12x"))
# Create a graph object from the data frame
graph <- graph_from_data_frame(agg_df2s, directed=TRUE)
#arrow_head<-rescale(agg_df$competitive_ability_sp1, to=c(0.5,10))
arrow_widths <- scales::rescale(agg_df2s$competitive_ability_sp1, to=c(0.8,6))
# Plot the graph using a circular layout
par(mar = c(0, 0, 0, 0))
plot(graph,
vertex.size=30,
vertex.color="#FFCCCC",
# vertex.color=col_p,
vertex.label.cex=2.1,
vertex.label.color="grey25",
edge.width=arrow_widths,
edge.label=round(E(graph)$competitive_ability_sp1, 2), # Add edge labels
edge.label.cex=1.6,
edge.label.color='black',
edge.arrow.size=1.8,
edge.color="black",
edge.label.y = c(-0.6, -0.45, 0.55, -0.6,0.55,0.10), # modifiy label position in y-axis
edge.label.x = c(0.70, 0.18, 0.7,-0.7,-0.7,0.4), # modifiy label position in x-axis
layout=layout_in_circle(graph,order=c("6x","2x","4x","12x")))
##############################################
# Interaction Coefficients Differences Graph #
##############################################
dataalpha <- read.table(here("MCT_analysis","AlphaDifferences_Graph.txt"), header = TRUE)
dataalpha$ID <- rep(1:16, 2)
df <- dataalpha
#Create 2 subset: HWA (Control) and LWA (Drought)
control_df <- df %>% filter(Treatment == "CONTROL")
drought_df <- df %>% filter(Treatment == "DROUGHT")
# Calculate difference between interaction coefficients (alpha values) Control-Droguth
# Assuming that each ID is unique within each treatment and exists in both
alpha_diff_df <- control_df %>%
dplyr::select(ID, alpha) %>%
left_join(drought_df %>% dplyr::select(ID, alpha), by = "ID", suffix = c("_CONTROL", "_DROUGHT")) %>%
mutate(alpha_diff = alpha_CONTROL - alpha_DROUGHT)
# Create a new variable combining ploidy and ploidy2 in the original dataframe
df <- df %>%
mutate(ploidy_combination = paste(ploidy, ploidy2, sep = "-"))
# We add the calculated alpha difference to the original dataframe, if necessary.
#This assumes you want to keep the original structure and add the calculated difference.
#Note: This will add NA in the rows where the difference is not calculated (for example, where Treatment is neither CONTROL nor DROUGHT).
df <- df %>%
left_join(alpha_diff_df %>% dplyr::select(ID, alpha_diff), by = "ID")
# Order ploidy combination
ploidy_order <- c("12x", "6x", "4x", "2x")
df <- df %>%
mutate(ploidy = factor(ploidy, levels = ploidy_order),
ploidy2 = factor(ploidy2, levels = ploidy_order),
ploidy_combination = factor(ploidy_combination, levels = unique(ploidy_combination[order(ploidy, ploidy2)])))
# Graph
ggplot(df, aes(x = ploidy_combination, y = alpha_diff)) +
geom_segment(aes(xend = ploidy_combination, yend = 0), color = "grey") + # L?neas hacia el eje x
geom_point(aes(color = alpha_diff), size = 6) + # Puntos con color seg?n la diferencia alpha
scale_color_gradient2(low = "#556b2f", mid = "#f5f5dc", high = "#daa520", midpoint = 0, limits = c(-1, 1)) + # Escala de color
coord_flip() + # Invertir ejes para hacerlo horizontal
theme_minimal() + # Tema minimalista
labs(x = "", y = expression(paste(Delta, alpha)),) +
theme(axis.title.x = element_text(size = 22), # Tama?o del t?tulo del eje x
axis.title.y = element_text(size = 1),
axis.text.x = element_text(size = 16, hjust = 0.5), # Ajustar texto del eje x
axis.text.y = element_text(size = 16, angle = 0, hjust = 0.5),
axis.line.x = element_line(color = "grey", size = 0.5),
axis.line.y = element_line(color = "grey", size = 0.5)) + # Ajustar texto del eje y
theme(legend.position = "none") # Quitar la leyenda de colo
##############################################
##### Intrinsic Growth Graph ################
##############################################
datos <- data.frame(
TREATMENT = c("CONTROL", "CONTROL", "CONTROL", "CONTROL", "DROUGHT", "DROUGHT", "DROUGHT", "DROUGHT"),
PLOIDY = factor(c("2X", "4X", "6X", "12X", "2X", "4X", "6X", "12X"), levels = c("2X", "4X", "6X", "12X")),
LAMDA = c(373.145, 201.459, 100.894, 219.540, 445.997, 555.786, 405.316, 388.746)
)
datos$LAMDA <- abs(datos$LAMDA)
# Redefine the PLOIDY factor to the desired order
datos$PLOIDY <- factor(gsub("X", "x", datos$PLOIDY), levels = c("12x", "6x", "4x", "2x"))
# Filter data for HWA (CONTROL) and LWA (DROUGHT)
control <- subset(datos, TREATMENT == "CONTROL")
drought <- subset(datos, TREATMENT == "DROUGHT")
# Combining data to form HWA (CONTROL) and LWA (DROUGHT) pairs
pairs <- merge(control, drought, by = "PLOIDY", suffixes = c("_CONTROL", "_DROUGHT"))
ggplot(datos, aes(x = LAMDA, y = PLOIDY, color = TREATMENT)) +
geom_point(data = control, size = 6) +
geom_point(data = drought, size = 6) +
geom_segment(data = pairs, aes(x = LAMDA_CONTROL, xend = LAMDA_DROUGHT, y = PLOIDY, yend = PLOIDY), color = "black") +
scale_color_manual(values = c("CONTROL" = "blue", "DROUGHT" = "red")) +
theme_minimal() +
labs(x = "Intrinsic growth (mg)", y = "Cytotype", color = "TREATMENT") +
theme(axis.title.x = element_text(size = 20), # Tama?o del t?tulo del eje x
axis.title.y = element_text(size = 20),
axis.text.y = element_text(size = 16, hjust = 1),
axis.text.x = element_text(size = 15, hjust = 0.5),
axis.line.x = element_line(color = "grey", size = 0.5))+
(theme(legend.position = "none")
) +
scale_x_continuous(limits = c(0, 600), breaks = seq(0, 600, by = 100))
#paired test
control <- subset(datos, TREATMENT == "CONTROL")
sequia <- subset(datos, TREATMENT == "DROUGHT")
test_pareado <- t.test(control$LAMDA, sequia$LAMDA, paired = TRUE)
print(test_pareado)