9 Case study 2

community_kegg_list <- list()
for (strategy in strategies){
  cluster_kegg_strategy <- cluster_kegg %>%
    filter(overall_strategy == strategy) %>%
    select(-overall_strategy) %>%
    column_to_rownames(var="secondary_cluster")
  
  bin_counts_norm_filt_strategy <- bin_counts_norm_filt %>%
    mutate(timepoint=substr(sample,6,7)) %>%
    filter(timepoint %in% c("HT","CT")) %>%
    filter(secondary_cluster %in% rownames(cluster_kegg_strategy)) %>%
    select(secondary_cluster,sample,read_count) %>%
    pivot_wider(names_from = "sample", values_from = "read_count", values_fn = sum) %>%
    mutate(across(where(is.double), ~ ./sum(.)))
  
   community_kegg_table <- c()
   samples <- colnames(bin_counts_norm_filt_strategy)[-1]
   
  for(sample in samples){
    community_kegg <- colSums(sweep(cluster_kegg_strategy[bin_counts_norm_filt_strategy$secondary_cluster,],1,bin_counts_norm_filt_strategy %>% pull({{sample}}), FUN="*"))
    community_kegg_table <- rbind(community_kegg_table,community_kegg)
  }
   rownames(community_kegg_table) <- samples
   
   community_kegg_list[[strategy]] <- community_kegg_table
}

community_kegg_comparison_list <- list()
for (strategy in strategies){
  community_kegg_comparison<- community_kegg_list[[strategy]] %>%
    as.data.frame() %>%
    rownames_to_column(var="sample") %>%
    pivot_longer(!sample,names_to="trait",values_to="value") %>%
    mutate(animal=substr(sample,1,5)) %>%
    mutate(treatment=substr(sample,6,7)) %>%
    group_by(trait) %>%
    mutate(meanvalue=mean(value)) %>%
    filter(meanvalue>0.3) %>%
    mutate(model_result = list(lmerTest::lmer(value ~ treatment + (1 | animal)))) %>%
    ungroup() %>%
    select(trait,model_result) %>%
    unique() %>%
    mutate(estimate = map_dbl(model_result, ~broom.mixed::tidy(.) %>% filter(term == "treatmentHT") %>% pull(estimate))) %>%
    mutate(p_value = map_dbl(model_result, ~broom.mixed::tidy(.) %>% filter(term == "treatmentHT") %>% pull(p.value))) %>%
    mutate(p_value_adj = p.adjust(p_value, method = "bonferroni")) %>%
    select(trait, estimate, p_value_adj)
  
  community_kegg_comparison_list[[strategy]] <- community_kegg_comparison
}

all_p_values <- imap_dfr(community_kegg_comparison_list, ~ .x %>%
  mutate(strategy = .y))

all_p_values %>%
  mutate(trend=ifelse(p_value_adj<0.05 & estimate>0,"hot","neutral")) %>%
  mutate(trend=ifelse(p_value_adj<0.05 & estimate<0,"cold",trend)) %>%
  group_by(trait, trend) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  pivot_wider(values_from="count",names_from="trend") %>%
  tt()

