Chapter 4 Community composition

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>%
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>% 
  filter(!is.na(day)) %>% 
  group_by(genome,treatment,day) %>% 
  summarise(count=mean(count)) %>% 
  left_join(gift_pcoa_vectors %>% rownames_to_column(var="genome"), by="genome") %>%
  left_join(genome_metadata, by="genome") %>%
  ggplot() +
      scale_color_manual(values=order_colors)+

      #genome positions
      geom_point(aes(x=Axis.1,y=Axis.2, color=order, size=count), 
                 alpha=0.8, shape=16) +
      #scale_color_manual(values=phylum_colors) +
      scale_size_continuous(range = c(0,10)) +
      #loading positions

     #Primary and secondary scale adjustments
     scale_x_continuous(name = paste0("PCoA1 (",round(gift_pcoa_rel_eigen[1]*100, digits = 2), " %)"),
                      sec.axis = sec_axis(~ . / scale, name = "Loadings on PCoA1")
            ) +
     scale_y_continuous(name = paste0("PCoA2 (",round(gift_pcoa_rel_eigen[2]*100, digits = 2), " %)"),
                      sec.axis = sec_axis(~ . / scale, name = "Loadings on PCoA2")
            ) +
    facet_grid(treatment ~ day) +
    theme_linedraw() + 
    theme(legend.position = "none")

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>%
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>% 
  filter(!is.na(day)) %>% 
  group_by(genome,treatment,day) %>% 
  summarise(count=mean(count)) %>% 
  left_join(gift_pcoa_vectors %>% rownames_to_column(var="genome"), by="genome") %>%
  left_join(genome_metadata, by="genome") %>%
  ggplot(aes(x=Axis.1, y=count, group=treatment, color=treatment)) +
        scale_color_manual(values=treatment_colors)+
        geom_smooth(method = lm, formula = y ~ splines::bs(x, 3), se = FALSE) +
        facet_grid(. ~ day) +
        theme_minimal() + 
        theme(legend.position = "none")

genome_gifts %>%
    to.elements(., GIFT_db) %>% 
    to.functions(., GIFT_db) %>% 
    to.community(., genome_counts_filt %>%
                    mutate_at(vars(-genome),~./sum(.)) %>%
                   column_to_rownames(var="genome"), GIFT_db) %>% 
    as.data.frame() %>% 
    rownames_to_column(var="sample") %>% 
    pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
    left_join(sample_metadata, by = join_by(sample == sample)) %>% 
    group_by(trait,treatment,day) %>% 
    summarise(value=mean(value)) %>% 
    mutate(trait=factor(trait)) %>% 
    filter(treatment != "TG0") %>% 
    ggplot(aes(x=fct_reorder(trait, desc(value)), y=value, fill=trait)) +
        geom_bar(stat="identity", width=1) +
    facet_grid(treatment ~ day) +
    theme_minimal() + 
    theme(legend.position = "none")

genome_gifts %>%
    to.elements(., GIFT_db) %>% 
    to.functions(., GIFT_db) %>% 
    to.community(., genome_counts_filt %>%
                    mutate_at(vars(-genome),~./sum(.)) %>%
                   column_to_rownames(var="genome"), GIFT_db) %>% 
    as.data.frame() %>% 
    rownames_to_column(var="sample") %>% 
    pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
    left_join(sample_metadata, by = join_by(sample == sample)) %>% 
    group_by(sample,treatment,day) %>% 
    summarise(value=mean(value)) %>% 
    filter(treatment != "TG0") %>% 
    ggplot(aes(x=day, y=value, group=day, fill=treatment, color=treatment)) +
        geom_boxplot() +
        geom_jitter() +
        scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
    facet_grid(treatment ~ .) +
    theme_linedraw() + 
    theme(legend.position = "none")

load("data/data.Rdata")

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

# Remove genomes with no hits
dist_filt <- dist[!rownames(dist) %in% c("GEXTRA:bin_000002", "GEXTRA:bin_000006"), 
                               !colnames(dist) %in% c("GEXTRA:bin_000002", "GEXTRA:bin_000006")]

functional <- genome_counts_filt %>%
  arrange(match(genome,rownames(dist))) %>%
  filter(!genome %in% c("GEXTRA:bin_000002", "GEXTRA:bin_000006")) %>% 
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist_filt) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

#Richness
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="richness") %>%
      ggplot(aes(y = value, x = day, group=day, color=treatment, fill=treatment)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      facet_grid(treatment ~ ., scales = "fixed") +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

#Neutral
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="neutral") %>%
      ggplot(aes(y = value, x = day, group=day, color=treatment, fill=treatment)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      facet_grid(treatment ~ ., scales = "fixed") +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      ) +
    theme_linedraw() + 
    theme(legend.position = "none")

#Phylogenetic
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="phylogenetic") %>%
      ggplot(aes(y = value, x = day, group=day, color=treatment, fill=treatment)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      facet_grid(treatment ~ ., scales = "fixed") +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

#Functional
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="functional") %>%
      ggplot(aes(y = value, x = treatment, group=treatment, color=treatment, fill=treatment)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      facet_wrap(. ~ day, scales = "fixed", ncol=5) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

#Richness
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="richness") %>%
      ggplot(aes(y = value, x = day, group=treatment, color=treatment, fill=treatment)) +
      geom_point() +
      geom_smooth(method="loess",se=TRUE)+
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      coord_cartesian(xlim = c(7, 35)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="neutral") %>%
      ggplot(aes(y = value, x = day, group=treatment, color=treatment, fill=treatment)) +
      geom_point() +
      geom_smooth(method="loess",se=TRUE)+
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      coord_cartesian(xlim = c(7, 35)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="phylogenetic") %>%
      ggplot(aes(y = value, x = day, group=treatment, color=treatment, fill=treatment)) +
      geom_point() +
      geom_smooth(method="loess",se=TRUE)+
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      coord_cartesian(xlim = c(7, 35)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="functional") %>%
      ggplot(aes(y = value, x = day, group=treatment, color=treatment, fill=treatment)) +
      geom_point() +
      geom_smooth(method="loess",se=TRUE)+
      scale_color_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      scale_fill_manual(name="Treatment",
          breaks=c("TG1","TG2","TG3","TG4","TG5"),
          values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
      coord_cartesian(xlim = c(7, 35)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

alpha <- alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "alpha") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(metric=="richness") %>%
  select(sample,alpha)

mci <- genome_gifts %>%
    to.elements(., GIFT_db) %>% 
    to.functions(., GIFT_db) %>% 
    to.community(., genome_counts_filt %>%
                    mutate_at(vars(-genome),~./sum(.)) %>%
                   column_to_rownames(var="genome"), GIFT_db) %>% 
    as.data.frame() %>% 
    rownames_to_column(var="sample") %>% 
    pivot_longer(-sample, names_to = "trait", values_to = "mci") %>%
    left_join(sample_metadata, by = join_by(sample == sample)) %>% 
    group_by(sample,treatment,day) %>% 
    summarise(mci=mean(mci)) %>% 
    filter(treatment != "TG0") 

inner_join(alpha,mci,by="sample") %>%
      ggplot(aes(y = alpha, x = mci)) +
      geom_smooth(method = "glm", 
                  formula = y~x, 
                  method.args = list(family = gaussian(link = 'log')),
                  color="#999999") +
      geom_point(aes(color=treatment)) +
      scale_color_manual(name="Treatment",
      breaks=c("TG1","TG2","TG3","TG4","TG5"),
      values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
      theme_classic() + 
        theme(legend.position = "none")