Chapter 10 Metabolic capacity index

10.1 Load MG data

load("data/data_mg.Rdata")

# Sequence depth
log_seq_depth <-
  read_counts %>%
  column_to_rownames(var = 'genome') %>%  
  {as.data.frame(log(colSums(.)))} %>%
  rename(seq_depth = 'log(colSums(.))') %>%
  rownames_to_column(var = 'animal_code')

10.2 GIFTs per MAG

10.2.1 Plot GIFTs of all bacteria by element

This plot is not included in Supplementary

gifts_elements %>%
  as.data.frame() %>% 
  rownames_to_column(var = 'genome') %>% 
  pivot_longer(!genome, names_to = 'Code_element', values_to = 'mci') %>%
  inner_join(GIFT_db2, by = 'Code_element') %>%
  ggplot(aes(x = genome,
             y = Code_element,
             fill = mci)) +
  geom_tile(color = '#ffffff') +
  scale_y_discrete(guide = guide_axis(check.overlap = TRUE)) +
  scale_x_discrete(guide = guide_axis(check.overlap = TRUE)) +
  facet_grid(Code_function ~ ., scales = 'free', space = 'free') +
  scale_fill_gradientn(colours = rev(terrain.colors(10))) +
  theme_void(base_size = 18) +
  theme(axis.text.y = element_text(),
        legend.position = 'top')

10.3 GIFTs per community

10.3.1 Calculate community GIFTs

elements_com <- to.community(gifts_elements, total, GIFT_db2) 

funcs_com <- to.community(gifts_functions, total, GIFT_db2)

domains_com <- to.community(gifts_domains, total, GIFT_db2)

save(elements_com,
     funcs_com,
     domains_com,
     file = "data/mci_com.Rdata")

10.3.2 Plot community MCI

The resulting plot corresponds to Figure 2f in the manuscript.

load("data/mci_com.Rdata")

funcs_com %>%
  as.data.frame() %>% 
  rownames_to_column('animal_code') %>%
  rowwise() %>%
  mutate(overall = mean(c_across(B01:D09))) %>%
  select(animal_code, overall) %>%
  left_join(chicken_metadata %>%
            select(animal_code, sampling_time),
            by = 'animal_code') %>%
  ggplot(aes(x = sampling_time,
             y = overall,
             group = sampling_time)) +
  geom_jitter(aes_string(colour = 'sampling_time'), size = 0.3, width = 0.2, alpha = 0.6) +
  geom_boxplot(aes_string(fill = 'sampling_time'), width = 0.4, alpha = 0.5, outlier.color = NA) +
  scale_fill_manual(values = c('#E69F00', '#CC6677', '#56B4E9')) +
  scale_color_manual(values = c('#E69F00', '#CC6677', '#56B4E9')) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.title.x = element_text(size = 8),
        axis.title.y = element_text(size = 8),
        axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 8),
        plot.title =  element_text(size = 10),
        legend.position = 'none')

mean_com_mci <-  
  funcs_com %>%
  as.data.frame() %>% 
  rownames_to_column('animal_code') %>%
  rowwise() %>%
  mutate(mean_mci = mean(c_across(B01:D09))) %>%
  select(animal_code, mean_mci) %>% 
  left_join(chicken_metadata, by = "animal_code") %>% 
  left_join(log_seq_depth, by = 'animal_code')

# Assuming each variable interacts with chicken age
com_mci <-
  lme(mean_mci ~ seq_depth + 
                 sex * age + 
                 breed * age +
                 treatment * age +
                 trial * age,
      random = ~1|pen,
      data = mean_com_mci)

# Best fitted model considered for the manuscript
com_mci <-
  lme(mean_mci ~ seq_depth + sex + breed + treatment + trial + age + trial * age,
      random = ~1|pen,
      data = mean_com_mci)

anova(com_mci)

            numDF denDF   F-value p-value
(Intercept)     1   337 165848.37  <.0001
seq_depth       1   337      2.01  0.1576
sex             1    42      2.97  0.0922
breed           1    42      6.54  0.0142
treatment       2    42      0.13  0.8788
trial           1    42      2.48  0.1226
age             1   337    282.35  <.0001
trial:age       1   337      2.65  0.1045

rm(list = ls())