Chapter 5 Community composition

# Merge data frames based on sample
# transplants_metadata <- sample_metadata %>%
#   mutate(Tube_code = str_remove_all(Tube_code, "_a"))
# transplants_metadata$newID <- paste(transplants_metadata$Tube_code,
#                                     "_",
#                                     transplants_metadata$individual)

merged_data <- genome_counts_filt %>%
  mutate_at(vars(-genome),  ~ . / sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  filter(individual!="NA")

ggplot(merged_data, aes(
  x = sample,
  y = count,
  fill = phylum,
  group = phylum
)) + #grouping enables keeping the same sorting of taxonomic units
  geom_bar(stat = "identity",
           colour = "white",
           linewidth = 0.1) + #plot stacked bars with white borders
  scale_fill_manual(values = phylum_colors) +
  facet_nested(. ~ factor (time_point, level=c("Wild", "Acclimation", "Antibiotics", "Inoculum1", "Inoculum2", "FMT1", "FMT2")) + type , scales = "free") + #facet per day and treatment
  guides(fill = guide_legend(ncol = 1)) +
  labs(fill = "Phylum", y = "Relative abundance", x = "Sample") +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      size = 0
    ),
    strip.text.x = element_text(
      size = 14,
      colour = "black",
      face = "bold"
    ),
    strip.background = element_rect(fill ="lightgrey"),
    axis.title = element_text(size = 18, face = "bold"),
    panel.background = element_blank(),
    legend.title = element_text(size = 20, face = "bold"),
    legend.text = element_text(size = 16)
  )

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Genus", x="Relative abundance", color="Phylum")