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(igraph)

Example

net <- null_network_make(n = 5, frac_of_coext = 1)
network_plot(net)
motif_count(net)

n = 5

```{bash eval=F} Rscript motif_null.R 5 1000

```r
df <- fread("simulation_result/motif_null_n5_1000.txt")
n=5
df <- mutate(df, countFrac = count/choose(n, 3))

pdf(file = paste0("figure/motif_null_n5_1000.pdf"), height=5, width=15)
# Boxplot plot
df %>%
  ggplot(aes(x = fracCoext, group = fracCoext, y = countFrac, fill = motif)) +
  geom_boxplot() +
  guides(fill = F) +
  facet_grid(.~motif) +
  theme_bw() +
  labs(x="fraction of coexistence", y="fraction of motif") +
  ggtitle(paste0("number of node =", n))

# Point
df %>%
  ggplot(aes(x = fracCoext, y = countFrac, color = motif)) +
  geom_point() +
  geom_smooth(, method = "lm") +
  guides(color = F) +
  facet_grid(.~motif) +
  theme_bw() +
  labs(x="fraction of coexistence", y="fraction of motif") +
  ggtitle(paste0("number of node =", n))

dev.off()

pdf("figure/motif_null_mean_n5_1000.pdf")

# Average; shade
df %>%
#  filter(motif %in% paste0("motif", c(1,3,5,7))) %>%
  group_by(fracCoext, motif)  %>%
  summarize(countMean = mean(countFrac), countSd = sd(countFrac)) %>%
  ggplot(aes(x = fracCoext, y = countMean, fill = motif)) +
  geom_area() +
  scale_fill_manual(values = grey(seq(1,0,length.out = 7))) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_bw() +
  labs(x="fraction of coexistence links", y = "faction of motif") +
  ggtitle(paste0("number of node =", n))

# Average
df %>%
  group_by(fracCoext, motif)  %>%
  summarize(countMean = mean(countFrac), countSd = sd(countFrac)) %>%
  ggplot(aes(x = fracCoext, y = countMean, color = motif)) +
  geom_point() + geom_line() +
  geom_segment(aes(x=fracCoext, xend=fracCoext, y=countMean+2*countSd, yend=countMean-2*countSd)) +
  theme_bw() +
  labs(x="fraction of coexistence links", y = "faction of motif") +
  ggtitle(paste0("number of node =", n))


dev.off()

n = 15

```{bash eval=F} Rscript motif_null.R 15 1000

```r
df <- fread("simulation_result/motif_null_n15_1000.txt")
n = 15
df <- mutate(df, countFrac = count/choose(n, 3))

pdf(file = paste0("figure/motif_null_n15_1000.pdf"), height=5, width=15)

# Boxplot plot
df %>%
  ggplot(aes(x = fracCoext, group = fracCoext, y = countFrac, fill = motif)) +
  geom_boxplot() +
  guides(fill = F) +
  facet_grid(.~motif) +
  theme_bw() +
  labs(x="fraction of coexistence", y="fraction of motif") +
  ggtitle(paste0("number of node =", n))

# Point
df %>%
  ggplot(aes(x = fracCoext, y = countFrac, color = motif)) +
  geom_point() +
  geom_smooth(, method = "lm") +
  guides(color = F) +
  facet_grid(.~motif) +
  theme_bw() +
  labs(x="fraction of coexistence", y="fraction of motif") +
  ggtitle(paste0("number of node =", n))
dev.off()

pdf("figure/motif_null_mean_n15_1000.pdf")

# Average; shade
df %>%
#  filter(motif %in% paste0("motif", c(1,3,5,7))) %>%
  group_by(fracCoext, motif)  %>%
  summarize(countMean = mean(countFrac), countSd = sd(countFrac)) %>%
  ggplot(aes(x = fracCoext, y = countMean, fill = motif)) +
  geom_area() +
  scale_fill_manual(values = grey(seq(1,0,length.out = 7))) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_bw() +
  theme(legend.position = "top") +
#  guides(fill = F) +
  labs(x="fraction of coexistence links", y = "faction of motif") +
  ggtitle(paste0("number of node =", n))

# Average
df %>%
  group_by(fracCoext, motif)  %>%
  summarize(countMean = mean(countFrac), countSd = sd(countFrac)) %>%
  ggplot(aes(x = fracCoext, y = countMean, color = motif)) +
  geom_point() + geom_line() +
  geom_segment(aes(x=fracCoext, xend=fracCoext, y=countMean+2*countSd, yend=countMean-2*countSd)) +
  theme_bw() +
  theme(legend.position = "top") +
#  guides(fill = F) +
  labs(x="fraction of coexistence links", y = "faction of motif") +
  ggtitle(paste0("number of node =", n))


dev.off()

Enrichment in motifs 1, 4, 6, 7 would suggest that pairwise interactions are preferred to maintain diversity, wherease enrichment in motifs would suggest HOIs.

pdf("figure/motif_null_mean_n15_1000_HOI.pdf")
df %>%
  mutate(suggInteraction = ifelse(motif %in% paste0("motif", c(1,4,6,7)), "pairwise interaction", "HOI"),
         replicate = rep(1:1000, 7*11)) %>%
  group_by(fracCoext, suggInteraction, replicate) %>%
  summarize(count = sum(count)) %>%
  group_by(fracCoext, suggInteraction) %>%
  mutate(countFrac = count/choose(n, 3)) %>%
  summarize(countMean = mean(countFrac), countSd = sd(countFrac)) %>%
  ggplot(aes(x = fracCoext, y = countMean, fill = suggInteraction)) +
  geom_area(color = "black") + 
  geom_vline(xintercept = 0.3, lty = 2) +
  scale_fill_manual(values = c("white", grey(.4))) +
  theme_bw() +
  theme(legend.title = element_blank(), legend.position = "top") +
  #guides(fill = F) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  labs(x="fraction of coexistence links", y = "faction of motif") +
  ggtitle(paste0("number of node =", n))
dev.off()