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

knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(data.table)

Introduction

This file is to analyse the data from simulation.

Treatment with "migration001" to "migration004"

Migration period vs. migration frequency migration_period vs. migration_period

Analysis

Read files

Read and merge files. Add model parameters.

# Check file availability
list.files("result/", pattern = "^migration") %>% substr(10, 12) %>%
  as.numeric() %>% factor(1:30) %>% table()

# Treatment name
treatment_name <- paste0("^migration001")

# Read file names
folder <- "result/"
f <- list.files(folder, pattern =  "(^migration01[0-9])|(^migration02[0-4])")

# Read data frame in result
# Add treatment and replicate 
x <- lapply(paste0("result/", f), function (file_name) {
  fread(file_name, head = T) %>%  
    mutate(treatment = gsub("[/.a-z]+", "", file_name) %>% substr(1, 3) %>% as.numeric(), 
           replicate = gsub("[/.a-z]+|(\\d+-)", "", file_name) %>% as.numeric())
})

# Function for reading model attributes----
attribute_read <- function (f) { # file name should be a list
  attribute_df <- rep(list(NA, length(f)))

  for (i in 1:length(f)) {
    # Read the first few lines of comment
    l <- readLines(paste0("result/", f[i]), 1, 1)
    counter <- 1
    while (tail(l, 1) %>% substr(1, 1) == "#") {
      counter <- counter + 1
      l <- readLines(paste0("result/", f[i]), counter)
    }
    l <- l[-counter]

    # Extract attributes
    attribute_extract <- function(x) {
      grep(x, l, value = TRUE) %>% strsplit(" ") %>% unlist() %>% `[`(length(.)) %>% as.numeric()
    }

    # Attribute name
    attribute_name <- c("P", "S_pool", "I", "S", "J", "mi",
                        "p", "cross_feeding", "function_group", 
                        "same_initial_community", "migration_size", 
                        "migration_period", "migration_event_times",
                        "same_migrant_community",
                        "time_step", "time_stabilization", "replicate")

    attribute_df[[i]] <- 
      strsplit(l, " ") %>% lapply(function(x) tail(x, 1)) %>% 
      unlist() %>% `[`(-8) %>%
      matrix(nrow = 1) %>% as.tibble() %>%
      setNames(attribute_name) %>%
      mutate(treatment = gsub("[/.a-z]+", "", f[i]) %>% substr(1, 3) %>% as.numeric())
  }


  if (length(f) >= 2) {
    attribute_df <- rbindlist(attribute_df)
  } else if (length(f) == 1) attribute_df <- attribute_df[[1]]

  # Return
  return(attribute_df)
}

model_attribute <- attribute_read(f)
model_attribute1 <- 
  model_attribute %>%
  mutate(treatment = as.numeric(treatment),
         replicate = as.numeric(treatment)) %>%
  select(J, migration_size, migration_period, migration_event_times, 
         treatment, replicate)

# Merge replicates----
x <- rbindlist(x)
x <- left_join(x, model_attribute1, 
               by = c("treatment", "replicate"))

Biodiversity indexes

Read and merge data. This chunk only deals with one treatment (20 replicates). The loop for all treatments is in another R script call migration_analysis.R

# Calculate relative abundance
x <- 
  x %>%
  group_by(replicate, time, type) %>%
  mutate(biomass_sum = sum(value),
         biomass_proportion = value / biomass_sum)

Plot alpha diveristy of one treatement, colored by replicates

x %>%
  filter(biomass_proportion > 0, type == "X") %>%
  summarize(alpha = sum(-biomass_proportion * log(biomass_proportion))) %>%
  group_by(time, type) %>% 
  mutate(replicate = factor(replicate, 1:20)) %>% 
  ggplot(aes(x = time, y = alpha, color = replicate)) + 
  geom_line() +
  geom_vline(xintercept = seq(model_attribute$migration_period,
                              model_attribute$migration_event_times * model_attribute$migration_period, 
                              model_attribute$migration_period)) +
  theme_bw()

Alpha diversity over time

x_summary <- 
  x %>%
  filter(biomass_proportion > 0) %>%
  summarize(alpha = sum(-biomass_proportion * log(biomass_proportion))) %>%
  group_by(time, type) %>%
  summarize(alpha_mean = mean(alpha),
            alpha_sd = sd(alpha))

