Chapter 8 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 = genome_tree)
beta_q1f <- genome_counts_filt %>%
filter(genome %in% labels(dist)[[1]]) %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, dist = dist)8.1 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)) %>%
mutate(Bd_status=factor(Bd_status,levels=c("Bd_positive", "Bd_negative")))%>%
group_by(Bd_status) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = Bd_status)) +
geom_point(size = 3, alpha = 0.9) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.5, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(face = "bold", size = 12),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Bd status')
beta_q0n_mat <- as.matrix(beta_q0n$S)
beta_q0n_filtered <- beta_q0n_mat[rownames(beta_q0n_mat) %in% sample_metadata$sample,
colnames(beta_q0n_mat) %in% sample_metadata$sample]
beta_q0n_dis <- as.dist(beta_q0n_filtered)
sample_metadata_ordered <- sample_metadata %>%
arrange(match(sample,labels(beta_q0n_dis))) %>%
column_to_rownames(., "sample")
all(sample_metadata_ordered$sample == labels(beta_q0n_dis))[1] 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 1 0.000347 0.00034692 0.2568 999 0.61
Residuals 37 0.049992 0.00135112
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Bd_negative Bd_positive
Bd_negative 0.608
Bd_positive 0.61536 adonis2(beta_q0n_dis ~ Bd_status+region+sex,
data = sample_metadata_ordered,
permutations = 999, by="margin") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Bd_status | 1 | 0.3829697 | 0.02583199 | 1.042168 | 0.362 |
| region | 1 | 1.2485628 | 0.08421778 | 3.397691 | 0.001 |
| sex | 1 | 0.3330269 | 0.02246326 | 0.906260 | 0.562 |
| Residual | 35 | 12.8615873 | 0.86753690 | NA | NA |
| Total | 38 | 14.8254065 | 1.00000000 | NA | NA |
8.2 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)) %>%
mutate(Bd_status=factor(Bd_status,levels=c("Bd_positive", "Bd_negative")))%>%
group_by(Bd_status) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = Bd_status)) +
geom_point(size = 3, alpha = 0.9) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.5, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(face = "bold", size = 12),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Bd status')
beta_q1n_mat <- as.matrix(beta_q1n$S)
beta_q1n_filtered <- beta_q1n_mat[rownames(beta_q1n_mat) %in% sample_metadata$sample,
colnames(beta_q1n_mat) %in% sample_metadata$sample]
beta_q1n_dis <- as.dist(beta_q1n_filtered)
sample_metadata_ordered <- sample_metadata %>%
arrange(match(sample,labels(beta_q1n_dis))) %>%
column_to_rownames(., "sample")
all(sample_metadata_ordered$sample == labels(beta_q1n_dis))[1] 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 1 0.01088 0.010878 0.1925 999 0.662
Residuals 37 2.09067 0.056505
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Bd_negative Bd_positive
Bd_negative 0.66
Bd_positive 0.66339 adonis2(beta_q1n_dis ~ Bd_status+region+sex,
data = sample_metadata_ordered,
permutations = 999, by="margin") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Bd_status | 1 | 0.2720914 | 0.02279371 | 0.9408119 | 0.459 |
| region | 1 | 1.2731225 | 0.10665235 | 4.4020830 | 0.001 |
| sex | 1 | 0.1900787 | 0.01592332 | 0.6572362 | 0.788 |
| Residual | 35 | 10.1223185 | 0.84796953 | NA | NA |
| Total | 38 | 11.9371252 | 1.00000000 | NA | NA |
8.3 Phylogenetic
beta_q1p$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)) %>%
mutate(Bd_status=factor(Bd_status,levels=c("Bd_positive", "Bd_negative")))%>%
group_by(Bd_status) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = Bd_status)) +
geom_point(size = 3, alpha = 0.9) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.5, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(face = "bold", size = 12),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Bd status')
beta_q1p_mat <- as.matrix(beta_q1p$S)
beta_q1p_filtered <- beta_q1p_mat[rownames(beta_q1p_mat) %in% sample_metadata$sample,
colnames(beta_q1p_mat) %in% sample_metadata$sample]
beta_q1p_dis <- as.dist(beta_q1p_filtered)
sample_metadata_ordered <- sample_metadata %>%
arrange(match(sample,labels(beta_q1p_dis))) %>%
column_to_rownames(., "sample")
all(sample_metadata_ordered$sample == labels(beta_q1p_dis))[1] 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 1 0.01925 0.019254 0.5297 999 0.487
Residuals 37 1.34478 0.036345
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Bd_negative Bd_positive
Bd_negative 0.502
Bd_positive 0.47129 adonis2(beta_q1p_dis ~ Bd_status+region+sex,
data = sample_metadata_ordered,
permutations = 999, by="margin") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Bd_status | 1 | 0.17324888 | 0.038445404 | 1.5323163 | 0.149 |
| region | 1 | 0.29210497 | 0.064820583 | 2.5835503 | 0.026 |
| sex | 1 | 0.04159565 | 0.009230429 | 0.3678967 | 0.900 |
| Residual | 35 | 3.95721881 | 0.878140603 | NA | NA |
| Total | 38 | 4.50636128 | 1.000000000 | NA | NA |
8.4 Functional
beta_q1f$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)) %>%
mutate(Bd_status=factor(Bd_status,levels=c("Bd_positive", "Bd_negative")))%>%
group_by(Bd_status) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = Bd_status)) +
geom_point(size = 3, alpha = 0.9) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.5, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(face = "bold", size = 12),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Bd status')
beta_q1f_mat <- as.matrix(beta_q1f$S)
beta_q1f_filtered <- beta_q1f_mat[rownames(beta_q1f_mat) %in% sample_metadata$sample,
colnames(beta_q1f_mat) %in% sample_metadata$sample]
beta_q1f_dis <- as.dist(beta_q1f_filtered)
sample_metadata_ordered <- sample_metadata %>%
arrange(match(sample,labels(beta_q1f_dis))) %>%
column_to_rownames(., "sample")
all(sample_metadata_ordered$sample == labels(beta_q1f_dis))[1] 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 1 0.10182 0.101818 2.5757 999 0.102
Residuals 37 1.46259 0.039529
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Bd_negative Bd_positive
Bd_negative 0.106
Bd_positive 0.11702 adonis2(beta_q1f_dis ~ Bd_status+region+sex,
data = sample_metadata_ordered,
permutations = 999, by="margin") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Bd_status | 1 | 0.09393142 | 0.01659706 | 0.7559309 | 0.450 |
| region | 1 | 0.88568078 | 0.15649389 | 7.1276838 | 0.007 |
| sex | 1 | 0.09744403 | 0.01721771 | 0.7841992 | 0.424 |
| Residual | 35 | 4.34907446 | 0.76845246 | NA | NA |
| Total | 38 | 5.65952309 | 1.00000000 | NA | NA |