Chapter 5 Alpha diversity

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

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  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(treatment != "T0") %>% 
  filter(metric=="richness") %>%
      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("T1","T3"),
          values=c("#6A9AC3","#AFD699")) +
      scale_fill_manual(name="Treatment",
          breaks=c("T1","T3"),
          values=c("#6A9AC350","#AFD69950")) +
      facet_wrap(. ~ trial, 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)
      )

#Neutral
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(treatment != "T0") %>% 
  filter(metric=="neutral") %>%
      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("T1","T3"),
          values=c("#6A9AC3","#AFD699")) +
      scale_fill_manual(name="Treatment",
          breaks=c("T1","T3"),
          values=c("#6A9AC350","#AFD69950")) +
      facet_wrap(. ~ trial, 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)
      )

#Phylogenetic
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(treatment != "T0") %>% 
  filter(metric=="phylogenetic") %>%
      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("T1","T3"),
          values=c("#6A9AC3","#AFD699")) +
      scale_fill_manual(name="Treatment",
          breaks=c("T1","T3"),
          values=c("#6A9AC350","#AFD69950")) +
      facet_wrap(. ~ trial, 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)
      )

#Functional
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  filter(treatment != "T0") %>% 
  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("T1","T3"),
          values=c("#6A9AC3","#AFD699")) +
      scale_fill_manual(name="Treatment",
          breaks=c("T1","T3"),
          values=c("#6A9AC350","#AFD69950")) +
      facet_wrap(. ~ trial, 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)
      )