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

# !diagnostics off
knitr::opts_chunk$set(cache = T, echo = F)
library(micro.crm)
library(tidyverse)
library(data.table)
library(vegan)
library(igraph)

Simulation

Run this chunk in bash. The parameter for simulation is save in the folder parameter_script/

```{bash eval = F}

Self-assembly

Rscript self-assembly.R self-assembly004

Isolate

Rscript self-assembly_isolate.R self-assembly004

Pairwise competition

Rscript self-assembly_pair.R self-assembly004 0.5 0.5

# Read simulated result

```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)

Plot

The analysis below takes only one treatment self-assembly004 for example.

Community

Richness

communities <- filter(result, time == max(time), type == "X") %>% dplyr::select(-time, -type)

# Richness
communities %>% 
  dplyr::group_by(treatment, replicate) %>%
  summarize(richness = n())

Community composition at equilibrium in bar chart

# 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)
p1
#pdf(paste0("figure/composition_", treatment_name,".pdf")); p1; dev.off()

Beta diversity

beta_diversity(communities)

Pair

Pairwise competitive outcome

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

p2
#pdf(paste0("figure/pair_", treatment_name,".pdf")); p2; dev.off()

Isolate

isolates <- filter(result_isolate, time == max(time), type == "X")

Write final time step data

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

Network analysis

Network analysis includes the following:

1. Motif count
2. Network randomization
3. Fraction of coexistence as a function of distance to the nearest diagonal tile

Motif count

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)

Elo ranking and coexistence at diagonal

net <- network_make(nodes = filter(communities, replicate == 9), edges = filter(pairs, replicate == 9))
network_plot(net)

# Ranking
elo_ranking <- elo_rank(net)

# Heatmap; to order the isolates by ranking, elo_ranking is required
adjacency_make(net, elo_ranking$Player) %>% heatmap_network_plot()

Number of coexistence as a function of distance to the nearest diagonal tile.

m <- adjacency_make(net, elo_ranking$Player)
dist_to_diag(m)

Randomization

b = 100
df_motif_list <- rep(list(as.data.frame(matrix(NA, b, 7))), 20)
df_distance <- rep(list(NA), 20)

for (i in 1:20) { # 20 replicates
  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 = F)
    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)) %>%
      mutate(shuffle = j)
  }
  df_distance[[i]] <- do.call("rbind", df_distance_list) %>% mutate(replicate = i) 
  df_motif_list[[i]] <- mutate(df_motif_list[[i]], shuffle = 1:b)
  print(i)
}

do.call("rbind", df_distance) %>% filter(replicate==2)

do.call("rbind", df_distance) %>%
  mutate(count = as.character(count)) %>%
  group_by(replicate, shuffle) %>%
  ggplot(aes(x = count, y = fracDraw)) +
  geom_violin() +
  facet_wrap(.~replicate) +
  theme_bw()

df_random <- df_motif_list %>% 
  do.call("rbind", .) %>%
  setNames(c(paste0("motif", 1:7), "shuffle")) %>%
  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) %>% 
#  mutate(CountFrac = Count/sum(Count)) %>%
  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("")

Chang-Yu-Chang/MigrationCommunity documentation built on Aug. 13, 2019, 9:41 p.m.