Chapter 5 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = keep.tip(genome_tree,tip=rownames(.)))
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method="gower")
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, dist = dist[rownames(.), rownames(.)])
save(beta_q0n, beta_q1n, beta_q1p, file="data/beta_diversity.Rdata")5.1 Permanova
#Richness
betadisper(beta_q0n$C, sample_metadata$treatment[sample_metadata$treatment != "TG0"]) %>% permutest(., pairwise = TRUE)
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 4 0.11176 0.027939 1.431 999 0.234
Residuals 107 2.08914 0.019525
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
TG1 TG2 TG3 TG4 TG5
TG1 0.489000 0.067000 0.223000 0.627
TG2 0.482966 0.185000 0.454000 0.756
TG3 0.080345 0.175034 0.448000 0.100
TG4 0.209775 0.478109 0.456447 0.316
TG5 0.643690 0.754593 0.110414 0.307687 adonis2(beta_q0n$C ~ treatment * day,
data = sample_metadata %>% filter(treatment != "TG0") %>% arrange(match(sample,labels(beta_q1n$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 9 | 7.484587 | 0.6208242 | 18.55606 | 0.001 |
| Residual | 102 | 4.571301 | 0.3791758 | NA | NA |
| Total | 111 | 12.055888 | 1.0000000 | NA | NA |
#Neutral diversity
betadisper(beta_q1n$C, sample_metadata$treatment[sample_metadata$treatment != "TG0"]) %>% permutest(., pairwise = TRUE)
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 4 0.03083 0.0077065 0.3238 999 0.865
Residuals 107 2.54696 0.0238034
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
TG1 TG2 TG3 TG4 TG5
TG1 0.88500 0.40900 0.50400 0.775
TG2 0.88020 0.42800 0.57000 0.871
TG3 0.36285 0.42428 0.74700 0.508
TG4 0.50197 0.59266 0.75648 0.694
TG5 0.76469 0.87585 0.52248 0.71357 adonis2(beta_q1n$C ~ treatment * day,
data = sample_metadata %>% filter(treatment != "TG0") %>% arrange(match(sample,labels(beta_q1n$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 9 | 8.812652 | 0.6062983 | 17.45327 | 0.001 |
| Residual | 102 | 5.722522 | 0.3937017 | NA | NA |
| Total | 111 | 14.535174 | 1.0000000 | NA | NA |
#Phylogenetic diversity
betadisper(beta_q1p$C, sample_metadata$treatment[sample_metadata$treatment != "TG0"]) %>% permutest(., pairwise = TRUE)
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 4 0.007294 0.0018234 0.6588 999 0.604
Residuals 107 0.296161 0.0027679
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
TG1 TG2 TG3 TG4 TG5
TG1 0.53300 0.14700 0.19000 0.265
TG2 0.52164 0.33700 0.40600 0.498
TG3 0.15794 0.32125 0.89900 0.891
TG4 0.21331 0.40937 0.90278 0.981
TG5 0.27885 0.48101 0.88912 0.97933 adonis2(beta_q1p$C ~ treatment * day,
data = sample_metadata %>% filter(treatment != "TG0") %>% arrange(match(sample,labels(beta_q1p$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 9 | 0.5927034 | 0.5435966 | 13.4985 | 0.001 |
| Residual | 102 | 0.4976333 | 0.4564034 | NA | NA |
| Total | 111 | 1.0903367 | 1.0000000 | NA | NA |
#Functional diversity
betadisper(beta_q1f$C, sample_metadata$treatment[sample_metadata$treatment != "TG0"]) %>% permutest(., pairwise = TRUE)
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 4 0.001863 0.00046574 0.9798 999 0.443
Residuals 107 0.050863 0.00047536
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
TG1 TG2 TG3 TG4 TG5
TG1 0.65600 0.22600 0.47200 0.368
TG2 0.66449 0.11700 0.21900 0.211
TG3 0.20079 0.11202 0.52200 0.694
TG4 0.43259 0.24866 0.51631 0.816
TG5 0.34760 0.20468 0.69546 0.80418 adonis2(beta_q1f$C ~ treatment * day,
data = sample_metadata %>% filter(treatment != "TG0") %>% arrange(match(sample,labels(beta_q1f$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 9 | 0.01178871 | 0.120672 | 1.555297 | 0.277 |
| Residual | 102 | 0.08590348 | 0.879328 | NA | NA |
| Total | 111 | 0.09769219 | 1.000000 | NA | NA |
5.2 Richness
nmds_richness <- beta_q0n$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment,day) %>%
mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
ungroup() 5.2.1 All time points
nmds_richness %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment, fill = treatment, shape = as.factor(day))) +
scale_color_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
scale_fill_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
scale_shape_manual(name="Day",
breaks=c(7,14,21,28,35),
values=c(21,22,23,24,25)) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(shape="Individual")
nmds_richness %>%
select(x_cen,y_cen,treatment,day) %>%
unique() %>%
arrange(treatment,day) %>%
ggplot(aes(x = x_cen, y = y_cen, color = treatment, fill = treatment)) +
scale_color_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
scale_fill_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
geom_point(size = 4) +
geom_path(aes(xend=c(tail(x_cen, n=-1), NA),
yend=c(tail(y_cen, n=-1), NA),
group=treatment),
size = 2,
alpha=0.5,
arrow=arrow(type = "closed", length=unit(0.5,"cm"))
) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
#geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(shape="Individual")
nmds_richness %>%
filter(day == 7) %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment, fill = treatment, shape = as.factor(day))) +
scale_color_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
scale_fill_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
scale_shape_manual(name="Day",
breaks=c(7,14,21,28,35),
values=c(21,22,23,24,25)) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(shape="Individual")
## Neutral diversity
nmds_neutral <- beta_q1n$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment,day) %>%
mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
ungroup() 5.2.2 All time points
nmds_neutral %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment, fill = treatment, shape = as.factor(day))) +
scale_color_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
scale_fill_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
scale_shape_manual(name="Day",
breaks=c(7,14,21,28,35),
values=c(21,22,23,24,25)) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(shape="Individual")
5.2.3 Trajectories
nmds_neutral %>%
select(x_cen,y_cen,treatment,day) %>%
unique() %>%
arrange(treatment,day) %>%
ggplot(aes(x = x_cen, y = y_cen, color = treatment, fill = treatment)) +
scale_color_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
scale_fill_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
geom_point(size = 4) +
geom_path(aes(xend=c(tail(x_cen, n=-1), NA),
yend=c(tail(y_cen, n=-1), NA),
group=treatment),
size = 2,
alpha=0.5,
arrow=arrow(type = "closed", length=unit(0.5,"cm"))
) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
#geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(shape="Individual")
5.3 Day 7
nmds_neutral %>%
filter(day == 7) %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment, fill = treatment, shape = as.factor(day))) +
scale_color_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE", "#6A9AC3","#97D8C4","#AFD699","#F3B942")) +
scale_fill_manual(name="Treatment",
breaks=c("TG1","TG2","TG3","TG4","TG5"),
values=c("#4059AE50", "#6A9AC350","#97D8C450","#AFD69950","#F3B94250")) +
scale_shape_manual(name="Day",
breaks=c(7,14,21,28,35),
values=c(21,22,23,24,25)) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(shape="Individual")