Chapter 4 Community composition

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(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  ggplot(., 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_grid(. ~ location,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x =element_blank(),
          axis.title.x = element_blank(),
          axis.text.y = element_text(size=8),
          axis.title.y = element_text(size=12),
          panel.background = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
          panel.border = element_rect(colour = "black", fill = NA),
          strip.background = element_rect(fill = "white", color = "black"),
          strip.text = element_text(size = 12, lineheight = 0.6)) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

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(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  ggplot(., 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(. ~ location + origin,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x =element_blank(),
          axis.title.x = element_blank(),
          axis.text.y = element_text(size=8),
          axis.title.y = element_text(size=12),
          panel.background = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
          panel.border = element_rect(colour = "black", fill = NA),
          strip.background = element_rect(fill = "white", color = "black"),
          strip.text = element_text(size = 12, lineheight = 0.6)) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

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(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  ggplot(., 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(. ~ origin + location,  scales="free_x") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x =element_blank(),
          axis.title.x = element_blank(),
          axis.text.y = element_text(size=8),
          axis.title.y = element_text(size=12),
          panel.background = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
          panel.border = element_rect(colour = "black", fill = NA),
          strip.background = element_rect(fill = "white", color = "black"),
          strip.text = element_text(size = 10, lineheight = 0.6)) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

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

phylum_summary %>%
    group_by(phylum) %>%
    summarise(total_mean=mean(relabun*100, na.rm=T),
              total_sd=sd(relabun*100, na.rm=T),
              domestic_mean=mean(relabun[origin=="Domestic"]*100, na.rm=T),
              domestic_sd=sd(relabun[origin=="Domestic"]*100, na.rm=T),
              feral_mean=mean(relabun[origin=="Feral"]*100, na.rm=T),
              feral_sd=sd(relabun[origin=="Feral"]*100, na.rm=T)) %>%
    mutate(total=str_c(round(total_mean,2),"±",round(total_sd,2)),
           domestic=str_c(round(domestic_mean,2),"±",round(domestic_sd,2)),
           feral=str_c(round(feral_mean,2),"±",round(feral_sd,2))) %>% 
    arrange(-total_mean) %>% 
    select(phylum,total,domestic,feral) %>% 
    tt()

phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    dplyr::select(phylum) %>%
    pull()

phylum_summary %>%
  filter(phylum %in% phylum_arrange) %>%
  mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
    filter(relabun > 0) %>%
  ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum, fill=phylum)) +
  scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        scale_fill_manual(values=phylum_colors[-8]) +
        geom_boxplot(alpha=0.2)+
        geom_jitter(alpha=0.5) + 
        facet_nested(. ~ location)+ 
        theme_minimal() + 
        theme(legend.position="none",
              strip.text.x = element_text(size = 14, color="black",face="bold"),
              axis.text.x = element_text(vjust = 0.5, size = 10),
              axis.text.y = element_text(size = 12),
              axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
              axis.title = element_text(size = 14, face = "bold"),
              axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
              axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0))
              ) +
        labs(y="Phylum",x="Relative abundance")

phylum_summary %>%
  filter(phylum %in% phylum_arrange) %>%
  mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
    filter(relabun > 0) %>%
  ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum, fill=phylum)) +
  scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        scale_fill_manual(values=phylum_colors[-8]) +
        geom_boxplot(alpha=0.2)+
        geom_jitter(alpha=0.5) + 
        theme_minimal() + 
        theme(legend.position="none") +
        labs(y="Phylum",x="Relative abundance")

family_summary %>%
  left_join(genome_metadata %>% dplyr::select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
  left_join(sample_metadata,by=join_by(sample==sample)) %>%
  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, fill=phylum)) +
      scale_color_manual(values=phylum_colors[-8]) +
      scale_fill_manual(values=phylum_colors[-8]) +
      geom_boxplot(alpha=0.2)+
      geom_jitter(alpha=0.5) + 
      facet_grid(.~location)+
      theme_minimal() +
      theme(legend.position="none",
                  strip.text.x = element_text(size = 14, color="black",face="bold"),
                  axis.text.x = element_text(vjust = 0.5, size = 6),
                  axis.text.y = element_text(size = 12),
                  axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
                  axis.title = element_text(size = 14, face = "bold"),
                  axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
                  axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)))+
      labs(y="Family", x="Relative abundance", color="Phylum")

family_summary %>%
    left_join(genome_metadata %>% dplyr::select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==sample)) %>%
    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, fill=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        scale_fill_manual(values=phylum_colors[-8]) +
        geom_boxplot(alpha=0.2)+
        geom_jitter(alpha=0.5) + 
        facet_grid(.~location)+
        theme_minimal() +
        theme(legend.position="none",
                    strip.text.x = element_text(size = 14, color="black",face="bold"),
                    axis.text.x = element_text(vjust = 0.5, size = 6),
                    axis.text.y = element_text(size = 12),
                    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
                    axis.title = element_text(size = 14, face = "bold"),
                    axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
                    axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)))+
        labs(y="Family", x="Relative abundance", color="Phylum")

genus_summary %>%
  left_join(genome_metadata %>% dplyr::select(genus,phylum) %>% unique(),by="genus") %>%
  left_join(sample_metadata,by=join_by(sample==sample)) %>%
  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, fill=phylum)) +
  scale_color_manual(values=phylum_colors[-8]) +
  scale_fill_manual(values=phylum_colors[-8]) +
  geom_boxplot(alpha=0.2)+
  geom_jitter(alpha=0.5) + 
  facet_grid(.~location)+
  theme_minimal() +
  theme(legend.position="none",
              strip.text.x = element_text(size = 14, color="black",face="bold"),
              axis.text.x = element_text(vjust = 0.5, size = 6),
              axis.text.y = element_text(size = 12),
              axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
              axis.title = element_text(size = 14, face = "bold"),
              axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
              axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)))+
  labs(y="Family", x="Relative abundance", color="Phylum")

genus_summary %>%
  left_join(genome_metadata %>% dplyr::select(genus,phylum) %>% unique(),by="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, fill=phylum)) +
  scale_color_manual(values=phylum_colors[-8]) +
  scale_fill_manual(values=phylum_colors[-8]) +
  geom_boxplot(alpha=0.2)+
  geom_jitter(alpha=0.5) + 
  facet_grid(.~origin)+
  theme_minimal() +
  theme(legend.position="none",
              strip.text.x = element_text(size = 14, color="black",face="bold"),
              axis.text.x = element_text(vjust = 0.5, size = 6),
              axis.text.y = element_text(size = 12),
              axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
              axis.title = element_text(size = 14, face = "bold"),
              axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
              axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)))+
  labs(y="Family", x="Relative abundance", color="Phylum")