pdf(paste0("plot/alpha_", treatment, ".pdf"), width = 10, height = 8)
x_summary %>%
  filter(type == "X") %>%
  ggplot(aes(x = time, y = alpha_mean)) +
  geom_line(lwd = .1) +
  geom_segment(aes(x = time, xend = time,
                   y = alpha_mean - alpha_sd,
                   yend = alpha_mean + alpha_sd)) + 
  geom_vline(xintercept = seq(model_attribute$migration_period,
                              model_attribute$migration_event_times * model_attribute$migration_period, 
                              model_attribute$migration_period)) +
       theme_bw() +
      ggtitle(paste0("Flow rate = ", model_attribute$J,
                     ", migration event times = ", model_attribute$migration_event_times,
                     ", migration size = ", model_attribute$migration_size,
                     ", migration period = ", model_attribute$migration_period,
                     ", cross-feeding = ", model_attribute$cross_feeding))
dev.off()

Beta diversity over time

Bray-Curtis dissimilarity. Use average pairwise BC dissimilarity Calculate by a function from package vegan.

# Split data frame by time
x_time <- filter(x, type == "X")
x_time <- split(x_time, x_time$time)

x_BC <- 
  x_time %>%
  sapply(function (df) {
    BC <- 
      dcast(df, formula = time + replicate ~ variable, 
            fill = 0) %>%
      select(starts_with("X"))
    BC[BC < 0] <- 0

    BC <- vegdist(BC, method = "bray")
    return(c(BC_mean = mean(BC),
             BC_sd = sd(BC)))
  })

x_BC2 <- 
  x_BC %>% t() %>% 
  as.data.frame() %>%
  mutate(time = unique(x$time))

Plot one treatment

pdf(paste0("plot/beta_", treatment, ".pdf"), width = 10, height = 8)
x_BC2 %>%
  ggplot(aes(x = time, y = BC_mean)) +
  geom_line(lwd = .1) +
  geom_segment(aes(x = time, xend = time,
                   y = BC_mean - BC_sd,
                   yend = BC_mean + BC_sd)) + 
  geom_vline(xintercept = seq(model_attribute$migration_period,
                              model_attribute$migration_event_times * model_attribute$migration_period, 
                              model_attribute$migration_period)) +
  theme_bw() +
  ggtitle(paste0("Flow rate = ", model_attribute$J,
                 ", migration event times = ", model_attribute$migration_event_times,
                 ", migration size = ", model_attribute$migration_size,
                 ", migration period = ", model_attribute$migration_period,
                 ", cross-feeding = ", model_attribute$cross_feeding))
dev.off()

Save Beta diversity and alpha diversity data for each treatment. Use this code chunk in the loop of script migration_analysis.R.

# Merge alpha and beta diversity data
x_summary2 <- x_summary %>% filter(type == "X")
x_diversity <- inner_join(x_summary2, x_BC2) %>%
  mutate(treatment = substr(treatment, 10, 12) %>% as.numeric,
         J = model_attribute$J,
         migration_event_times = model_attribute$migration_event_times,
         migration_size = model_attribute$migration_size,
         migration_period = model_attribute$migration_period,
         cross_feeding = model_attribute$cross_feeding)

#
file_name <- file("result/diversity_time.txt")
fwrite(x_diversity, file = paste0("result/diversity_time.txt"), 
       append = T, col.names = ifelse((substr(treatment, 10, 12) %>% as.numeric) == 1,
                                      TRUE, FALSE))
close(file_name)

Diversity in all treatments

Merge file

# Merge files
div_time <- 
  list.files("result/", pattern = "diversity_time_migration") %>%
  lapply(function(x) {
    fread(paste0("result/", x))
  }) %>%
  rbindlist()

fwrite(div_time, "result/diversity_time.txt")

# Read merged file
div_time <- fread("result/diversity_time.txt")

pdf("plot/beta_time.pdf", width = 10, height = 8)
div_time %>%
  mutate(treatment = factor(treatment, 1:20)) %>%
  filter(treatment %in% 1:20) %>%
  ggplot(aes(x = time, y = BC_mean)) +
  geom_line(lwd = .3) +
  facet_grid(migration_size ~ J, labeller = label_both) +
  geom_vline(xintercept = 10 * (1:20), color = "red", lwd = .05) + 
  theme_bw() +
  ylim(0, 1) +
  labs(y = "Bray-Curtis distance") +
  ggtitle(paste0("Species pool = ", 10000))
dev.off()

Compare diversitieis at the beginning and the last time point

pdf("plot/beta_J_migration.pdf")#, width = 6, height = 10)
div_time %>%
  filter(treatment %in% 1:20) %>%
  filter(time %in% c(min(time), max(time))) %>%
  mutate(time = factor(time)) %>%
  ggplot(aes(x = migration_size, y = BC_mean, color = time)) +
  geom_point() + 
  geom_segment(aes(x = migration_size, xend = migration_size,
                   y = BC_mean - BC_sd,
                   yend = BC_mean + BC_sd)) +
  facet_grid(J ~ ., labeller = label_both) +
  theme_bw() +
  ylim(.7, 1) +
  labs(y = "Bray-Curtis Dissimilarity")

