Chapter 7 Beta diversity
load("data/data.Rdata")
behaviour_colors <- c(
tame = "#1f77b4",
aggr = "#d62728",
unsel = "#2ca02c"
)library(hilldiv2)
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0) #hill pair (beta divesity)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
genome_counts_filt_beta <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0))%>%
rownames_to_column(., "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt_beta$genome)
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)
genome_gifts1 <- genome_gifts[genome_counts_filt_beta$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
dist <- genome_gifts1 %>%
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)7.0.1 Permanova
7.0.1.1 Richness
The permanova test assumes within-group dispersions are equal, so before doing this test, we must check the dispersion (variances) of each of the variables.
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 2 0.05431 0.027155 2.3627 999 0.107
Residuals 43 0.49420 0.011493
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 3 1.0776 0.35921 14.463 999 0.001 ***
Residuals 42 1.0432 0.02484
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Richness
adonis2(beta_q0n$S ~ fox_behaviour+gut_location,
by="terms",
data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| fox_behaviour | 2 | 0.7825353 | 0.04917873 | 1.739559 | 0.062 |
| gut_location | 3 | 6.1325920 | 0.38540508 | 9.088411 | 0.001 |
| Residual | 40 | 8.9969411 | 0.56541619 | NA | NA |
| Total | 45 | 15.9120684 | 1.00000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 tame vs aggr 1 0.3562095 0.9847508 0.02662586 0.375 1.000
2 tame vs unsel 1 0.5551591 1.5945065 0.05778348 0.126 0.378
3 aggr vs unsel 1 0.2660081 0.7800316 0.03147823 0.503 1.000 pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 K_ileum vs D_ileum 1 2.0999028 8.728723 0.2936125 0.001 0.006 *
2 K_ileum vs F_colon 1 4.4051839 46.802248 0.7006089 0.001 0.006 *
3 K_ileum vs M_colon 1 1.5958589 6.863418 0.2463236 0.001 0.006 *
4 D_ileum vs F_colon 1 1.6943456 7.242188 0.2564315 0.001 0.006 *
5 D_ileum vs M_colon 1 0.9359521 2.602023 0.1057646 0.010 0.060
6 F_colon vs M_colon 1 1.6519235 7.312706 0.2582835 0.001 0.006 *7.0.1.2 Neutral
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 2 0.07802 0.039012 0.7095 999 0.496
Residuals 43 2.36440 0.054986
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 3 1.53752 0.51251 26.845 999 0.001 ***
Residuals 42 0.80183 0.01909
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Neutral diversity
adonis2(beta_q1n$S ~ fox_behaviour +gut_location,
by="terms",
data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| fox_behaviour | 2 | 0.5878185 | 0.05087496 | 2.329149 | 0.065 |
| gut_location | 3 | 5.9188680 | 0.51227064 | 15.635129 | 0.001 |
| Residual | 40 | 5.0474952 | 0.43685441 | NA | NA |
| Total | 45 | 11.5541817 | 1.00000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 tame vs aggr 1 0.2903289 1.0883705 0.02934533 0.303 0.909
2 tame vs unsel 1 0.4021091 1.6575473 0.05993110 0.168 0.504
3 aggr vs unsel 1 0.1826192 0.7277943 0.02943224 0.485 1.000 pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 K_ileum vs D_ileum 1 1.8947421 11.029882 0.3443622 0.001 0.006 *
2 K_ileum vs F_colon 1 4.5677089 140.238751 0.8751862 0.001 0.006 *
3 K_ileum vs M_colon 1 0.4938235 7.424449 0.2611994 0.004 0.024 .
4 D_ileum vs F_colon 1 1.5545564 7.689604 0.2680276 0.001 0.006 *
5 D_ileum vs M_colon 1 1.2010433 5.294345 0.1939722 0.001 0.006 *
6 F_colon vs M_colon 1 2.2584522 23.308509 0.5260504 0.001 0.006 *7.0.1.3 Phylogenetic
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 2 0.06047 0.030237 0.9492 999 0.398
Residuals 43 1.36981 0.031856
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 3 0.95298 0.31766 18.558 999 0.001 ***
Residuals 42 0.71892 0.01712
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Phylogenetic diversity
adonis2(beta_q1p$S ~ fox_behaviour+ gut_location,
by="terms",
data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| fox_behaviour | 2 | 0.3794287 | 0.05178089 | 2.492071 | 0.055 |
| gut_location | 3 | 3.9030657 | 0.53265400 | 17.090110 | 0.001 |
| Residual | 40 | 3.0450873 | 0.41556511 | NA | NA |
| Total | 45 | 7.3275817 | 1.00000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 tame vs aggr 1 0.16506690 0.9826933 0.02657171 0.367 1.000
2 tame vs unsel 1 0.30829662 2.0858230 0.07426605 0.112 0.336
3 aggr vs unsel 1 0.09803357 0.5872774 0.02388542 0.582 1.000 pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 K_ileum vs D_ileum 1 1.207405 11.239217 0.3486194 0.001 0.006 *
2 K_ileum vs F_colon 1 3.040599 281.693909 0.9337076 0.001 0.006 *
3 K_ileum vs M_colon 1 0.446189 9.615557 0.3140742 0.002 0.012 .
4 D_ileum vs F_colon 1 1.114288 9.540696 0.3123929 0.002 0.012 .
5 D_ileum vs M_colon 1 0.791897 5.425242 0.1978193 0.003 0.018 .
6 F_colon vs M_colon 1 1.298322 23.280729 0.5257531 0.001 0.006 *7.0.1.4 Functional
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 2 0.02976 0.014878 0.3215 999 0.735
Residuals 43 1.99003 0.046280
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 3 0.41181 0.13727 4.1123 999 0.009 **
Residuals 42 1.40196 0.03338
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Functional diversity
adonis2(beta_q1f$S ~ fox_behaviour+ gut_location,
by="terms",
data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| fox_behaviour | 2 | 0.1378307 | 0.01798362 | 1.066254 | 0.344 |
| gut_location | 3 | 4.9410817 | 0.64469315 | 25.482706 | 0.001 |
| Residual | 40 | 2.5853255 | 0.33732323 | NA | NA |
| Total | 45 | 7.6642379 | 1.00000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 tame vs aggr 1 0.06764083 0.3798683 0.01044172 0.632 1
2 tame vs unsel 1 0.07835157 0.4880544 0.01842545 0.561 1
3 aggr vs unsel 1 0.06061608 0.3255638 0.01338361 0.614 1 pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 K_ileum vs D_ileum 1 1.4924134 21.729688 0.5085384 0.001 0.006 *
2 K_ileum vs F_colon 1 3.9491603 204.127340 0.9107650 0.001 0.006 *
3 K_ileum vs M_colon 1 0.1228807 2.490336 0.1060154 0.130 0.780
4 D_ileum vs F_colon 1 0.9256356 11.421302 0.3522777 0.003 0.018 .
5 D_ileum vs M_colon 1 0.7091245 6.635195 0.2317147 0.008 0.048 .
6 F_colon vs M_colon 1 2.7938493 45.276168 0.6831440 0.001 0.006 *7.0.2 NMDS: Non-metric Multidimensional Scaling
7.0.2.1 Richness
richness_behaviour <-beta_q0n$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>% #2D ordination
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(fox_behaviour) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>% #centroids
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = fox_behaviour, shape = as.factor(gut_location))) +
geom_point(size = 4) +
scale_color_manual(values = c(
tame = "#1f77b4",
aggr = "#d62728",
unsel = "#2ca02c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(color="Fox behaviour", shape="Gut Location")+geom_text_repel(aes(label = sample), size=3)
richness_behaviour
7.0.2.2 Richness between gut locations
richness_gutlocation <- 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(gut_location) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = gut_location, shape = as.factor(fox_behaviour))) +
geom_point(size = 4) +
scale_color_manual(values = c(
F_colon = "#d6604d",
M_colon = "#542788",
D_ileum = "#fdb863",
K_ileum = "#bb99d8")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(title = "Richness",
color="Gut location",
shape="Fox Behaviour")
#+geom_text_repel(aes(label = sample), size=3)
richness_gutlocation
7.0.2.3 Neutral diversity
neutral_behaviour <- 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(fox_behaviour) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = fox_behaviour, shape = as.factor(gut_location))) +
geom_point(size = 4) +
scale_color_manual(values = behaviour_colors,labels=c("tame" = "Tame", "aggr" = "Aggressive", "unsel" = "Unselected")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(color="Fox Behaviour", shape="Gut Location")
neutral_behaviour
####Neutral diversity by gut location
neutral_gutlocation <- 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(gut_location) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = gut_location, shape = as.factor(fox_behaviour))) +
geom_point(size = 4) +
scale_color_manual(values = c(
F_colon = "#d6604d",
M_colon = "#542788",
D_ileum = "#fdb863",
K_ileum = "#bb99d8")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(title = "Neutral",color="Gut location", shape="Fox Behaviour")
neutral_gutlocation
7.0.2.4 Phylogenetic diversity
phylogenetic_behaviour <-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)) %>%
group_by(fox_behaviour) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = fox_behaviour, shape = as.factor(gut_location))) +
geom_point(size = 4) +
scale_color_manual(values = behaviour_colors,labels=c("tame" = "Tame", "aggr" = "Aggressive", "unsel" = "Unselected")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(title = "Phylogenetic",color="Fox Behaviour", shape="Gut Location")+geom_text_repel(aes(label = sample), size=3)
phylogenetic_behaviour 
7.0.2.5 Phylogenetic diversity by gut location
phylogenetic_gutlocation <- 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)) %>%
group_by(gut_location) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = gut_location, shape = as.factor(fox_behaviour))) +
geom_point(size = 4) +
scale_color_manual(values = c(
F_colon = "#d6604d",
M_colon = "#542788",
D_ileum = "#fdb863",
K_ileum = "#bb99d8")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(title = "Phylogenetic",color="Gut location", shape="Fox Behaviour")
#+geom_text_repel(aes(label = sample), size=3)
phylogenetic_gutlocation
7.0.2.6 Functional diversity
functional_behaviour <- 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)) %>%
group_by(fox_behaviour) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = fox_behaviour, shape = as.factor(gut_location))) +
geom_point(size = 4) +
scale_color_manual(values = behaviour_colors,labels=c("tame" = "Tame", "aggr" = "Aggressive", "unsel" = "Unselected")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(color="Fox Behaviour", shape="Gut Location")+geom_text_repel(aes(label = sample), size=3)
functional_behaviour
7.0.2.7 Functional diversity by gut location
functional_gutlocation <- 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)) %>%
group_by(gut_location) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = gut_location, shape = as.factor(fox_behaviour))) +
geom_point(size = 4) +
scale_color_manual(values = c(
F_colon = "#d6604d",
M_colon = "#542788",
D_ileum = "#fdb863",
K_ileum = "#bb99d8")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
) +
labs(title = "Functional",color="Gut location", shape="Fox Behaviour")
#+geom_text_repel(aes(label = sample), size=3)
functional_gutlocation
library(patchwork)
(neutral_gutlocation + phylogenetic_gutlocation + phylogenetic_gutlocation + functional_gutlocation)+
plot_layout(
guides = "collect",
ncol = 2
) 