# !diagnostics off
library(tidyverse)
library(deSolve)
library(data.table)
library(micro.crm)
library(vegan)
rm(list=ls())
# Augument passed from command line
args = commandArgs(trailingOnly=TRUE)
treatment_name = args[1]
# Read files from the community ----
parameter_load(parameter_text = paste0("parameter_script/", treatment_name, ".txt"))
## Community
result <- result_read(dir = "simulation_result/", pattern = treatment_name, threshold = threshold)
## Isolate
#result_isolate <- result_isolate_read(dir = "simulation_result/isolates/", pattern = treatment_name, threshold = threshold)
## Pair
#result_pair <- result_pair_read(dir = "simulation_result/pairs/", pattern = treatment_name, threshold = threshold)
# Plot ----
p_title <- paste0("u_main=", group_mu, ", d_general=",d_general, ", family=", function_group, ", same_stoi=", same_stoi)
## Community ----
### Dynamics
pdf(paste0("figure/dynamics_", treatment_name, ".pdf"))
for (r in 1:replicate) {
p0 <- result %>%
filter(replicate == r) %>%
ggplot(aes(x = time, y = value, color = variable)) +
geom_line() +
facet_grid(type~.) +
theme_bw() +
guides(color = F) +
labs(x = "time") +
ggtitle("replicate=", r)
print(p0)
}
dev.off()
### Final composition
communities <- filter(result, time == max(time), type == "X") %>% select(-time, -type)
p1 <- communities %>%
arrange(value) %>%
ggplot(aes(x = replicate, y = value, fill = function_group)) +
geom_bar(stat="identity", position = 'fill', color = "black") +
facet_grid(.~treatment) +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
scale_fill_discrete() +
theme_bw() +
guides(fill = guide_legend(title = "Family"), color = F) +
labs(x = "replicate", y = "relative abundancee") +
ggtitle(p_title)
pdf(paste0("figure/composition_", treatment_name,".pdf")); p1; dev.off()
### Beta diversity
#beta_diversity(communities)
if (FALSE) {
## Pair----
if (length(result_pair$value) != 0) {
# All pairwise interactions
pairs <- result_pair %>% filter(time == max(time), type == "X") %>%
dplyr::select(-time, -type, -mixing_ratio, -variable, -function_group) %>%
group_by(treatment, replicate, pair) %>%
mutate(isolate = rep(c("abundance1", "abundance2"))) %>%
spread(isolate, value)
# Isolate ID
pairs$isolate1 <- sapply(strsplit(pairs$pair, "_"), function(x)`[`(x, 1))
pairs$isolate2 <- sapply(strsplit(pairs$pair, "_"), function(x)`[`(x, 2))
# Functional groups
pairs$family1 <- cut(as.numeric(substr(pairs$isolate1, 2, 6)), breaks = seq(0, S_pool, S_pool/function_group), labels = 1:function_group)
pairs$family2 <- cut(as.numeric(substr(pairs$isolate2, 2, 6)), breaks = seq(0, S_pool, S_pool/function_group), labels = 1:function_group)
# Competition outcome
pairs <- pairs %>% mutate(interaction = ifelse((abundance1!=0 & abundance2==0) | (abundance1==0 & abundance2!=0), "exclusion", ifelse(all(c(abundance1, abundance2) != 0), "coexistence", "non-growth"))) %>%
mutate(interaction = factor(interaction, levels = c("coexistence", "exclusion", "non-growth")))
# From-to for network
pairs <- pairs %>% ungroup %>%
mutate(from = ifelse(abundance1 > abundance2, isolate1, isolate2),
to = ifelse(abundance1 > abundance2, isolate2, isolate1))
# Fraction of coexistence pairs
pairs1 <- filter(pairs, family1 %in% c(1,2), family2 %in% c(1,2))
pairs1$two_families <- factor(ifelse((pairs1$family1 != pairs1$family2), "F-N", ifelse(pairs1$family1 == 1, "F-F", "N-N")), levels = c("F-F", "F-N", "N-N"))
## Plot
p2 <- pairs1 %>%
select(interaction, two_families) %>%
table() %>% data.frame() %>%
ggplot(aes(x = two_families, y = Freq, fill = interaction)) +
geom_col(position = "dodge", color = "black") +
theme_bw() +
theme(legend.position = "top", legend.title = element_blank()) +
labs(x = "") +
ggtitle(p_title)
pdf(paste0("figure/pair_", treatment_name,".pdf")); p2; dev.off()
}
## Isolate ----
if (length(result_isolate$value) != 0) isolates <- filter(result_isolate, time == max(time), type == "X")
}
## Save final time step ---
write.table(communities, paste0("simulation_result/equilibrium/", treatment_name, "_communities.txt"), row.names = F)
if (FALSE) {
write.table(pairs, paste0("simulation_result/equilibrium/", treatment_name, "_pairs.txt"), row.names = F)
write.table(isolates, paste0("simulation_result/equilibrium/", treatment_name, "_isolates.txt"), row.names = F)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.