div_time %>%
  filter(treatment %in% 1:20) %>%
  filter(time %in% c(min(time), max(time))) %>%
  mutate(time = factor(time)) %>%
  ggplot(aes(x = J, y = BC_mean, color = time)) +
  geom_point() + 
  #geom_segment(aes(x = J, xend = J,
  #                 y = BC_mean - BC_sd,
  #                 yend = BC_mean + BC_sd)) +
  facet_grid(migration_size ~ ., labeller = label_both) +
  theme_bw() +
  ylim(.7, 1) +
  labs(y = "Bray-Curtis Dissimilarity")
dev.off()

Richness

# Read merged file
div_time <- fread("result/diversity_time.txt")

pdf("plot/richness_time.pdf", width = 10, height = 8)
div_time %>%
  mutate(treatment = factor(treatment, 1:20)) %>%
  filter(treatment %in% 1:20) %>%
  ggplot(aes(x = time, y = richness_mean)) +
  geom_line(lwd = .3) +
  facet_grid(migration_size ~ J, labeller = label_both) +
  geom_vline(xintercept = 10 * (1:20), color = "red", lwd = .05) + 
  theme_bw() +
  labs(y = "Species richness") +
  ggtitle(paste0("Species pool = ", 10000))
dev.off()

Sign of increased/decreased BC distance at migration

div_time_mig <- 
div_time %>%
  mutate(treatment = factor(treatment, 1:20)) %>%
  filter(treatment %in% 1:20) %>%
  filter(time %in% c(seq(9, 199, 10), seq(10, 200, 10))) 
#%>% mutate(mig = rep(rep(1:20, each = 2), times = 10))

div_time_mig <- 
  div_time_sub %>%
  filter(time %in% seq(9, 199, 10)) %>%
  select(treatment, J, migration_event_times, 
         migration_size, migration_period, cross_feeding) %>%
  mutate(BC_sign = div_time_mig$BC_mean[seq(2, 400, 2)] - div_time_mig$BC_mean[seq(1, 399, 2)], mig = rep(1:20, 10)) %>%
  mutate(BC_sign = ifelse(BC_sign > 0, 1, 0)) 

pdf("plot/BC_sign.pdf", width = 10, height = 8)
div_time_mig %>%
  ggplot(aes(x= mig, y = BC_sign)) +
  geom_line() +
  facet_grid(migration_size ~ J, labeller = label_both) +
  theme_bw() +
  labs(x = "migration event", y = "Sign of BC increase")
dev.off()

Biomass

# Read merged file
div_time <- fread("result/diversity_time.txt")

pdf("plot/biomass_time.pdf", width = 10, height = 8)
div_time %>%
  mutate(treatment = factor(treatment, 1:20)) %>%
  filter(treatment %in% 1:20) %>%
  ggplot(aes(x = time, y = biomass_sum_mean)) +
  geom_line(lwd = .3) +
  facet_grid(migration_size ~ J, labeller = label_both) +
  geom_vline(xintercept = 10 * (1:20), color = "red", lwd = .05) + 
  theme_bw() +
  labs(y = "Species richness") +
  ggtitle(paste0("Species pool = ", 10000))
dev.off()

Plot

# Plot 
if(FALSE){
  result %>%
    as.data.frame() %>%
    melt(id.vars = "time") %>%
    mutate(type = substr(variable, 1, 1)) %>%
    ggplot(aes(time, value, color = variable)) +
    geom_line() +
    facet_grid(type ~ .) +
    xlim(0, time_limit) +
    labs(x = "time", y = "biomass") +
    theme_bw() +
    guides(color = FALSE) # Remove legends
}

Benchmark

mbm <- microbenchmark(
  a = {
    sp <- community_generate(S_pool = 10000, S = S, r = 1)
    sp
    result <- CR_model(sp, P = 5,
                       P_biomass = setNames(rep(1, 5),
                                            paste0("R", 1:5)),
                       S_biomass = setNames(rep(1, length(sp)),
                                            paste0("X", sprintf("%05d", sp))),
                       time_limit = 10, time_step = 1) %>%
      as.data.frame()
  },
  b = {
    sp <- community_generate(S_pool = 10000, S = S, r = 1)
    sp <- c(sp, 3, 4, 5)
    result <- CR_model(sp, P = 5,
                       P_biomass = setNames(rep(1, 5),
                                            paste0("R", 1:5)),
                       S_biomass = setNames(c(rep(1, length(sp)-3), rep(0, 3)),
                                            paste0("X", sprintf("%05d", sp))),
                       time_limit = 10, time_step = 1) %>%
      as.data.frame()

  },
  times = 100
)
autoplot(mbm)