output/script/migration-parent_analysis.R

In Chang-Yu-Chang/MigrationCommunity: Microbial Consumer-Resource Model

# !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)
}