Chapter 14 What are the differences at FMT1 between groups?
14.1 Is the GM of the intervention group more similar to warm population after one week from FMT?
14.1.1 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 =="FMT1") %>%
filter(metric=="richness") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.7046013alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point =="FMT1") %>%
filter(metric=="neutral") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.4100881alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point =="FMT1") %>%
filter(metric=="phylogenetic") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.932363814.1.2 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(method="t.test",size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_blank(),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
***p<0.05 when comparing neutral beta diversity between Cold-control and Warm-control
14.1.3 Beta diversity
Number of samples used
samples_to_keep_post1 <- sample_metadata %>%
filter(time_point == "FMT1") %>%
select(Tube_code) %>%
pull()
subset_meta_post1 <- sample_metadata %>%
filter(time_point == "FMT1")
length(samples_to_keep_post1)[1] 26Richness
richness_post1 <- as.matrix(beta_q0n$S)
richness_post1 <- as.dist(richness_post1[rownames(richness_post1) %in% samples_to_keep_post1,
colnames(richness_post1) %in% samples_to_keep_post1])
betadisper(richness_post1, subset_meta_post1$type) %>% 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 2 0.010965 0.0054824 1.1892 999 0.321
Residuals 23 0.106037 0.0046103
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.653000 0.044
Cold_intervention 0.611267 0.453
Warm_control 0.051747 0.423493 adonis2(richness_post1 ~ type,
data = subset_meta_post1 %>% arrange(match(Tube_code,labels(richness_post1))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.193445 | 0.1448598 | 1.948088 | 0.001 |
| Residual | 23 | 7.045172 | 0.8551402 | NA | NA |
| Total | 25 | 8.238617 | 1.0000000 | NA | NA |
subset_meta_post1_arrange <- column_to_rownames(subset_meta_post1, "Tube_code")
subset_meta_post1_arrange<-subset_meta_post1_arrange[labels(richness_post1),]
pairwise <- pairwise.adonis(richness_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_control vs Cold_intervention | 1 | 0.5608952 | 1.729368 | 0.10337316 | 0.010 | 0.030 | . |
| Cold_control vs Warm_control | 1 | 0.8405234 | 2.787341 | 0.14836272 | 0.001 | 0.003 | * |
| Cold_intervention vs Warm_control | 1 | 0.3744068 | 1.276240 | 0.07841121 | 0.107 | 0.321 |
Neutral
neutral_post1 <- as.matrix(beta_q1n$S)
neutral_post1 <- as.dist(neutral_post1[rownames(neutral_post1) %in% samples_to_keep_post1,
colnames(neutral_post1) %in% samples_to_keep_post1])
betadisper(neutral_post1, subset_meta_post1$type) %>% 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 2 0.003113 0.0015566 0.1424 999 0.882
Residuals 23 0.251460 0.0109331
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.87300 0.540
Cold_intervention 0.86225 0.781
Warm_control 0.54602 0.76361 adonis2(neutral_post1 ~ type,
data = subset_meta_post1 %>% arrange(match(Tube_code,labels(neutral_post1))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.395592 | 0.1901271 | 2.69976 | 0.001 |
| Residual | 23 | 5.944717 | 0.8098729 | NA | NA |
| Total | 25 | 7.340309 | 1.0000000 | NA | NA |
subset_meta_post1_arrange <- column_to_rownames(subset_meta_post1, "Tube_code")
subset_meta_post1_arrange<-subset_meta_post1_arrange[labels(neutral_post1),]
pairwise <- pairwise.adonis(neutral_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_control vs Cold_intervention | 1 | 0.6048786 | 2.250708 | 0.13047046 | 0.017 | 0.051 | |
| Cold_control vs Warm_control | 1 | 1.0518347 | 4.143449 | 0.20569710 | 0.001 | 0.003 | * |
| Cold_intervention vs Warm_control | 1 | 0.4158490 | 1.643022 | 0.09872138 | 0.057 | 0.171 |
Phylogenetic
phylo_post1 <- as.matrix(beta_q1p$S)
phylo_post1 <- as.dist(phylo_post1[rownames(phylo_post1) %in% samples_to_keep_post1,
colnames(phylo_post1) %in% samples_to_keep_post1])
betadisper(phylo_post1, subset_meta_post1$type) %>% 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 2 0.00398 0.001988 0.1209 999 0.913
Residuals 23 0.37833 0.016449
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.73000 0.917
Cold_intervention 0.65311 0.782
Warm_control 0.90970 0.74945 adonis2(phylo_post1 ~ type,
data = subset_meta_post1 %>% arrange(match(Tube_code,labels(phylo_post1))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 0.07537671 | 0.0713385 | 0.8834142 | 0.534 |
| Residual | 23 | 0.98122963 | 0.9286615 | NA | NA |
| Total | 25 | 1.05660634 | 1.0000000 | NA | NA |
subset_meta_post1_arrange <- column_to_rownames(subset_meta_post1, "Tube_code")
subset_meta_post1_arrange<-subset_meta_post1_arrange[labels(phylo_post1),]
pairwise <- pairwise.adonis(phylo_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_control vs Cold_intervention | 1 | 0.01874964 | 0.4216291 | 0.02734011 | 0.764 | 1.000 | |
| Cold_control vs Warm_control | 1 | 0.06078920 | 1.7647327 | 0.09933911 | 0.105 | 0.315 | |
| Cold_intervention vs Warm_control | 1 | 0.03216734 | 0.6483005 | 0.04142945 | 0.697 | 1.000 |
NMDS
richness_post1 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post1, by = join_by(sample == Tube_code)) %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c","#76b183")) +
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(legend.position="none")
neutral_post1 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post1, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c","#76b183")) +
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(legend.position="none")
phylo_post1 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post1, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c","#76b183")) +
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(legend.position="none")
14.1.3.1 Functional differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Cold_control 0.373 0.0247
2 Cold_intervention 0.353 0.0186
3 Warm_control 0.372 0.0370MCI_tm5_ele <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "FMT1")
shapiro.test(MCI_tm5_ele$value)#normality test
Shapiro-Wilk normality test
data: MCI_tm5_ele$value
W = 0.9246, p-value = 0.05771 Df Sum Sq Mean Sq F value Pr(>F)
type 2 0.002156 0.0010781 1.359 0.277
Residuals 23 0.018251 0.0007935 Functional differences
CC and CI
element_gift %>%
filter(type %in% c("Cold_control", "Cold_intervention") & time_point == "FMT1") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 6)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, type),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames()# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>CI and WC
element_gift %>%
filter(type %in% c("Warm_control", "Cold_intervention") &
time_point == "FMT1") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where(~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 6)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
remove_rownames()# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>CC and WC
element_gift %>%
filter(type %in% c("Warm_control", "Cold_control") &
time_point == "FMT1") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where(~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 6)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
remove_rownames()# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>14.2 Is the GM of the intervention group still more similar to warm population after one week from FMT?
14.2.1 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 =="FMT2") %>%
filter(metric=="richness") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.1206887alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point =="FMT2") %>%
filter(metric=="neutral") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.5856489alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point =="FMT2") %>%
filter(metric=="phylogenetic") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.350017114.2.2 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(method="t.test",size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_blank(),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
***p<0.05 when comparing richness and neutral beta diversity between Cold-control and Cold-intervention
14.2.3 Beta diversity
Number of samples used
samples_to_keep_FMT2 <- sample_metadata %>%
filter(time_point == "FMT2") %>%
select(Tube_code) %>%
pull()
subset_meta_FMT2 <- sample_metadata %>%
filter(time_point == "FMT2")
length(samples_to_keep_FMT2)[1] 27Richness
richness_FMT2 <- as.matrix(beta_q0n$S)
richness_FMT2 <- as.dist(richness_FMT2[rownames(richness_FMT2) %in% samples_to_keep_FMT2,
colnames(richness_FMT2) %in% samples_to_keep_FMT2])
betadisper(richness_FMT2, subset_meta_FMT2$type) %>% 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 2 0.00320 0.0016002 0.2688 999 0.787
Residuals 24 0.14286 0.0059526
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.51500 0.937
Cold_intervention 0.49519 0.564
Warm_control 0.93090 0.53472 adonis2(richness_FMT2 ~ type,
data = subset_meta_FMT2 %>% arrange(match(Tube_code,labels(richness_FMT2))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.516568 | 0.1987688 | 2.97695 | 0.001 |
| Residual | 24 | 6.113243 | 0.8012312 | NA | NA |
| Total | 26 | 7.629812 | 1.0000000 | NA | NA |
subset_meta_FMT2_arrange <- column_to_rownames(subset_meta_FMT2, "Tube_code")
subset_meta_FMT2_arrange<-subset_meta_FMT2_arrange[labels(richness_FMT2),]
pairwise <- pairwise.adonis(richness_FMT2, subset_meta_FMT2_arrange$type, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_control vs Cold_intervention | 1 | 0.6523627 | 2.595811 | 0.1395912 | 0.001 | 0.003 | * |
| Cold_control vs Warm_control | 1 | 1.1383473 | 4.313909 | 0.2123623 | 0.002 | 0.006 | * |
| Cold_intervention vs Warm_control | 1 | 0.4841426 | 1.944634 | 0.1083685 | 0.001 | 0.003 | * |
Neutral
neutral_FMT2 <- as.matrix(beta_q1n$S)
neutral_FMT2 <- as.dist(neutral_FMT2[rownames(neutral_FMT2) %in% samples_to_keep_FMT2,
colnames(neutral_FMT2) %in% samples_to_keep_FMT2])
betadisper(neutral_FMT2, subset_meta_FMT2$type) %>% 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 2 0.005409 0.0027044 0.3101 999 0.728
Residuals 24 0.209280 0.0087200
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.63400 0.726
Cold_intervention 0.64372 0.465
Warm_control 0.75385 0.46141 adonis2(neutral_FMT2 ~ type,
data = subset_meta_FMT2 %>% arrange(match(Tube_code,labels(neutral_FMT2))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.927421 | 0.2610114 | 4.23841 | 0.001 |
| Residual | 24 | 5.457010 | 0.7389886 | NA | NA |
| Total | 26 | 7.384430 | 1.0000000 | NA | NA |
subset_meta_FMT2_arrange <- column_to_rownames(subset_meta_FMT2, "Tube_code")
subset_meta_FMT2_arrange<-subset_meta_FMT2_arrange[labels(neutral_FMT2),]
pairwise <- pairwise.adonis(neutral_FMT2, subset_meta_FMT2_arrange$type, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_control vs Cold_intervention | 1 | 1.0248413 | 4.665114 | 0.2257483 | 0.002 | 0.006 | * |
| Cold_control vs Warm_control | 1 | 1.3640304 | 5.787044 | 0.2656186 | 0.001 | 0.003 | * |
| Cold_intervention vs Warm_control | 1 | 0.5022591 | 2.215131 | 0.1216094 | 0.001 | 0.003 | * |
Phylogenetic
phylo_FMT2 <- as.matrix(beta_q1p$S)
phylo_FMT2 <- as.dist(phylo_FMT2[rownames(phylo_FMT2) %in% samples_to_keep_FMT2,
colnames(phylo_FMT2) %in% samples_to_keep_FMT2])
betadisper(phylo_FMT2, subset_meta_FMT2$type) %>% 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 2 0.002914 0.0014571 0.3524 999 0.699
Residuals 24 0.099228 0.0041345
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.59600 0.784
Cold_intervention 0.61692 0.410
Warm_control 0.77655 0.41199 adonis2(phylo_FMT2 ~ type,
data = subset_meta_FMT2 %>% arrange(match(Tube_code,labels(phylo_FMT2))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 0.1596594 | 0.204565 | 3.086085 | 0.001 |
| Residual | 24 | 0.6208232 | 0.795435 | NA | NA |
| Total | 26 | 0.7804826 | 1.000000 | NA | NA |
subset_meta_FMT2_arrange <- column_to_rownames(subset_meta_FMT2, "Tube_code")
subset_meta_FMT2_arrange<-subset_meta_FMT2_arrange[labels(phylo_FMT2),]
pairwise <- pairwise.adonis(phylo_FMT2, subset_meta_FMT2_arrange$type, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_control vs Cold_intervention | 1 | 0.05942531 | 2.398475 | 0.1303627 | 0.031 | 0.093 | |
| Cold_control vs Warm_control | 1 | 0.11091887 | 4.047943 | 0.2019131 | 0.004 | 0.012 | . |
| Cold_intervention vs Warm_control | 1 | 0.06914494 | 2.719532 | 0.1452778 | 0.005 | 0.015 | . |
NMDS
richness_FMT2 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_FMT2, by = join_by(sample == Tube_code)) %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c","#76b183")) +
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(legend.position="none")
neutral_FMT2 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_FMT2, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c","#76b183")) +
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(legend.position="none")
phylo_FMT2 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_FMT2, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("#4477AA","#d57d2c","#76b183")) +
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(legend.position="none")
14.2.3.1 Functional differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2") %>%
group_by(type) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
type MCI sd
<chr> <dbl> <dbl>
1 Cold_control 0.352 0.0221
2 Cold_intervention 0.346 0.0255
3 Warm_control 0.350 0.0294MCI_FMT2 <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "FMT2")
shapiro.test(MCI_FMT2$value)
Shapiro-Wilk normality test
data: MCI_FMT2$value
W = 0.96179, p-value = 0.4057 Df Sum Sq Mean Sq F value Pr(>F)
type 2 0.000176 0.0000878 0.132 0.877
Residuals 24 0.016019 0.0006674 Functional differences
CC and CI
significant_elements<-element_gift %>%
filter(type %in% c("Cold_control", "Cold_intervention") & time_point == "FMT2") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 6)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, type),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))element_gift_sig <- element_gift %>%
select(Tube_code, all_of(intersect(
significant_elements$trait,
colnames(element_gift)
))) %>%
left_join(., sample_metadata[c(1, 6)], by = join_by(Tube_code == Tube_code)) %>%
filter(type!="Warm_control")
difference_table <- element_gift_sig %>%
select(-Tube_code) %>%
group_by(type) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Cold_control-Cold_intervention)%>%
mutate(group_color = ifelse(Difference <0, "Cold_control","Cold_intervention")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c("#4477AA","#76b183")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
legend.title = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Function") +
ylab("Mean difference")
CI and WC
element_gift %>%
filter(type %in% c("Warm_control", "Cold_intervention") &
time_point == "FMT2") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where(~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 6)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
remove_rownames()# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>CC and WC
element_gift %>%
filter(type %in% c("Warm_control", "Cold_control") &
time_point == "FMT2") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where(~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 6)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ type)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
remove_rownames()# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>