Chapter 7 Bacterial temporal trends

7.1 Load MG data

load("data/data_mg.Rdata")

7.2 Functional ordination - Temporal PCoA with PCA loadings

The resulting plot corresponds to Figure 2c.

gift_pcoa <- 
  gifts_elements %>% 
  as.data.frame() %>%
  vegdist(method="euclidean") %>%
  pcoa()

gift_pcoa_rel_eigen <- gift_pcoa$values$Relative_eig[1:10]

# Get genome positions
gift_pcoa_vectors <- 
  gift_pcoa$vectors %>% #extract vectors
  as.data.frame() %>% 
  select(Axis.1,Axis.2) # keep the first 2 axes


total %>%
  rownames_to_column("genome") %>% 
  pivot_longer(-genome, names_to = "animal_code", values_to = "count") %>%
  left_join(chicken_metadata, by = join_by(animal_code == animal_code)) %>% 
  filter(!is.na(sampling_time)) %>% 
  group_by(genome, trial, sampling_time) %>% 
  summarise(count = mean(count)) %>% 
  left_join(gift_pcoa_vectors %>% rownames_to_column(var = "genome"), by = "genome") %>%
  left_join(genome_taxonomy %>% select(genome, phylum, order), by = "genome") %>%
  ggplot() +
      theme_minimal() +
      #genome positions
      geom_point(aes(x = Axis.1, y = Axis.2, color = order, size = count), 
                 alpha = 0.8, shape = 16) +
      # colors
      scale_color_manual(values = order_colors) +    
      # scale_color_manual(values=phylum_colors) +
      scale_size_continuous(range = c(1,10)) +
      # scale adjustments
      scale_x_continuous(name = paste0("PCoA1 (",round(gift_pcoa_rel_eigen[1]*100, digits = 2), " %)")) +
      scale_y_continuous(name = paste0("PCoA2 (",round(gift_pcoa_rel_eigen[2]*100, digits = 2), " %)")) +
      facet_grid(trial ~ sampling_time) +
      theme(legend.position = "none",
            axis.line = element_line(color = "grey30"),
            panel.border = element_rect(color = "grey80", fill = NA),
            axis.ticks = element_line(color = "grey30"),
            panel.grid = element_line(color = "grey80"),
            strip.background = element_rect(fill = "grey80", color = NA))

7.3 Compositional barplots

7.3.1 Compositional barplot at order level

physeq_comp_rare_order <-
  aggregate_rare(
    phylo_seq,
    level = "order",
    detection = 0.01 / 100,
    prevalence = 50 / 100,
    include.lowest = TRUE)
  
phy_order_df <- 
  psmelt(physeq_comp_rare_order) %>%
  select(OTU, Sample, trial_age, Abundance) %>%
  rename(order = 'OTU') %>%
  separate_wider_delim(trial_age, "_", names = c("trial", "sampling_time"))

# separating by sampling time and trial
phy_order_df %>%
  mutate(order = factor(order, 
                        levels = c('Veillonellales', 'Acidaminococcales',
                                   'Selenomonadales', 'UBA4068', 'UBA1212',
                                   'Monoglobales', 'Peptostreptococcales',
                                   'Monoglobales_A', 'Clostridiales',
                                   'Lachnospirales', 'Christensenellales',
                                   'TANB77', 'Oscillospirales', 'RFN20', 
                                   'Bacillales','Acholeplasmatales', 
                                   'Erysipelotrichales', 'ML615J-28', 
                                   'Lactobacillales','RF39','Deferribacterales',
                                   'Actinomycetales', 'Coriobacteriales',
                                   'Gastranaerophilales','Flavobacteriales',
                                   'Bacteroidales', 'Rs-D84', 'Burkholderiales',
                                   'RF32','Enterobacterales', 
                                   'Desulfovibrionales', 'Verrucomicrobiales',
                                   'Opitutales','Victivallales', 
                                   'Methanomassiliicoccales', 'Synergistales',
                                   'Campylobacterales','Saccharimonadales',
                                   'Methanobacteriales', 'Methanomicrobiales'))) %>%
  mutate(sampling_time = str_replace(sampling_time,
                                     pattern = "7",
                                     replacement = "Day 7"),
         sampling_time = str_replace(sampling_time,
                                     pattern = "21",
                                     replacement = "Day 21"),
         sampling_time = str_replace(sampling_time,
                                     pattern = "35",
                                     replacement = "Day 35")) %>%
  mutate(sampling_time = factor(sampling_time,
                                levels = c("Day 7", "Day 21", "Day 35"))) %>%
  ggplot(aes(x = Abundance, y = Sample,
             fill = order, group = order)) +
  geom_bar(stat = "identity", position = "stack", width = 1,
         color ="darkgrey", size = 0.005) +
  facet_nested(sampling_time * trial ~., scales = "free_y") +
  scale_color_manual('Order', values = order_colors) +
  scale_fill_manual('Order', values = order_colors) +
  ylab('Relative abundance') +
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        panel.spacing.x = unit(0.1, "cm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 10),
        strip.text = element_text(size = 8),
        legend.position = "right")  +
  guides(fill = guide_legend(ncol = 1))

rm(list = ls())