trait	neutral	hot	cold
M00001	10	NA	NA
M00002	10	NA	NA
M00003	6	4	NA
M00004	NA	NA	10
M00005	6	4	NA
M00007	NA	NA	10
M00009	NA	10	NA
M00011	NA	10	NA
M00015	10	NA	NA
M00016	3	NA	7
M00017	8	NA	2
M00018	9	NA	1
M00019	10	NA	NA
M00020	10	NA	NA
M00021	3	NA	7
M00022	NA	NA	10
M00023	10	NA	NA
M00024	10	NA	NA
M00025	10	NA	NA
M00026	10	NA	NA
M00028	NA	NA	10
M00029	NA	NA	10
M00035	10	NA	NA
M00045	10	NA	NA
M00048	10	NA	NA
M00049	10	NA	NA
M00050	3	7	NA
M00051	10	NA	NA
M00052	NA	10	NA
M00053	NA	10	NA
M00060	NA	10	NA
M00061	8	NA	2
M00063	NA	10	NA
M00082	4	NA	6
M00083	10	NA	NA
M00086	2	8	NA
M00089	NA	NA	10
M00093	9	1	NA
M00096	1	NA	9
M00115	NA	10	NA
M00116	NA	10	NA
M00119	NA	10	NA
M00120	10	NA	NA
M00121	10	NA	NA
M00122	10	NA	NA
M00123	NA	10	NA
M00124	NA	10	NA
M00125	NA	10	NA
M00126	NA	10	NA
M00127	NA	10	NA
M00140	10	NA	NA
M00141	NA	10	NA
M00165	10	NA	NA
M00166	8	2	NA
M00167	9	NA	1
M00169	1	9	NA
M00172	1	9	NA
M00173	NA	10	NA
M00178	10	NA	NA
M00179	10	NA	NA
M00183	10	NA	NA
M00260	10	NA	NA
M00307	10	NA	NA
M00308	4	NA	6
M00345	10	NA	NA
M00346	6	4	NA
M00368	7	3	NA
M00373	NA	10	NA
M00377	NA	NA	10
M00432	10	NA	NA
M00525	5	NA	5
M00526	10	NA	NA
M00527	10	NA	NA
M00532	1	9	NA
M00535	10	NA	NA
M00549	10	NA	NA
M00552	10	NA	NA
M00565	NA	NA	10
M00570	10	NA	NA
M00572	10	NA	NA
M00573	NA	10	NA
M00577	NA	10	NA
M00579	1	NA	9
M00609	7	NA	3
M00627	10	NA	NA
M00631	10	NA	NA
M00632	8	2	NA
M00651	10	NA	NA
M00741	NA	10	NA
M00793	10	NA	NA
M00841	NA	10	NA
M00842	NA	10	NA
M00843	NA	10	NA
M00844	NA	NA	10
M00845	NA	NA	10
M00846	10	NA	NA
M00854	NA	NA	10
M00855	5	NA	5
M00866	NA	10	NA
M00874	10	NA	NA

#Maximum number of significances
all_p_values %>%
  mutate(trend=ifelse(p_value_adj<0.05 & estimate>0,"hot","neutral")) %>%
  mutate(trend=ifelse(p_value_adj<0.05 & estimate<0,"cold",trend)) %>%
  group_by(strategy, trend) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  pivot_wider(values_from="count",names_from="trend") %>%
  tt()

strategy	cold	hot	neutral
coassembly_animal	17	31	52
coassembly_timepoint_all	18	27	55
coassembly_timepoint_cage	15	33	52
multicoverage_animal	16	32	52
multicoverage_timepoint_all	17	27	56
multicoverage_timepoint_cage	19	27	54
multisplit_animal	21	32	47
multisplit_timepoint_all	17	33	50
multisplit_timepoint_cage	16	29	55
single_coverage	17	31	52

functional_trends <- all_p_values %>%
  mutate(trend=ifelse(p_value_adj<0.05 & estimate>0,"hot","neutral")) %>%
  mutate(trend=ifelse(p_value_adj<0.05 & estimate<0,"cold",trend)) %>%
  select(trait,strategy,trend) %>%
  pivot_wider(names_from="trait",values_from = "trend") %>%
  column_to_rownames(var="strategy")

functional_tree <- all_p_values %>%
  select(-p_value_adj) %>%
  pivot_wider(names_from="trait",values_from="estimate") %>%
  column_to_rownames(var="strategy") %>%
  dist() %>%
  hclust() %>%
  ggtree()

gheatmap(functional_tree, functional_trends, offset=0.15, colnames=FALSE, width=8) +
  geom_tiplab(size=2, align=TRUE) +
  scale_fill_manual(values=c("#78bacf","#e0766e","#f4f4f4"))