Chapter 13 Faecal Microbiota Transplantation (FMT)
13.1 Does FMT change the microbial community over time?
13.1.1 Alpha diversity
sample_metadata <- sample_metadata %>%
mutate(treatment=time_point) %>%
mutate(treatment=case_when(
treatment %in% c("Inoculum1", "Inoculum2") ~ "Inoculum",
TRUE ~ treatment
))
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Cold_intervention" & treatment %in% c("FMT1","FMT2") ) 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", "phylogenetic"))) %>%
filter(type=="Cold_intervention" & treatment %in% c("FMT1", "FMT2")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
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("#76b183", '#40714b50')) +
scale_fill_manual(values=c("#76b183", '#40714b50')) +
facet_wrap( ~ metric, scales = "free")+
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.y = element_text(size = 14),
axis.title.x = element_blank(),
strip.text = element_text(size = 16, lineheight = 0.6)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.3876) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik -2*log(L) df.resid
171.0 174.3 -81.5 163.0 13
Scaled residuals:
Min 1Q Median 3Q Max
-1.73495 -0.71265 -0.07086 0.40744 1.88240
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 1.073e-11 3.275e-06
Number of obs: 17, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9343 0.1759 22.369 <2e-16 ***
treatmentFMT2 0.4052 0.2402 1.687 0.0917 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tretmntFMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
128.7946 131.6268 -60.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.497472 11.25384
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 24.65947 4.164756 8 5.920988 0.0004
treatmentFMT2 14.87494 5.482531 7 2.713151 0.0301
Correlation:
(Intr)
treatmentFMT2 -0.7
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.76462323 -0.55701822 -0.04763166 0.55267588 1.30333436
Number of Observations: 17
Number of Groups: 9 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
56.17121 59.00341 -24.0856
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.6979711 0.8383981
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 4.163606 0.3820859 8 10.897043 0.0000
treatmentFMT2 0.916226 0.4122640 7 2.222426 0.0617
Correlation:
(Intr)
treatmentFMT2 -0.583
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.1362905 -0.5779277 -0.1466831 0.4812394 2.3357705
Number of Observations: 17
Number of Groups: 9 13.1.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("FMT1", "FMT2")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("FMT1", "FMT2"))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$treatment) %>% 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.017610 0.017610 2.8225 999 0.094 .
Residuals 15 0.093585 0.006239
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3560084 | 0.07968734 | 1.298809 | 0.00390625 |
| Residual | 15 | 4.1115570 | 0.92031266 | NA | NA |
| Total | 16 | 4.4675655 | 1.00000000 | 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$treatment) %>% 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.009828 0.0098278 0.7113 999 0.42
Residuals 15 0.207263 0.0138175 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3149954 | 0.08110541 | 1.323962 | 0.00390625 |
| Residual | 15 | 3.5687823 | 0.91889459 | NA | NA |
| Total | 16 | 3.8837776 | 1.00000000 | 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$treatment) %>% 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.010955 0.010955 0.6602 999 0.534
Residuals 15 0.248892 0.016593 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06391536 | 0.09449338 | 1.565312 | 0.0703125 |
| Residual | 15 | 0.61248501 | 0.90550662 | NA | NA |
| Total | 16 | 0.67640037 | 1.00000000 | NA | NA |
NMDS
richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="Time point",values=c("#76b183", '#40714b50')) +
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"
)
neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="Time point",values=c("#76b183", '#40714b50')) +
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"
)
phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="Time point",values=c("#76b183", '#40714b50')) +
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"
)
13.2 Do FMT change the microbial community compared to antibiotics and acclimation?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Cold_intervention" & treatment %in% c("Antibiotics", "Acclimation","FMT1","FMT2") ) 13.2.1 Alpha diversity
Richness
# Modelq0GLMMNB <- glmer.nb(richness ~ treatment + (1|individual), data = alpha_div_meta)
# summary(Modelq0GLMMNB)
#emmeans(Modelq0GLMMNB, pairwise ~ treatment)
model_nb <- MASS::glm.nb(richness ~ treatment, data = alpha_div_meta)
summary(model_nb)
Call:
MASS::glm.nb(formula = richness ~ treatment, data = alpha_div_meta,
init.theta = 2.622371751, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.8067 0.2118 17.977 < 2e-16 ***
treatmentAntibiotics -0.9651 0.3281 -2.941 0.00327 **
treatmentFMT1 0.1276 0.3081 0.414 0.67878
treatmentFMT2 0.5328 0.2978 1.789 0.07355 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.6224) family taken to be 1)
Null deviance: 54.752 on 32 degrees of freedom
Residual deviance: 35.688 on 29 degrees of freedom
AIC: 314.49
Number of Fisher Scoring iterations: 1
Theta: 2.622
Std. Err.: 0.685
2 x log-likelihood: -304.492 $emmeans
treatment emmean SE df asymp.LCL asymp.UCL
Acclimation 3.81 0.212 Inf 3.39 4.22
Antibiotics 2.84 0.251 Inf 2.35 3.33
FMT1 3.93 0.224 Inf 3.50 4.37
FMT2 4.34 0.209 Inf 3.93 4.75
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - Antibiotics 0.965 0.328 Inf 2.941 0.0172
Acclimation - FMT1 -0.128 0.308 Inf -0.414 0.9761
Acclimation - FMT2 -0.533 0.298 Inf -1.789 0.2784
Antibiotics - FMT1 -1.093 0.336 Inf -3.252 0.0063
Antibiotics - FMT2 -1.498 0.327 Inf -4.587 <.0001
FMT1 - FMT2 -0.405 0.306 Inf -1.322 0.5488
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 4 estimates Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
235.3449 243.5487 -111.6725
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.485828 9.015585
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.310151 3.356642 21 5.156985 0.0000
treatmentAntibiotics -8.962870 4.586486 21 -1.954191 0.0641
treatmentFMT1 7.203346 4.402070 21 1.636354 0.1167
treatmentFMT2 22.224251 4.249987 21 5.229251 0.0000
Correlation:
(Intr) trtmnA trFMT1
treatmentAntibiotics -0.587
treatmentFMT1 -0.611 0.457
treatmentFMT2 -0.633 0.463 0.483
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.17355918 -0.51312522 0.04290664 0.35530718 1.54377193
Number of Observations: 33
Number of Groups: 9 $emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 17.31 3.36 8 9.570 25.1
Antibiotics 8.35 3.77 8 -0.355 17.0
FMT1 24.51 3.55 8 16.334 32.7
FMT2 39.53 3.36 8 31.794 47.3
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 8.96 4.59 21 1.954 0.2367
Acclimation - FMT1 -7.20 4.40 21 -1.636 0.3809
Acclimation - FMT2 -22.22 4.25 21 -5.229 0.0002
Antibiotics - FMT1 -16.17 4.69 21 -3.448 0.0119
Antibiotics - FMT2 -31.19 4.59 21 -6.800 <.0001
FMT1 - FMT2 -15.02 4.40 21 -3.412 0.0129
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
235.3449 243.5487 -111.6725
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.485828 9.015585
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.310151 3.356642 21 5.156985 0.0000
treatmentAntibiotics -8.962870 4.586486 21 -1.954191 0.0641
treatmentFMT1 7.203346 4.402070 21 1.636354 0.1167
treatmentFMT2 22.224251 4.249987 21 5.229251 0.0000
Correlation:
(Intr) trtmnA trFMT1
treatmentAntibiotics -0.587
treatmentFMT1 -0.611 0.457
treatmentFMT2 -0.633 0.463 0.483
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.17355918 -0.51312522 0.04290664 0.35530718 1.54377193
Number of Observations: 33
Number of Groups: 9 $emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 17.31 3.36 8 9.570 25.1
Antibiotics 8.35 3.77 8 -0.355 17.0
FMT1 24.51 3.55 8 16.334 32.7
FMT2 39.53 3.36 8 31.794 47.3
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 8.96 4.59 21 1.954 0.2367
Acclimation - FMT1 -7.20 4.40 21 -1.636 0.3809
Acclimation - FMT2 -22.22 4.25 21 -5.229 0.0002
Antibiotics - FMT1 -16.17 4.69 21 -3.448 0.0119
Antibiotics - FMT2 -31.19 4.59 21 -6.800 <.0001
FMT1 - FMT2 -15.02 4.40 21 -3.412 0.0129
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
123.367 131.5707 -55.68348
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.6707465 1.301744
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 5.503182 0.4881300 21 11.274009 0.0000
treatmentAntibiotics -1.243682 0.6625044 21 -1.877243 0.0744
treatmentFMT1 -1.370236 0.6357519 21 -2.155300 0.0429
treatmentFMT2 -0.423350 0.6136480 21 -0.689891 0.4978
Correlation:
(Intr) trtmnA trFMT1
treatmentAntibiotics -0.582
treatmentFMT1 -0.607 0.457
treatmentFMT2 -0.629 0.463 0.483
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.09114445 -0.42055042 -0.08791605 0.37048911 2.04753704
Number of Observations: 33
Number of Groups: 9 $emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 5.50 0.488 8 4.38 6.63
Antibiotics 4.26 0.548 8 3.00 5.52
FMT1 4.13 0.516 8 2.94 5.32
FMT2 5.08 0.488 8 3.95 6.21
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 1.244 0.663 21 1.877 0.2675
Acclimation - FMT1 1.370 0.636 21 2.155 0.1686
Acclimation - FMT2 0.423 0.614 21 0.690 0.8998
Antibiotics - FMT1 0.127 0.677 21 0.187 0.9976
Antibiotics - FMT2 -0.820 0.663 21 -1.238 0.6104
FMT1 - FMT2 -0.947 0.636 21 -1.489 0.4612
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates 13.2.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("Acclimation", "Antibiotics", "FMT1", "FMT2")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("Acclimation", "Antibiotics", "FMT1", "FMT2"))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$treatment) %>% 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 3 0.068214 0.0227380 3.5642 999 0.022 *
Residuals 29 0.185007 0.0063796
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 3 | 2.646967 | 0.2239942 | 2.790285 | 0.001 |
| Residual | 29 | 9.170156 | 0.7760058 | NA | NA |
| Total | 32 | 11.817124 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8287871 2.293722 0.14077335 0.006 0.036 .
2 Acclimation vs FMT1 1 0.9572484 2.939042 0.16383497 0.001 0.006 *
3 Acclimation vs FMT2 1 0.9209464 3.233991 0.16813938 0.001 0.006 *
4 Antibiotics vs FMT1 1 0.9764237 2.751191 0.17466559 0.002 0.012 .
5 Antibiotics vs FMT2 1 1.2837856 4.194747 0.23054713 0.002 0.012 .
6 FMT1 vs FMT2 1 0.3560084 1.298809 0.07968734 0.136 0.816 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$treatment) %>% 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 3 0.06048 0.020160 1.4603 999 0.247
Residuals 29 0.40036 0.013806 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 3 | 3.357422 | 0.29403 | 4.026078 | 0.001 |
| Residual | 29 | 8.061213 | 0.70597 | NA | NA |
| Total | 32 | 11.418635 | 1.00000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8814932 2.747044 0.16403157 0.002 0.012 .
2 Acclimation vs FMT1 1 1.3133104 4.634354 0.23603291 0.001 0.006 *
3 Acclimation vs FMT2 1 1.3568927 5.355947 0.25079416 0.001 0.006 *
4 Antibiotics vs FMT1 1 1.3427497 4.355528 0.25095913 0.001 0.006 *
5 Antibiotics vs FMT2 1 1.5277817 5.613270 0.28619756 0.001 0.006 *
6 FMT1 vs FMT2 1 0.3149954 1.323962 0.08110541 0.152 0.912 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$treatment) %>% 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 3 0.26778 0.089261 3.9211 999 0.011 *
Residuals 29 0.66016 0.022764
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 3 | 1.821048 | 0.4220637 | 7.059515 | 0.001 |
| Residual | 29 | 2.493580 | 0.5779363 | NA | NA |
| Total | 32 | 4.314628 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.68859259 5.124832 0.26796742 0.006 0.036 .
2 Acclimation vs FMT1 1 0.25480266 3.292748 0.18000292 0.047 0.282
3 Acclimation vs FMT2 1 0.36626555 6.418219 0.28629479 0.003 0.018 .
4 Antibiotics vs FMT1 1 1.10004729 9.048069 0.41037921 0.004 0.024 .
5 Antibiotics vs FMT2 1 1.28533838 13.501094 0.49092934 0.001 0.006 *
6 FMT1 vs FMT2 1 0.06391536 1.565312 0.09449338 0.152 0.912 13.3 Is the gut microbiota similar to the inoculum after FMT?
13.3.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Cold_intervention" & treatment %in% c("Acclimation", "Antibiotics", "Inoculum", "FMT1", "FMT2") )
alpha_div_meta$treatment <- factor(alpha_div_meta$treatment, levels=c("Acclimation","Antibiotics","Inoculum", "FMT1", "FMT2"))
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","phylogenetic"))) %>%
filter(type=="Cold_intervention" & treatment %in% c("Acclimation", "Antibiotics", "Inoculum", "FMT1", "FMT2")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(name="Time point",values=c('#008080',"#003a45", "#d5b52c", "#76b183",'#40714b50')) +
scale_fill_manual(name="Time point",values=c('#008080',"#003a45", "#d5b52c", "#76b183",'#40714b50')) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text = element_text(size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.title.y = element_text(size = 14),
axis.title.x = element_blank(),
strip.text = element_text(size = 16, lineheight = 0.6)
) +
ylab("Alpha diversity")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Cold_intervention" & treatment %in% c("Inoculum", "FMT1", "FMT2") )Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(6.4056) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik -2*log(L) df.resid
301.6 308.6 -145.8 291.6 25
Scaled residuals:
Min 1Q Median 3Q Max
-2.05920 -0.54162 -0.04501 0.44663 2.23418
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 1.905e-12 1.38e-06
Number of obs: 30, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9343 0.1482 26.549 < 2e-16 ***
treatmentFMT2 0.4052 0.2019 2.007 0.044737 *
treatmentInoculum 0.7123 0.1863 3.824 0.000131 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trFMT2
tretmntFMT2 -0.734
trtmntInclm -0.795 0.584
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
230.8255 237.3046 -110.4127
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 5.858906 11.65548
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 24.76620 4.590358 19 5.395267 0.0000
treatmentFMT2 14.76820 5.687862 19 2.596441 0.0177
treatmentInoculum 21.54734 5.367099 19 4.014708 0.0007
Correlation:
(Intr) trFMT2
treatmentFMT2 -0.661
treatmentInoculum -0.706 0.570
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.19432220 -0.45597977 -0.05252155 0.50896675 1.77735772
Number of Observations: 30
Number of Groups: 9 Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
88.86545 95.34463 -39.43272
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.8599767 0.7111049
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 4.107135 0.3838878 19 10.698791 0.0000
treatmentFMT2 0.972697 0.3483993 19 2.791904 0.0116
treatmentInoculum 1.771914 0.3366142 19 5.263932 0.0000
Correlation:
(Intr) trFMT2
treatmentFMT2 -0.487
treatmentInoculum -0.514 0.566
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.58729803 -0.56622082 -0.06324004 0.47595562 2.31993950
Number of Observations: 30
Number of Groups: 9 13.3.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("Inoculum", "FMT1", "FMT2") ) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("Inoculum", "FMT1", "FMT2") )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$treatment) %>% 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 2 0.096689 0.048345 6.2474 999 0.008 **
Residuals 27 0.208934 0.007738
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.540586 | 0.2013816 | 3.404193 | 0.001 |
| Residual | 27 | 6.109499 | 0.7986184 | NA | NA |
| Total | 29 | 7.650086 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Inoculum vs FMT1 1 1.0476254 4.718687 0.19894384 0.001 0.003 *
2 Inoculum vs FMT2 1 0.8256961 4.246172 0.17512753 0.001 0.003 *
3 FMT1 vs FMT2 1 0.3560084 1.298809 0.07968734 0.128 0.384 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$treatment) %>% 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 2 0.03454 0.017270 1.2605 999 0.344
Residuals 27 0.36992 0.013701 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.462072 | 0.2050566 | 3.482342 | 0.001 |
| Residual | 27 | 5.668017 | 0.7949434 | NA | NA |
| Total | 29 | 7.130089 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Inoculum vs FMT1 1 1.0059569 4.799565 0.20166606 0.001 0.003 *
2 Inoculum vs FMT2 1 0.7900975 4.174913 0.17269608 0.001 0.003 *
3 FMT1 vs FMT2 1 0.3149954 1.323962 0.08110541 0.153 0.459 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$treatment) %>% 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 2 0.01120 0.0056009 0.3676 999 0.705
Residuals 27 0.41137 0.0152361 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 0.193734 | 0.1520363 | 2.420493 | 0.004 |
| Residual | 27 | 1.080528 | 0.8479637 | NA | NA |
| Total | 29 | 1.274262 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Inoculum vs FMT1 1 0.11151373 2.359092 0.11044910 0.037 0.111
2 Inoculum vs FMT2 1 0.10834872 3.331524 0.14279069 0.013 0.039 .
3 FMT1 vs FMT2 1 0.06391536 1.565312 0.09449338 0.140 0.420 NMDS
samples_to_keep <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("Acclimation","Antibiotics","Inoculum", "FMT1", "FMT2") ) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("Acclimation","Antibiotics","Inoculum", "FMT1", "FMT2") )
alpha_div_meta$treatment <- factor(alpha_div_meta$treatment, levels=c("Acclimation","Antibiotics","Inoculum", "FMT1", "FMT2"))
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 3) +
scale_color_manual(name="Time point",values=c('#008080',"#003a45", "#d5b52c", "#76b183",'#40714b50')) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.5, 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"
)
neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 3) +
scale_color_manual(name="Time point",values=c('#008080',"#003a45", "#d5b52c", "#76b183",'#40714b50')) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.5, 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"
)
phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 3) +
scale_color_manual(name="Time point",values=c('#008080',"#003a45", "#d5b52c", "#76b183",'#40714b50'))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.5, 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"
)
13.3.2.1 Differences in the functional capacities
sample_sub <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("FMT1","FMT2", "Inoculum"))
genome_counts_filtered <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample)))
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
genome_counts_filtered_row <- genome_counts_filtered %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 3 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT1 0.353 0.0186
2 FMT2 0.346 0.0255
3 Inoculum 0.354 0.0375MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.87208, p-value = 0.001864
Kruskal-Wallis rank sum test
data: value by treatment
Kruskal-Wallis chi-squared = 0.28797, df = 2, p-value = 0.865913.3.3 Inoculum vs FMT1
13.3.3.1 Differences in the functional capacities
sample_sub <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("FMT1", "Inoculum"))
genome_counts_filtered <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample)))
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
genome_counts_filtered_row <- genome_counts_filtered %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)13.3.3.1.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT1 0.359 0.0214
2 Inoculum 0.351 0.0426MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.93307, p-value = 0.1586
Wilcoxon rank sum exact test
data: value by treatment
W = 65, p-value = 0.3738
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(10,11)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(10,11)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
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=FMT1-Inoculum)%>%
mutate(group_color = ifelse(Difference <0, "Inoculum","FMT1")) 13.3.3.1.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT1 0.353 0.0186
2 Inoculum 0.354 0.0375MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.83139, p-value = 0.002064
Wilcoxon rank sum exact test
data: value by treatment
W = 59, p-value = 0.6452
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(10,11)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
significant_functional# A tibble: 2 × 4
trait p_value p_adjust Function
<chr> <dbl> <dbl> <chr>
1 B04 0.00341 0.0454 SCFA biosynthesis
2 B10 0.00454 0.0454 Antibiotic biosynthesis13.3.4 Inoculum vs FMT2
13.3.4.1 Differences in the functional capacities
sample_sub <- sample_metadata %>%
filter(type == "Cold_intervention" & treatment %in% c("FMT2", "Inoculum"))
genome_counts_filtered <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample)))
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
genome_counts_filtered_row <- genome_counts_filtered %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)13.3.4.1.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT2 0.350 0.0293
2 Inoculum 0.351 0.0426MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94558, p-value = 0.2575
Wilcoxon rank sum exact test
data: value by treatment
W = 63, p-value = 0.7938
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(10,11)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(10,11)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
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=FMT2-Inoculum)%>%
mutate(group_color = ifelse(Difference <0, "Inoculum","FMT2")) 13.3.4.1.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT2 0.346 0.0255
2 Inoculum 0.354 0.0375MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.86728, p-value = 0.00697
Wilcoxon rank sum exact test
data: value by treatment
W = 61, p-value = 0.896
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(10,11)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
significant_functional# A tibble: 1 × 4
trait p_value p_adjust Function
<chr> <dbl> <dbl> <chr>
1 B04 0.00193 0.0387 SCFA biosynthesis13.4 What is the trend of the microbiota in each type after FMT?
13.4.1 Alpha diversity
label_map <- c(
"Control" = "Cold-control",
"Treatment" = "Cold-intervention",
"Hot_control" = "Warm-control",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic 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 %in% c("FMT1","Acclimation", "FMT2")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
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")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2") Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.839) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik -2*log(L) df.resid
755.3 781.4 -366.6 733.3 68
Scaled residuals:
Min 1Q Median 3Q Max
-2.01947 -0.49054 0.05733 0.57010 1.94355
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.02455 0.1567
Number of obs: 79, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.8960 0.1684 23.132 <2e-16 ***
typeCold_intervention -0.1165 0.2281 -0.511 0.6096
typeWarm_control 0.5571 0.2266 2.459 0.0139 *
time_pointFMT1 -0.3001 0.2156 -1.392 0.1639
time_pointFMT2 0.1528 0.2161 0.707 0.4795
typeCold_intervention:time_pointFMT1 0.4174 0.3052 1.368 0.1714
typeWarm_control:time_pointFMT1 -0.1276 0.2967 -0.430 0.6673
typeCold_intervention:time_pointFMT2 0.4015 0.2974 1.350 0.1769
typeWarm_control:time_pointFMT2 -0.3927 0.2972 -1.321 0.1864
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typCl_ typWr_ t_FMT1 t_FMT2 tC_:_FMT1 tW_:_FMT1 tC_:_FMT2
typCld_ntrv -0.721
typWrm_cntr -0.742 0.535
tim_pntFMT1 -0.678 0.498 0.504
tim_pntFMT2 -0.703 0.508 0.522 0.532
typC_:_FMT1 0.483 -0.666 -0.358 -0.707 -0.378
typW_:_FMT1 0.497 -0.362 -0.672 -0.727 -0.389 0.515
typC_:_FMT2 0.499 -0.688 -0.371 -0.384 -0.719 0.511 0.280
typW_:_FMT2 0.516 -0.371 -0.685 -0.388 -0.730 0.276 0.516 0.524 $emmeans
type emmean SE df asymp.LCL asymp.UCL
Cold_control 3.85 0.1040 Inf 3.64 4.05
Cold_intervention 4.00 0.1020 Inf 3.80 4.20
Warm_control 4.23 0.0985 Inf 4.04 4.42
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Cold_control - Cold_intervention -0.157 0.143 Inf -1.093 0.5184
Cold_control - Warm_control -0.384 0.142 Inf -2.702 0.0189
Cold_intervention - Warm_control -0.227 0.141 Inf -1.609 0.2419
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates $emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.04 0.0931 Inf 3.86 4.23
FMT1 3.84 0.0942 Inf 3.65 4.02
FMT2 4.20 0.0888 Inf 4.02 4.37
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - FMT1 0.204 0.122 Inf 1.665 0.2187
Acclimation - FMT2 -0.156 0.121 Inf -1.292 0.3998
FMT1 - FMT2 -0.359 0.121 Inf -2.966 0.0085
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates Neutral
Modelq1n <- nlme::lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Modelq1n)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
560.9723 585.7058 -269.4862
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.411316 9.033337
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 21.069668 3.541939 46 5.948626 0.0000
typeCold_intervention -3.759517 4.875891 24 -0.771042 0.4482
typeWarm_control 23.020914 4.875891 24 4.721375 0.0001
time_pointFMT1 -3.710335 4.410208 46 -0.841306 0.4045
time_pointFMT2 2.453665 4.410208 46 0.556360 0.5807
typeCold_intervention:time_pointFMT1 10.918511 6.236977 46 1.750610 0.0867
typeWarm_control:time_pointFMT1 -14.895231 6.130541 46 -2.429676 0.0191
typeCold_intervention:time_pointFMT2 19.770586 6.130541 46 3.224933 0.0023
typeWarm_control:time_pointFMT2 -14.705214 6.130541 46 -2.398681 0.0206
Correlation:
(Intr) typCl_ typWr_ t_FMT1 t_FMT2 tC_:_FMT1 tW_:_FMT1 tC_:_FMT2
typeCold_intervention -0.726
typeWarm_control -0.726 0.528
time_pointFMT1 -0.665 0.483 0.483
time_pointFMT2 -0.665 0.483 0.483 0.534
typeCold_intervention:time_pointFMT1 0.470 -0.640 -0.341 -0.707 -0.377
typeWarm_control:time_pointFMT1 0.478 -0.347 -0.651 -0.719 -0.384 0.509
typeCold_intervention:time_pointFMT2 0.478 -0.651 -0.347 -0.384 -0.719 0.509 0.276
typeWarm_control:time_pointFMT2 0.478 -0.347 -0.651 -0.384 -0.719 0.272 0.518 0.518
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.13398738 -0.71904596 0.05318925 0.64850781 1.95572592
Number of Observations: 79
Number of Groups: 27 $emmeans
type emmean SE df lower.CL upper.CL
Cold_control 20.7 2.31 26 15.9 25.4
Cold_intervention 27.1 2.31 24 22.4 31.9
Warm_control 33.8 2.28 24 29.1 38.5
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_control - Cold_intervention -6.47 3.27 24 -1.982 0.1385
Cold_control - Warm_control -13.15 3.24 24 -4.057 0.0013
Cold_intervention - Warm_control -6.68 3.24 24 -2.061 0.1195
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 27.5 1.97 24 23.4 31.6
FMT1 22.5 1.97 24 18.4 26.5
FMT2 31.6 1.93 24 27.6 35.6
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 5.04 2.52 46 2.000 0.1235
Acclimation - FMT2 -4.14 2.49 46 -1.665 0.2296
FMT1 - FMT2 -9.18 2.49 46 -3.689 0.0017
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates Phylogenetic
Model_phylo <- nlme::lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
268.2771 293.0106 -123.1386
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.490394 1.134442
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 5.215724 0.4359152 46 11.964998 0.0000
typeCold_intervention 0.287458 0.5997818 24 0.479272 0.6361
typeWarm_control 1.265818 0.5997818 24 2.110465 0.0454
time_pointFMT1 -0.810372 0.5534409 46 -1.464243 0.1499
time_pointFMT2 -0.468071 0.5534409 46 -0.845747 0.4021
typeCold_intervention:time_pointFMT1 -0.550956 0.7826836 46 -0.703932 0.4850
typeWarm_control:time_pointFMT1 -1.415793 0.7696024 46 -1.839642 0.0723
typeCold_intervention:time_pointFMT2 0.044721 0.7696024 46 0.058109 0.9539
typeWarm_control:time_pointFMT2 -0.592994 0.7696024 46 -0.770520 0.4449
Correlation:
(Intr) typCl_ typWr_ t_FMT1 t_FMT2 tC_:_FMT1 tW_:_FMT1 tC_:_FMT2
typeCold_intervention -0.727
typeWarm_control -0.727 0.528
time_pointFMT1 -0.677 0.492 0.492
time_pointFMT2 -0.677 0.492 0.492 0.533
typeCold_intervention:time_pointFMT1 0.479 -0.652 -0.348 -0.707 -0.377
typeWarm_control:time_pointFMT1 0.487 -0.354 -0.664 -0.719 -0.383 0.508
typeCold_intervention:time_pointFMT2 0.487 -0.664 -0.354 -0.383 -0.719 0.508 0.276
typeWarm_control:time_pointFMT2 0.487 -0.354 -0.664 -0.383 -0.719 0.271 0.517 0.517
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.82031941 -0.51508546 -0.06170957 0.56728854 2.26098194
Number of Observations: 79
Number of Groups: 27 $emmeans
type emmean SE df lower.CL upper.CL
Cold_control 4.79 0.277 26 4.22 5.36
Cold_intervention 4.91 0.277 24 4.34 5.48
Warm_control 5.39 0.273 24 4.82 5.95
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_control - Cold_intervention -0.119 0.392 24 -0.303 0.9507
Cold_control - Warm_control -0.596 0.389 24 -1.534 0.2933
Cold_intervention - Warm_control -0.478 0.389 24 -1.229 0.4483
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 5.73 0.243 24 5.23 6.23
FMT1 4.27 0.243 24 3.77 4.77
FMT2 5.08 0.238 24 4.59 5.57
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 1.466 0.316 46 4.639 0.0001
Acclimation - FMT2 0.651 0.312 46 2.083 0.1045
FMT1 - FMT2 -0.815 0.312 46 -2.609 0.0321
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates 13.4.1.1 Beta diversity
Number of samples used
samples_to_keep_post7 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2") %>%
select(Tube_code) %>%
pull()
subset_meta_post7 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2")
subset_meta_post7$time_point<-as.factor(subset_meta_post7$time_point)
subset_meta_post7$type<-as.factor(subset_meta_post7$type)
length(samples_to_keep_post7)[1] 79Richness
richness_post7 <- as.matrix(beta_q0n$S)
richness_post7 <- as.dist(richness_post7[rownames(richness_post7) %in% samples_to_keep_post7,
colnames(richness_post7) %in% samples_to_keep_post7])
betadisper(richness_post7, subset_meta_post7$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.06472 0.032359 4.9377 999 0.006 **
Residuals 76 0.49806 0.006553
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.6710000 0.007
Cold_intervention 0.6794313 0.002
Warm_control 0.0122783 0.0036083 adonis2(richness_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 2 | 3.133065 | 0.12029225 | 5.585520 | 0.001 |
| time_point | 2 | 1.700292 | 0.06528175 | 3.031222 | 0.001 |
| type:time_point | 4 | 1.579667 | 0.06065041 | 1.408088 | 0.001 |
| Residual | 70 | 19.632418 | 0.75377559 | NA | NA |
| Total | 78 | 26.045441 | 1.00000000 | NA | NA |
Neutral
neutral_post7 <- as.matrix(beta_q1n$S)
neutral_post7 <- as.dist(neutral_post7[rownames(neutral_post7) %in% samples_to_keep_post7,
colnames(neutral_post7) %in% samples_to_keep_post7])
betadisper(neutral_post7, subset_meta_post7$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.06775 0.033873 3.8305 999 0.035 *
Residuals 76 0.67206 0.008843
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.1560000 0.214
Cold_intervention 0.1609162 0.010
Warm_control 0.1983916 0.0070193 adonis2(neutral_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 2 | 3.708525 | 0.14974144 | 7.624504 | 0.001 |
| time_point | 2 | 2.216763 | 0.08950765 | 4.557532 | 0.001 |
| type:time_point | 4 | 1.817057 | 0.07336847 | 1.867880 | 0.001 |
| Residual | 70 | 17.023843 | 0.68738244 | NA | NA |
| Total | 78 | 24.766189 | 1.00000000 | NA | NA |
Phylogenetic
phylo_post7 <- as.matrix(beta_q1p$S)
phylo_post7 <- as.dist(phylo_post7[rownames(phylo_post7) %in% samples_to_keep_post7,
colnames(phylo_post7) %in% samples_to_keep_post7])
betadisper(phylo_post7, subset_meta_post7$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.05474 0.027368 2.1613 999 0.118
Residuals 76 0.96240 0.012663
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_control Cold_intervention Warm_control
Cold_control 0.138000 0.717
Cold_intervention 0.151643 0.088
Warm_control 0.713310 0.083385 adonis2(phylo_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| type | 2 | 0.2951181 | 0.07519436 | 3.605402 | 0.001 |
| time_point | 2 | 0.5543931 | 0.14125613 | 6.772916 | 0.001 |
| type:time_point | 4 | 0.2103208 | 0.05358852 | 1.284725 | 0.145 |
| Residual | 70 | 2.8649048 | 0.72996099 | NA | NA |
| Total | 78 | 3.9247368 | 1.00000000 | NA | NA |
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_richness_nmds_post7 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
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")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_neutral_nmds_post7 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
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")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_phylogenetic_nmds_post7 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
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")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_nmds_post7, beta_neutral_nmds_post7, beta_phylogenetic_nmds_post7, ncol=3, nrow=1, common.legend = TRUE, legend="right")
13.5 Differences between acclimation and FMT1 across the three experimental groups
13.5.1 Beta diversity
CI from acclimation to FMT1
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Cold_intervention") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Cold_intervention")
length(samples_to_keep)[1] 17richness <- 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$time_point) %>% 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 1 0.001277 0.0012765 0.1315 999 0.694
Residuals 15 0.145582 0.0097055
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.697
FMT1 0.72191 adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.9572484 | 0.163835 | 2.939042 | 0.00390625 |
| Residual | 15 | 4.8855118 | 0.836165 | NA | NA |
| Total | 16 | 5.8427602 | 1.000000 | NA | NA |
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$time_point) %>% 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 1 0.004448 0.004448 0.3031 999 0.601
Residuals 15 0.220103 0.014674
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.581
FMT1 0.59003 adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 1.313310 | 0.2360329 | 4.634354 | 0.00390625 |
| Residual | 15 | 4.250788 | 0.7639671 | NA | NA |
| Total | 16 | 5.564099 | 1.0000000 | NA | NA |
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$time_point) %>% 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 1 0.02237 0.022372 0.8219 999 0.367
Residuals 15 0.40828 0.027219
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.347
FMT1 0.37895 adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
tt()| Df | SumOfSqs | R2 | F | Pr(>F) |
|---|---|---|---|---|
| 1 | 0.2548027 | 0.1800029 | 3.292748 | 0.03125 |
| 15 | 1.1607447 | 0.8199971 | NA | NA |
| 16 | 1.4155473 | 1.0000000 | NA | NA |
CC from acclimation to FMT1
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Cold_control") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Cold_control")
length(samples_to_keep)[1] 17richness <- 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$time_point) %>% 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 1 0.002483 0.0024832 0.2766 999 0.621
Residuals 15 0.134656 0.0089771
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.609
FMT1 0.60662 adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.6736847 | 0.1218982 | 2.082302 | 0.00390625 |
| Residual | 15 | 4.8529311 | 0.8781018 | NA | NA |
| Total | 16 | 5.5266158 | 1.0000000 | NA | NA |
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$time_point) %>% 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 1 0.002588 0.0025881 0.1669 999 0.72
Residuals 15 0.232575 0.0155050
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.714
FMT1 0.68863 adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.8134372 | 0.1708689 | 3.091228 | 0.00390625 |
| Residual | 15 | 3.9471558 | 0.8291311 | NA | NA |
| Total | 16 | 4.7605930 | 1.0000000 | NA | NA |
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$time_point) %>% 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 1 0.00219 0.0021905 0.1944 999 0.654
Residuals 15 0.16898 0.0112652
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.675
FMT1 0.66553 adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
tt()| Df | SumOfSqs | R2 | F | Pr(>F) |
|---|---|---|---|---|
| 1 | 0.1182834 | 0.1678279 | 3.025117 | 0.015625 |
| 15 | 0.5865068 | 0.8321721 | NA | NA |
| 16 | 0.7047902 | 1.0000000 | NA | NA |
WC from acclimation to FMT1
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Warm_control") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Warm_control")
length(samples_to_keep)[1] 18richness <- 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$time_point) %>% 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 1 0.022301 0.0223012 14.28 999 0.002 **
Residuals 16 0.024988 0.0015617
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.004
FMT1 0.0016445 adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.5832849 | 0.1336578 | 2.468453 | 0.00390625 |
| Residual | 16 | 3.7807317 | 0.8663422 | NA | NA |
| Total | 17 | 4.3640166 | 1.0000000 | NA | NA |
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$time_point) %>% 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 1 0.014714 0.0147136 2.6348 999 0.104
Residuals 16 0.089350 0.0055844
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.114
FMT1 0.12408 adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.697518 | 0.1715318 | 3.31275 | 0.00390625 |
| Residual | 16 | 3.368889 | 0.8284682 | NA | NA |
| Total | 17 | 4.066407 | 1.0000000 | NA | NA |
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$time_point) %>% 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 1 0.001474 0.0014742 0.1708 999 0.805
Residuals 16 0.138103 0.0086314
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Acclimation FMT1
Acclimation 0.837
FMT1 0.68489 adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
tt()| Df | SumOfSqs | R2 | F | Pr(>F) |
|---|---|---|---|---|
| 1 | 0.0924906 | 0.1569444 | 2.978582 | 0.01953125 |
| 16 | 0.4968302 | 0.8430556 | NA | NA |
| 17 | 0.5893208 | 1.0000000 | NA | NA |
13.5.2 Functional differences
CI from acclimation to FMT1
significant_element<-element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Cold_intervention") %>%
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, 8)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$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_element$trait,
colnames(element_gift)
))) %>%
left_join(., sample_metadata[c(1, 6, 8 )], by = join_by(Tube_code == Tube_code)) %>%
filter(time_point %in% c("FMT1","Acclimation"))%>%
filter(type=="Cold_intervention")
difference_table <- element_gift_sig %>%
select(-Tube_code, -type) %>%
group_by(time_point) %>%
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=Acclimation-FMT1)%>%
mutate(group_color = ifelse(Difference <0,"FMT1", "Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c("#76b183",'#008080')) +
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")
CC from acclimation to FMT1
significant_element<-element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Cold_control") %>%
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, 8)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$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_element$trait,
colnames(element_gift)
))) %>%
left_join(., sample_metadata[c(1, 6, 8)], by = join_by(Tube_code == Tube_code)) %>%
filter(time_point %in% c("FMT1","Acclimation"))%>%
filter(type=="Cold_control")
difference_table <- element_gift_sig %>%
select(-Tube_code, -type) %>%
group_by(time_point) %>%
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=Acclimation-FMT1)%>%
mutate(group_color = ifelse(Difference <0,"FMT1", "Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c("#4477AA",'#008080')) +
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")
WC from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Warm_control") %>%
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, 8)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$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"))# A tibble: 2 × 6
trait p_value p_adjust Domain Function Element
<chr> <dbl> <dbl> <chr> <chr> <chr>
1 D0801 0.000288 0.0210 Degradation Xenobiotic degradation_Toluene Toluene
2 D0802 0.000288 0.0210 Degradation Xenobiotic degradation_Xylene Xylene 13.6 Differences between FMT1 and FMT2 across the three experimental groups
13.6.1 Functional differences
CI from FMT1 to FMT2
significant_element<-element_gift %>%
filter(time_point %in% c("FMT2","FMT1") & type == "Cold_intervention") %>%
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, 8)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$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_element$trait,
colnames(element_gift)
))) %>%
left_join(., sample_metadata[c(1, 6, 8)], by = join_by(Tube_code == Tube_code)) %>%
filter(time_point %in% c("FMT2","FMT1"))%>%
filter(type=="Cold_intervention")
difference_table <- element_gift_sig %>%
select(-Tube_code, -type) %>%
group_by(time_point) %>%
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=FMT1-FMT2)%>%
mutate(group_color = ifelse(Difference <0,"FMT2", "FMT1")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c('#008080',"#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")
CC from FMT1 to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","FMT1") & type == "Cold_control") %>%
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, 8)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$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"))# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>WC from FMT1 to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","FMT1") & type == "Warm_control") %>%
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, 8)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$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"))# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>