# !diagnostics off knitr::opts_chunk$set(cache = T, echo = F) library(micro.crm) library(igraph) library(tidyverse) library(data.table)
This script includes the following steps:
Isolation. This step reads the result from self-assebled communities and subset the communities at final time step.
Pairwise competition within communities. For each self-assembled communities, conduct pairwise competition.
Run the chunk below in bash.
```{bash eval = F}
Rscript self-assembly.R self-assembly001
Rscript self-assembly-isolate.R self-assembly001
Rscript self-assembly-pair.R self-assembly001 0.5 0.5
# Read files from the community ```r rm(list=ls()) treatment_name <- "self-assembly004" 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)
The analysis below takes only one treatment self-assembly004
for example.
# Richness communities %>% dplyr::group_by(treatment, replicate) %>% summarize(richness = n())
communities <- filter(result, time == max(time), type == "X") %>% dplyr::select(-time, -type) # Final composition p_title <- paste0("d_general=",d_general, ", group_mu=", group_mu, ", function_group=", function_group) 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(communities)
# 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 %>% dplyr::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()
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) 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)
df_obv <- as.data.frame(matrix(NA, 20, 7)) %>% setNames(paste0("motif", 1:7)) # Motif count for (i in 1:20) { df_obv[i,] <- motif_count(network_make(nodes = filter(communities, replicate == i),edges = filter(pairs, replicate == i))) } df_obv_sum <- df_obv %>% mutate(Replicate = 1:20) %>% gather("Motif", "Count", 1:7) %>% group_by(Motif) %>% summarize(Count = sum(Count)) # Sum across all communities df_obv_sum
net <- network_make(nodes = filter(communities, replicate == 9), edges = filter(pairs, replicate == 9)) network_plot(net)
net <- network_make(nodes = filter(communities, replicate == 9), edges = filter(pairs, replicate == 9)) net_random <- network_randomize(net, keep_direction = F) adjacency_make(net, isolate_rank = elo_rank(net)$Player) adjacency_make(net_random, isolate_rank = elo_rank(net_random)$Player) dist_to_diag(adjacency_make(net_random, isolate_rank = elo_rank(net_random)$Player)) df_distance_list[[2]]
# Ranking elo_ranking <- elo_rank(net) # Heatmap adjacency_make(net, elo_ranking$Player) %>% heatmap_network_plot()
Number of coexistence to diagonal
m <- adjacency_make(net, elo_ranking$Player) dist_to_diag(m)
df_motif_list <- rep(list(as.data.frame(matrix(NA, b, 7))), 20) df_distance_list <- rep(list(NA), 20) b = 100 for (i in 1:20) { # 20 replicates i = 1 net <- network_make(nodes = filter(communities, replicate == i), edges = filter(pairs, replicate == i)) df_distance_list <- rep(list(NA), b) for (j in 1:b) { net_random <- network_randomize(net, keep_direction = T) df_motif_list[[i]][j,] <- motif_count(net_random) df_distance_list[[j]] <- dist_to_diag(adjacency_make(net_random, isolate_rank = elo_rank(net_random)$Player)) } print(i) }
df_random <- df_list %>% do.call("rbind", .) %>% setNames(paste0("motif", 1:7)) %>% mutate(Replicate = rep(1:20, each = b), Shuffle = rep(1:b, 20)) %>% gather("Motif", "Count", 1:7) %>% group_by(Motif, Shuffle) %>% summarize(Count = sum(Count)) %>% # group_by(Motif) %>% summarize(CountMean = mean(Count), CountSd = sd(Count)) %>% ungroup() {.} df_random_sum <- df_random %>% group_by(Motif) %>% summarize(CountMean = mean(Count), CountSd = sd(Count)) df_random_percentile <- df_random %>% arrange(Motif, Count) %>% select(-Shuffle) %>% split.data.frame(f = df_random$Motif) %>% lapply(function(x) x[c(0.95*b, 0.05*b),]) %>% do.call("rbind", .) %>% mutate(Percentile = c("p95", "p5")) %>% spread(Percentile, Count) %>%
left_join(df_obv_sum, df_random_percentile, by = "Motif") %>% left_join(df_random_sum, by = "Motif") %>% ggplot() + geom_point(aes(x = Motif, y = CountMean, shape = "Random network"), position = position_dodge(width = 1)) + geom_errorbar(aes(x = Motif, ymax = p95, ymin = p5), position = position_dodge(width = 1), width = .2) + geom_point(aes(x = Motif, y = Count, shape = "Observation"), position = position_dodge(width = 1), size = 1) + # scale_y_continuous(limits = c(0, 70)) + theme_bw() + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank(), legend.position = "top") + facet_grid(.~Motif, scale = "free") + scale_shape_manual(name = NULL, values = c(16, 3)) + xlab("Motif") + ylab("Count") + ggtitle("")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.