Chapter 11 Does antibiotic administration remove the differences found in the two populations?
11.1 Alpha diversity
Shapiro test
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="richness") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.001165318alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="neutral") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.0003462674alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="phylogenetic") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.06610892alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point == "Antibiotics" )alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness", "neutral"))) %>%
filter(metric!="phylogenetic") %>%
filter(time_point == "Antibiotics" ) %>%
ggplot(aes(y = value, x = species, color=species, fill=species)) +
geom_boxplot(width = 0.5,alpha=0.5 ,outlier.shape = NA, show.legend = FALSE) +
geom_jitter(width = 0.2, show.legend = FALSE) +
scale_color_manual(values=c('#008080', "#d57d2c")) +
scale_fill_manual(values=c('#008080', "#d57d2c")) +
stat_compare_means(method = "t.test", show.legend = F, size = 3, label.x = c(1.5))+
facet_grid( ~ metric, scales = "free_y")+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
strip.text = element_text(size = 16, lineheight = 0.6)
) +
ylab("Alpha diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
filter(metric=="phylogenetic") %>%
filter(time_point == "Antibiotics" ) %>%
ggplot(aes(y = value, x = species, color=species, fill=species)) +
geom_boxplot(width = 0.5,alpha=0.5 ,outlier.shape = NA, show.legend = FALSE) +
geom_jitter(width = 0.2, show.legend = FALSE) +
scale_color_manual(values=c('#008080', "#d57d2c")) +
scale_fill_manual(values=c('#008080', "#d57d2c")) +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
facet_grid( ~ metric, scales = "free_y")+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
strip.text = element_text(size = 16, lineheight = 0.6)
) +
ylab("Alpha diversity")
11.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point == "Antibiotics") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point == "Antibiotics")Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$species) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.015377 0.0153774 6.8934 999 0.014 *
Residuals 21 0.046845 0.0022307
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ species,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.361743 | 0.1533845 | 3.80465 | 0.001 |
| Residual | 21 | 7.516224 | 0.8466155 | NA | NA |
| Total | 22 | 8.877966 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$species) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.030649 0.0306488 3.8694 999 0.068 .
Residuals 21 0.166339 0.0079209
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ species,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.787698 | 0.208772 | 5.541022 | 0.001 |
| Residual | 21 | 6.775221 | 0.791228 | NA | NA |
| Total | 22 | 8.562919 | 1.000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$species) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.012165 0.012165 0.9979 999 0.33
Residuals 21 0.256012 0.012191 adonis2(phylo ~ species,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.8970751 | 0.1890846 | 4.896659 | 0.001 |
| Residual | 21 | 3.8472307 | 0.8109154 | NA | NA |
| Total | 22 | 4.7443057 | 1.0000000 | NA | NA |
richness %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta, by = join_by(sample == Tube_code)) %>%
group_by(species, time_point) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species)) +
scale_color_manual(values=c('#008080',"#d57d2c")) +
geom_point(size = 2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme()
neutral %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta, by = join_by(sample == Tube_code))%>%
group_by(species, time_point) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species)) +
scale_color_manual(values=c('#008080',"#d57d2c")) +
geom_point(size = 2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x = "NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme()
phylo %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta, by = join_by(sample == Tube_code))%>%
group_by(species, time_point) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species)) +
scale_color_manual(values = c('#008080', "#d57d2c"))+
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (),x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme() 