Chapter 6 Compositional analysis

load("data/data_filtered.Rdata")
sample_metadata %>% 
  group_by(region, Bd_status, sex) %>% 
  summarise(count=n())
`summarise()` has grouped output by 'region', 'Bd_status'. You can override using the `.groups` argument.
# A tibble: 7 × 4
# Groups:   region, Bd_status [4]
  region          Bd_status   sex    count
  <chr>           <chr>       <chr>  <int>
1 Matra mountains Bd_negative Female    10
2 Matra mountains Bd_negative Male      11
3 Matra mountains Bd_positive Female     2
4 Matra mountains Bd_positive Male       8
5 None            Bd_negative Female     1
6 None            Bd_negative Male       5
7 None            Bd_positive Male       2
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(!is.na(count)) %>%
  ggplot(aes(y=count,x=sample, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + 
    scale_fill_manual(values=phylum_colors) +
  facet_grid(~ Bd_status, scale="free", space="free")+
    guides(fill = guide_legend(ncol = 1)) +
    theme(
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
    strip.text = element_text(size = 8, lineheight = 0.6),
    strip.placement = "outside",
    axis.text.x = element_blank(), axis.ticks.x = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

6.1 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  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,Bd_status) %>%
  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),
    positive_mean = mean(relabun[Bd_status == "Bd_positive"] * 100, na.rm = T),
    positive_sd = sd(relabun[Bd_status == "Bd_positive"] * 100, na.rm = T),
    negative_mean = mean(relabun[Bd_status == "Bd_negative"] * 100, na.rm = T),
    negative_sd = sd(relabun[Bd_status == "Bd_negative"] * 100, na.rm = T)
  ) %>%
  mutate(
    total = str_c(round(total_mean, 3), "±", round(total_sd, 3)),
    Positive = str_c(round(positive_mean, 3), "±", round(positive_sd, 3)),
    Negative = str_c(round(negative_mean, 3), "±", round(negative_sd, 3))
    ) %>%
  arrange(-total_mean) %>%
  dplyr::select(phylum, total, Positive, Negative) %>%
  tt()
phylum total Positive Negative
p__Pseudomonadota 73.437±23.223 70.928±18.007 74.552±25.433
p__Actinomycetota 6.249±13.878 3.852±7.991 7.314±15.832
p__Bacteroidota 5.885±7.341 9.221±10.595 4.402±4.883
p__Verrucomicrobiota 5.635±16 2.96±5.359 6.824±18.9
p__Patescibacteria 5.435±14.284 7.971±13.383 4.308±14.77
p__Cyanobacteriota 1.422±5.766 2.147±7.203 1.1±5.126
p__Desulfobacterota_F 0.416±1.031 0.671±1.538 0.303±0.713
p__Campylobacterota 0.336±1.012 0.714±1.617 0.168±0.545
p__Bacillota 0.311±1.118 0.573±1.879 0.195±0.536
p__Chloroflexota 0.247±0.875 0.196±0.68 0.269±0.96
p__Gemmatimonadota 0.216±0.603 0.055±0.19 0.288±0.706
p__Thermoproteota 0.208±0.908 0.328±1.138 0.155±0.805
p__Nitrospirota 0.071±0.443 0.23±0.798 0±0
p__Halobacteriota 0.046±0.232 0±0 0.067±0.277
p__Bdellovibrionota 0.037±0.234 0±0 0.054±0.281
p__Acidobacteriota 0.025±0.157 0.082±0.284 0±0
p__Bacillota_A 0.022±0.136 0.071±0.245 0±0