Alpha diversity
alpha_div <- tibble()
for (strategy in strategies){
bin_counts_norm_filt_strategy <- bin_counts_norm_filt %>%
filter(overall_strategy == {{strategy}}) %>%
select(secondary_cluster,sample,read_count) %>%
pivot_wider(names_from = "sample", values_from = "read_count", values_fn = sum) %>%
column_to_rownames(var="secondary_cluster")
cluster_tree_strategy <- keep.tip(cluster_tree, tip=rownames(bin_counts_norm_filt_strategy))
q0n <- bin_counts_norm_filt_strategy %>%
hilldiv(.,q=0) %>% as.numeric()
q1n <- bin_counts_norm_filt_strategy %>%
hilldiv(.,q=1) %>% as.numeric()
q1p <- bin_counts_norm_filt_strategy %>%
hilldiv(.,q=1,tree=cluster_tree_strategy) %>% as.numeric()
alpha_div_strategy <- tibble(strategy=strategy,sample=colnames(bin_counts_norm_filt_strategy),richness=q0n,neutral=round(q1n,3),phylo=round(q1p,3))
alpha_div<- bind_rows(alpha_div,alpha_div_strategy)
}
save(alpha_div,file="results/alpha_diversity.Rdata")
load("results/alpha_diversity.Rdata")
alpha_div %>%
pivot_longer(!c(strategy,sample),names_to = "metric",values_to = "value") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylo"))) %>%
mutate(strategy=factor(strategy,levels=rev(strategies))) %>%
ggplot(aes(y=strategy, x=value, color=strategy))+
geom_boxplot(outlier.shape = NA)+
geom_jitter(alpha=0.1)+
scale_color_manual(values=strategy_colors)+
facet_nested(. ~ metric, scales = "free") +
theme_minimal() +
labs(x="Diversity",y="Strategy") +
theme(legend.position = "none")
load("results/alpha_diversity.Rdata")
alpha_test <- alpha_div %>%
pivot_longer(!c(strategy,sample),names_to = "metric",values_to = "value") %>%
mutate(strategy=factor(strategy,levels=strategies)) %>%
filter(metric=="richness") %>%
lmerTest::lmer(value ~ strategy + (1 | sample), data = .)
alpha_test %>%
broom.mixed::tidy()
# A tibble: 12 × 8
effect group term estimate std.error statistic df p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 65.8 0.872 75.5 325. 3.42e-208
2 fixed <NA> strategymulticoverage_animal 0.867 1.23 0.703 325. 4.82e- 1
3 fixed <NA> strategymulticoverage_timepoint_all 0.787 1.23 0.638 325. 5.24e- 1
4 fixed <NA> strategymulticoverage_timepoint_cage -7.82 1.23 -6.34 325. 7.48e- 10
5 fixed <NA> strategycoassembly_animal 12.5 0.290 43.2 1192. 3.27e-246
6 fixed <NA> strategycoassembly_timepoint_all 13.3 0.290 45.8 1192. 4.45e-265
7 fixed <NA> strategycoassembly_timepoint_cage 11.9 0.290 40.9 1192. 1.78e-229
8 fixed <NA> strategymultisplit_animal -8.67 1.23 -7.04 325. 1.17e- 11
9 fixed <NA> strategymultisplit_timepoint_all -1.93 1.23 -1.56 325. 1.19e- 1
10 fixed <NA> strategymultisplit_timepoint_cage -3.08 1.23 -2.50 325. 1.30e- 2
11 ran_pars sample sd__(Intercept) 10.4 NA NA NA NA
12 ran_pars Residual sd__Observation 2.51 NA NA NA NA
alpha_test %>%
anova() %>%
broom::tidy()
# A tibble: 1 × 7
term sumsq meansq NumDF DenDF statistic p.value
<chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 strategy 31409. 3490. 9 893. 553. 0
load("results/alpha_diversity.Rdata")
alpha_test <- alpha_div %>%
pivot_longer(!c(strategy,sample),names_to = "metric",values_to = "value") %>%
mutate(strategy=factor(strategy,levels=strategies)) %>%
filter(metric=="neutral") %>%
lmerTest::lmer(value ~ strategy + (1 | sample), data = .)
alpha_test %>%
broom.mixed::tidy()
# A tibble: 12 × 8
effect group term estimate std.error statistic df p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 29.6 0.630 47.0 303. 3.51e-141
2 fixed <NA> strategymulticoverage_animal -0.610 0.892 -0.685 303. 4.94e- 1
3 fixed <NA> strategymulticoverage_timepoint_all -0.519 0.892 -0.582 303. 5.61e- 1
4 fixed <NA> strategymulticoverage_timepoint_cage -3.18 0.892 -3.56 303. 4.23e- 4
5 fixed <NA> strategycoassembly_animal 1.61 0.0921 17.5 1192. 3.02e- 61
6 fixed <NA> strategycoassembly_timepoint_all 1.59 0.0921 17.2 1192. 1.67e- 59
7 fixed <NA> strategycoassembly_timepoint_cage 1.61 0.0921 17.4 1192. 7.53e- 61
8 fixed <NA> strategymultisplit_animal -3.95 0.892 -4.43 303. 1.32e- 5
9 fixed <NA> strategymultisplit_timepoint_all -0.768 0.892 -0.861 303. 3.90e- 1
10 fixed <NA> strategymultisplit_timepoint_cage -1.56 0.892 -1.74 303. 8.21e- 2
11 ran_pars sample sd__(Intercept) 7.68 NA NA NA NA
12 ran_pars Residual sd__Observation 0.798 NA NA NA NA
alpha_test %>%
anova() %>%
broom::tidy()
# A tibble: 1 × 7
term sumsq meansq NumDF DenDF statistic p.value
<chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 strategy 1900. 211. 9 893. 332. 1.17e-277
load("results/alpha_diversity.Rdata")
alpha_test <- alpha_div %>%
pivot_longer(!c(strategy,sample),names_to = "metric",values_to = "value") %>%
mutate(strategy=factor(strategy,levels=strategies)) %>%
filter(metric=="phylo") %>%
lmerTest::lmer(value ~ strategy + (1 | sample), data = .)
alpha_test %>%
broom.mixed::tidy()
# A tibble: 12 × 8
effect group term estimate std.error statistic df p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 4.74 0.0572 82.8 313. 6.02e-215
2 fixed <NA> strategymulticoverage_animal -0.0705 0.0809 -0.871 313. 3.84e- 1
3 fixed <NA> strategymulticoverage_timepoint_all -0.0425 0.0809 -0.526 313. 6.00e- 1
4 fixed <NA> strategymulticoverage_timepoint_cage -0.226 0.0809 -2.80 313. 5.46e- 3
5 fixed <NA> strategycoassembly_animal 0.118 0.0142 8.28 1192. 3.38e- 16
6 fixed <NA> strategycoassembly_timepoint_all 0.189 0.0142 13.3 1192. 1.15e- 37
7 fixed <NA> strategycoassembly_timepoint_cage 0.0591 0.0142 4.15 1192. 3.53e- 5
8 fixed <NA> strategymultisplit_animal -0.256 0.0809 -3.16 313. 1.74e- 3
9 fixed <NA> strategymultisplit_timepoint_all -0.0337 0.0809 -0.417 313. 6.77e- 1
10 fixed <NA> strategymultisplit_timepoint_cage -0.0757 0.0809 -0.935 313. 3.50e- 1
11 ran_pars sample sd__(Intercept) 0.690 NA NA NA NA
12 ran_pars Residual sd__Observation 0.123 NA NA NA NA
alpha_test %>%
anova() %>%
broom::tidy()
# A tibble: 1 × 7
term sumsq meansq NumDF DenDF statistic p.value
<chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 strategy 10.2 1.13 9 893. 74.5 2.05e-102