3.1.1 LOAD

dat <- read_excel("./data/Data_collection.xlsx", sheet="pH_and_CFU", na = "NA")

3.1.2 LOOK

# Take a glimpse
glimpse(dat)

## Rows: 154
## Columns: 7
## $ Sample  <chr> "pAK80-D0", "Wg2-D0", "SK11-D0", "33-D0", "pAK80-P1-D1", "pAK8~
## $ Strain  <chr> "pAK80", "Wg2", "SK11", "33", "pAK80", "pAK80", "pAK80", "Wg2"~
## $ Protein <chr> "Inoculation", "Inoculation", "Inoculation", "Inoculation", "P~
## $ Day     <dbl> 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,~
## $ pH      <dbl> 7.20, 7.20, 7.20, 7.20, 6.96, 6.96, 6.96, 6.84, 6.85, 6.84, 6.~
## $ CFU     <dbl> 6.20e+06, 6.54e+06, 6.00e+06, 6.43e+06, 1.70e+08, 2.42e+08, 1.~
## $ logN    <dbl> 6.792392, 6.815578, 6.778151, 6.808211, 8.230449, 8.383815, 8.~

# Make a explorative summary 
skimr::skim(dat)

Variable type: character

Variable type: numeric

3.1.3 SAVE

# Save cleaned data
save(dat, file = "R_objects/stat_test_data_growth.RData")

# clear the environment and release memory
rm(list = ls(all.names = TRUE)[ls(all.names = TRUE) != "params"])
invisible(gc())

3.2.1 CATEGORICAL VARIABLES

Chi squared test of independence.

# load data 
load("R_objects/stat_test_data_growth.RData")

# This converts all character columns to factor (if only specific columns should be converted replace "where(is_character)" with "c(<col1>,<col2>,...)"
dat <- dat %>% mutate(across(where(is_character), as_factor))


# Create vector of categorical variables with more than two categories and less than half of the number of samples.
cat_var <- dat %>% 
  select_if(is.factor) %>% 
  select(where(~n_distinct(.) > 1)) %>% 
  select(where(~n_distinct(.) < (nrow(dat)/2))) %>%
  colnames()

length(cat_var) # If less than 2 variables the following code should not be run

## [1] 2

# Test the variables pairwise 
cat_var %>% 
  combn(2) %>% 
  t() %>% 
  as_tibble() %>% 
  rowwise %>% 
  mutate(chisq_test = list(
    table(dat[[V1]], dat[[V2]]) %>% chisq.test()
    ),
    chisq_pval = chisq_test$p.value
    )

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## i Using compatibility `.name_repair`.

## Warning: There was 1 warning in `mutate()`.
## i In argument: `chisq_test = list(table(dat[[V1]], dat[[V2]]) %>%
##   chisq.test())`.
## i In row 1.
## Caused by warning in `chisq.test()`:
## ! Chi-squared approximation may be incorrect

## # A tibble: 1 x 4
## # Rowwise: 
##   V1     V2      chisq_test chisq_pval
##   <chr>  <chr>   <list>          <dbl>
## 1 Strain Protein <htest>          1.00

3.2.2 QUANTITATIVE VARIABLES

Spearman’s test without first testing the normality of the data

# create vector with the relevant variables 
CON.VARS <- dat %>% select(where(is.numeric)) %>% colnames()

# Plot all variables against each other
pairs(dat[,CON.VARS], pch = 19,  cex = 0.5,
      lower.panel=NULL)

# Run Spearman test
corrmat <- cor(dat[,CON.VARS], method = "spearman", use = "pairwise.complete.obs")

# Create heatmap
corrmat_rounded <- round(corrmat, 2)

melted_corrmat_rounded <- tibble(Var1 = rep(row.names(corrmat_rounded), length(row.names(corrmat_rounded))),
                                 Var2 = rep(row.names(corrmat_rounded), each = length(row.names(corrmat_rounded))),
                                 dist = as.numeric(matrix(corrmat_rounded)))

ggplot(melted_corrmat_rounded, aes(x = Var1, y = Var2, fill = dist)) + 
  labs(title = "Correlation between continous variables") + 
  geom_tile(color = "white") + 
  scale_fill_gradient2(low = "blue", mid = "white", 
                       high = "red", midpoint = 0, limit = c(-1,1), 
                       space = "Lab", name = "Correlation coefficient") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, size = 12, 
                                   hjust = 1, vjust = 1), 
        axis.text.y = element_text(size = 12, hjust = 1, vjust = 0.5),
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank()) + 
  geom_text(aes(x = Var1, y = Var2, label = dist), 
            color = "black", size = 4) + 
  coord_fixed()

# clear the environment and release memory
rm(list = ls(all.names = TRUE)[ls(all.names = TRUE) != "params"])
invisible(gc())

3.3.1 PREPARE VARIABLES

# Set names of variables
PREDICTOR <- c("Day","Strain", "Protein")
OUTCOME <- "pH"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 52 x 7
##    Strain Protein     Day   variable     n  mean    sd
##    <chr>  <chr>       <fct> <fct>    <dbl> <dbl> <dbl>
##  1 33     Casein      1     pH           3  4.99 0.015
##  2 33     Hydrolysate 1     pH           2  6.97 0    
##  3 33     Potato      1     pH           3  6.74 0.032
##  4 33     Water       1     pH           2  7.02 0.007
##  5 No     Casein      1     pH           3  6.97 0.01 
##  6 No     Hydrolysate 1     pH           2  6.98 0.028
##  7 No     Potato      1     pH           3  7.00 0.023
##  8 No     Water       1     pH           2  7.04 0    
##  9 SK11   Casein      1     pH           3  6.96 0    
## 10 SK11   Hydrolysate 1     pH           2  7    0    
## # i 42 more rows

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

3.3.2 PRELIMINARY TEST

Identify outliers

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 x 9
##   Strain Protein Day   Sample         pH       CFU  logN is.outlier is.extreme
##   <chr>  <chr>   <fct> <chr>       <dbl>     <dbl> <dbl> <lgl>      <lgl>     
## 1 Wg2    Potato  6     Wg2-P3-D6    5.02   6040000  6.78 TRUE       FALSE     
## 2 pAK80  Casein  6     pAK80-C4-D6  6.66 580000000  8.76 TRUE       FALSE

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.The Shapiro-Wilk test is performed on all data groups when possible with n > 2 and only with unique values.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)

# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 x 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.281 2.47e-23

# Compute Shapiro-Wilk test of normality for each data group
sets <- dat.clean %>%
  group_by(Day, Protein, Strain) %>% 
  summarise(n_unique = n_distinct(pH),
            n = n()) %>%
  filter(n > 2, n_unique > 1) %>% 
  transmute(day_protein_strain = paste(Day, Protein, Strain))

## `summarise()` has grouped output by 'Day', 'Protein'. You can override using
## the `.groups` argument.

dat.clean %>% 
  mutate(day_protein_strain = paste(Day, Protein, Strain)) %>%
  filter(day_protein_strain %in% sets$day_protein_strain) %>%
  group_by(Day, Protein, Strain) %>% 
  shapiro_test(vars = "pH") %>% 
  adjust_pvalue(method = "fdr") %>%
  filter(p.adj < 0.05)

## # A tibble: 7 x 7
##   Strain Protein Day   variable statistic       p   p.adj
##   <chr>  <chr>   <fct> <chr>        <dbl>   <dbl>   <dbl>
## 1 No     Potato  1     pH           0.75  0       0      
## 2 Wg2    Potato  1     pH           0.75  0       0      
## 3 SK11   Casein  2     pH           0.75  0       0      
## 4 Wg2    Casein  2     pH           0.75  0       0      
## 5 Wg2    Potato  2     pH           0.75  0       0      
## 6 pAK80  Potato  2     pH           0.75  0       0      
## 7 Wg2    Potato  6     pH           0.630 0.00124 0.00443

Test homogneity of variance assumption

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 x 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1    51    89      46.8 6.76e-47

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

3.3.3 PERFORM TEST

The data are not normal distributed. The Shapiro-Wilk normality test found that at least one group was not normal distributed. As a result, for statistical data analysis, the Kruskal-Wallis one-way analysis of variance test is used, followed by a Wilcoxon test for pairwise group comparisons.

# for each protein source
res.protein <- dat.clean %>% 
  group_by(Strain, Day) %>%
  kruskal_test(pH ~ Protein) %>%
  adjust_pvalue(method = "fdr")
res.protein %>% filter(p < 0.05)

## # A tibble: 10 x 9
##    Strain Day   .y.       n statistic    df      p method          p.adj
##    <chr>  <fct> <chr> <int>     <dbl> <int>  <dbl> <chr>           <dbl>
##  1 33     1     pH       10      8.51     3 0.0366 Kruskal-Wallis 0.0683
##  2 33     2     pH        8      6.25     2 0.0439 Kruskal-Wallis 0.0683
##  3 33     6     pH       12      8.81     3 0.032  Kruskal-Wallis 0.0683
##  4 No     6     pH       12      7.85     3 0.0493 Kruskal-Wallis 0.0690
##  5 SK11   1     pH       10      8.18     3 0.0424 Kruskal-Wallis 0.0683
##  6 SK11   6     pH       12     10.2      3 0.0167 Kruskal-Wallis 0.0683
##  7 Wg2    1     pH       10      8.56     3 0.0358 Kruskal-Wallis 0.0683
##  8 Wg2    2     pH        8      6.40     2 0.0407 Kruskal-Wallis 0.0683
##  9 Wg2    6     pH       12      9.71     3 0.0212 Kruskal-Wallis 0.0683
## 10 pAK80  6     pH       12     10.2      3 0.0173 Kruskal-Wallis 0.0683

# for each strain
res.strain <- dat.clean %>% 
  group_by(Day, Protein) %>%
  kruskal_test(pH ~ Strain) %>%
  adjust_pvalue(method = "fdr")
res.strain %>% filter(p < 0.05)

## # A tibble: 6 x 9
##   Protein Day   .y.       n statistic    df       p method           p.adj
##   <chr>   <fct> <chr> <int>     <dbl> <int>   <dbl> <chr>            <dbl>
## 1 Casein  1     pH       15     11.9      4 0.0182  Kruskal-Wallis 0.0400 
## 2 Potato  1     pH       15     13.4      4 0.00943 Kruskal-Wallis 0.0346 
## 3 Casein  2     pH       12     10.5      3 0.0151  Kruskal-Wallis 0.0400 
## 4 Potato  2     pH       11      9.54     3 0.0229  Kruskal-Wallis 0.0420 
## 5 Casein  6     pH       20     17.6      4 0.00151 Kruskal-Wallis 0.00830
## 6 Potato  6     pH       20     18.4      4 0.00105 Kruskal-Wallis 0.00830

pwc.strain <- dat.clean %>%
  filter(Protein %in% c("Casein", "Potato"), Day==6) %>%
  group_by(Protein) %>%
  pairwise_wilcox_test(pH ~ Strain) %>%
  adjust_pvalue(method = "fdr")
pwc.strain %>% filter(p.adj < 0.05)

## # A tibble: 19 x 10
##    Protein .y.   group1 group2    n1    n2 statistic     p  p.adj p.adj.signif
##    <chr>   <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>  <dbl> <chr>       
##  1 Casein  pH    33     No         4     4       0   0.028 0.0322 ns          
##  2 Casein  pH    33     pAK80      4     4       0   0.029 0.0322 ns          
##  3 Casein  pH    33     SK11       4     4       0   0.028 0.0322 ns          
##  4 Casein  pH    33     Wg2        4     4       0.5 0.041 0.0432 ns          
##  5 Casein  pH    No     pAK80      4     4      16   0.028 0.0322 ns          
##  6 Casein  pH    No     SK11       4     4      16   0.028 0.0322 ns          
##  7 Casein  pH    No     Wg2        4     4      16   0.028 0.0322 ns          
##  8 Casein  pH    pAK80  Wg2        4     4      16   0.029 0.0322 ns          
##  9 Casein  pH    SK11   Wg2        4     4      16   0.029 0.0322 ns          
## 10 Potato  pH    33     No         4     4       0   0.029 0.0322 ns          
## 11 Potato  pH    33     pAK80      4     4       0   0.029 0.0322 ns          
## 12 Potato  pH    33     SK11       4     4       0   0.029 0.0322 ns          
## 13 Potato  pH    33     Wg2        4     4       0   0.026 0.0322 ns          
## 14 Potato  pH    No     pAK80      4     4      16   0.029 0.0322 ns          
## 15 Potato  pH    No     SK11       4     4      16   0.028 0.0322 ns          
## 16 Potato  pH    No     Wg2        4     4      16   0.026 0.0322 ns          
## 17 Potato  pH    pAK80  SK11       4     4      16   0.029 0.0322 ns          
## 18 Potato  pH    pAK80  Wg2        4     4      16   0.026 0.0322 ns          
## 19 Potato  pH    SK11   Wg2        4     4      16   0.026 0.0322 ns

3.3.4 VISUALIZE

#Name order in plot
dat.clean$Day <- as.factor(dat.clean$Day)
dat.clean$Strain <- factor(dat.clean$Strain, 
                           levels = c("No", "pAK80", "SK11", "Wg2", "33"), 
                           labels =c ("Non-inoculated", 
                                      "MS22418 (PrtP-)", 
                                      "MS22421 (PrtP/SK11)", 
                                      "MS22425 (PrtP/Wg2)", 
                                      "MS22427 (PrtP/MS22333)"))
                           
dat.clean$Protein <- factor(dat.clean$Protein, levels = c("Casein", "Potato", "Hydrolysate", "Water"))

#Boxplot
bxp.pH <- dat.clean %>%
  filter(., Protein  %in% c("Casein", "Potato", "Hydrolysate", "Water")) %>%
  filter (., Day %in% c(1, 2, 6)) %>%
  ggboxplot(.,
            x = "Day",
            y = "pH",
            color = "Strain", 
            facet.by = "Protein")+
  labs(x="Time [days]", y="pH")+
  scale_y_continuous(limits = c(4.8, 7.2))+
  facet_wrap(~Protein, 
             nrow=2, 
             labeller = as_labeller(c("Casein"="CCDM", 
                                      "Potato"="PCDM", 
                                      "Hydrolysate"="HCDM", 
                                      "Water"="CDM")))+
  scale_color_manual(values = c("Non-inoculated"="#E89242FF", 
                                  "MS22418 (PrtP-)"="#82491EFF", 
                                  "MS22421 (PrtP/SK11)"="#24325FFF", 
                                  "MS22425 (PrtP/Wg2)"="#B7E4F9FF", 
                                  "MS22427 (PrtP/MS22333)"="#FB6467FF")
                     )+
  geom_hline(yintercept=7.0, linetype="dashed", color = "#7d7d7dff")+
  theme_light()+
  theme(legend.position="bottom")
 
save(bxp.pH, file = "R_objects/Fig1A.RData")

3.4.1 DATA

# load data 
load("R_objects/stat_test_data_growth.RData")

# Add replicate identifier
dat$replicate <- str_split(dat$Sample, "-", simplify = T)[,2]

#Calculate pH change
dat.wide <- dat %>%
  pivot_wider(values_from = pH, names_from = Day, names_prefix = "Day",id_cols = c(replicate, Strain, Protein), values_fill = NA)

dat.wide <- dat.wide %>%
  mutate(pH_1 = Day1 - 7.2,
         pH_2 = Day2 - Day1,
         pH_6 = Day6 - Day2,
         norm_1 = 100*Day1/7.2,
         norm_2 = 100*Day2/7.2,
         norm_6 = 100*Day6/7.2
         )

ggplot(dat.wide, aes(x = Protein, y = pH_1, color = Strain)) +
  geom_boxplot(outlier.colour = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1,
                                             dodge.width = 0.7))

## Warning: Removed 14 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 14 rows containing missing values (`geom_point()`).

ggplot(dat.wide, aes(x = Protein, y = pH_2, color = Strain)) +
  geom_boxplot(outlier.colour = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1,
                                             dodge.width = 0.7))

## Warning: Removed 33 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 33 rows containing missing values (`geom_point()`).

ggplot(dat.wide, aes(x = Protein, y = pH_6, color = Strain)) +
  geom_boxplot(outlier.colour = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1,
                                             dodge.width = 0.7))

## Warning: Removed 33 rows containing non-finite values (`stat_boxplot()`).
## Removed 33 rows containing missing values (`geom_point()`).

# Create subset
dat.clean.wide <- dat.wide %>% pivot_longer(names_to = "Day",values_to = "pH_change",cols = c(pH_1, pH_2, pH_6))

ggplot(dat.clean.wide, aes(x = Day, y = pH_change, color = Protein, group = replicate)) +
  geom_point() + geom_line() +
  facet_wrap("Strain")

## Warning: Removed 80 rows containing missing values (`geom_point()`).

## Warning: Removed 46 rows containing missing values (`geom_line()`).

Save data

# Save cleaned data
save(dat.clean.wide, file = "R_objects/stat_test_data_pH_change.RData")

# clear the environment and release memory
rm(list = ls(all.names = TRUE)[ls(all.names = TRUE) != "params"])
invisible(gc())

3.4.2 ALL - Strain & Protein - pH_change

# load data 
load("R_objects/stat_test_data_pH_change.RData")

# Create subset
dat.clean <- 
  dat.clean.wide %>% 
  filter(!is.na(pH_change))

# Set names of variables
PREDICTOR <- c("Protein","Strain", "Day")
OUTCOME <- "pH_change"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 44 x 7
##    Strain Protein Day   variable      n   mean    sd
##    <chr>  <chr>   <chr> <fct>     <dbl>  <dbl> <dbl>
##  1 33     Casein  pH_1  pH_change     3 -2.21  0.015
##  2 33     Casein  pH_2  pH_change     3  0.207 0.021
##  3 33     Casein  pH_6  pH_change     3 -0.2   0    
##  4 No     Casein  pH_1  pH_change     3 -0.23  0.01 
##  5 SK11   Casein  pH_1  pH_change     3 -0.24  0    
##  6 SK11   Casein  pH_2  pH_change     3  0.117 0.006
##  7 SK11   Casein  pH_6  pH_change     3 -0.283 0.015
##  8 Wg2    Casein  pH_1  pH_change     3 -2.17  0.026
##  9 Wg2    Casein  pH_2  pH_change     3  0.083 0.023
## 10 Wg2    Casein  pH_6  pH_change     3 -0.073 0.012
## # i 34 more rows

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

##  [1] Strain     Protein    Day        replicate  Day0       Day1      
##  [7] Day2       Day6       norm_1     norm_2     norm_6     pH_change 
## [13] is.outlier is.extreme
## <0 rækker> (eller 0-længde row.names)

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)

# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 x 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.363 6.80e-20

# Compute Shapiro-Wilk test of normality for each data group
sets <- dat.clean %>%
  group_by(Day, Protein, Strain) %>% 
  summarise(n_unique = n_distinct(pH_change),
            n = n()) %>%
  filter(n > 2, n_unique > 1) %>% 
  transmute(day_protein_strain = paste(Day, Protein, Strain))

## `summarise()` has grouped output by 'Day', 'Protein'. You can override using
## the `.groups` argument.

dat.clean %>% 
  mutate(day_protein_strain = paste(Day, Protein, Strain)) %>%
  filter(day_protein_strain %in% sets$day_protein_strain) %>%
  group_by(Day, Protein, Strain) %>% 
  shapiro_test(vars = "pH_change") %>% 
  adjust_pvalue(method = "fdr") %>%
  filter(p.adj < 0.05)

## # A tibble: 8 x 7
##   Strain Protein Day   variable  statistic     p p.adj
##   <chr>  <chr>   <chr> <chr>         <dbl> <dbl> <dbl>
## 1 No     Potato  pH_1  pH_change      0.75     0     0
## 2 Wg2    Potato  pH_1  pH_change      0.75     0     0
## 3 SK11   Casein  pH_2  pH_change      0.75     0     0
## 4 Wg2    Casein  pH_2  pH_change      0.75     0     0
## 5 Wg2    Potato  pH_2  pH_change      0.75     0     0
## 6 pAK80  Potato  pH_2  pH_change      0.75     0     0
## 7 Wg2    Casein  pH_6  pH_change      0.75     0     0
## 8 pAK80  Casein  pH_6  pH_change      0.75     0     0

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 x 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1    43    68      75.4 1.48e-43

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

# for each protein source
res.protein <- dat.clean %>% 
  group_by(Strain, Day) %>%
  kruskal_test(pH_change ~ Protein) %>%
  adjust_pvalue(method = "fdr")
res.protein %>% filter(p < 0.05)

## # A tibble: 8 x 9
##   Strain Day   .y.           n statistic    df      p method          p.adj
##   <chr>  <chr> <chr>     <int>     <dbl> <int>  <dbl> <chr>           <dbl>
## 1 33     pH_1  pH_change    10      8.51     3 0.0366 Kruskal-Wallis 0.0713
## 2 33     pH_2  pH_change     8      6.25     2 0.0439 Kruskal-Wallis 0.0713
## 3 33     pH_6  pH_change     8      6.56     2 0.0376 Kruskal-Wallis 0.0713
## 4 SK11   pH_1  pH_change    10      8.18     3 0.0424 Kruskal-Wallis 0.0713
## 5 Wg2    pH_1  pH_change    10      8.56     3 0.0358 Kruskal-Wallis 0.0713
## 6 Wg2    pH_2  pH_change     8      6.40     2 0.0407 Kruskal-Wallis 0.0713
## 7 Wg2    pH_6  pH_change     8      6.33     2 0.0423 Kruskal-Wallis 0.0713
## 8 pAK80  pH_6  pH_change     8      6.33     2 0.0423 Kruskal-Wallis 0.0713

# for each strain
res.strain <- dat.clean %>% 
  group_by(Day, Protein) %>%
  kruskal_test(pH_change ~ Strain) %>%
  adjust_pvalue(method = "fdr")
res.strain %>% filter(p < 0.05)

## # A tibble: 6 x 9
##   Protein Day   .y.           n statistic    df       p method          p.adj
##   <chr>   <chr> <chr>     <int>     <dbl> <int>   <dbl> <chr>           <dbl>
## 1 Casein  pH_1  pH_change    15     11.9      4 0.0182  Kruskal-Wallis 0.0468
## 2 Potato  pH_1  pH_change    15     13.4      4 0.00943 Kruskal-Wallis 0.0468
## 3 Casein  pH_2  pH_change    12     10.3      3 0.0164  Kruskal-Wallis 0.0468
## 4 Potato  pH_2  pH_change    11      9.50     3 0.0234  Kruskal-Wallis 0.0468
## 5 Casein  pH_6  pH_change    12      9.84     3 0.02    Kruskal-Wallis 0.0468
## 6 Potato  pH_6  pH_change    11      8.77     3 0.0326  Kruskal-Wallis 0.0543

pwc.strain <- dat.clean %>%
  filter(Protein %in% c("Casein", "Potato") & Day %in% c("pH_1")) %>%
  group_by(Day, Protein) %>%
  pairwise_wilcox_test(pH_change ~ Strain) %>%
  adjust_pvalue(method = "fdr")
pwc.strain %>% filter (p < 0.5)

## # A tibble: 19 x 11
##    Protein Day   .y.       group1 group2    n1    n2 statistic     p p.adj
##    <chr>   <chr> <chr>     <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
##  1 Casein  pH_1  pH_change 33     No         3     3       0   0.1   0.133
##  2 Casein  pH_1  pH_change 33     pAK80      3     3       0   0.1   0.133
##  3 Casein  pH_1  pH_change 33     SK11       3     3       0   0.064 0.133
##  4 Casein  pH_1  pH_change 33     Wg2        3     3       1   0.2   0.211
##  5 Casein  pH_1  pH_change No     SK11       3     3       7.5 0.197 0.211
##  6 Casein  pH_1  pH_change No     Wg2        3     3       9   0.1   0.133
##  7 Casein  pH_1  pH_change pAK80  SK11       3     3       7.5 0.197 0.211
##  8 Casein  pH_1  pH_change pAK80  Wg2        3     3       9   0.1   0.133
##  9 Casein  pH_1  pH_change SK11   Wg2        3     3       9   0.064 0.133
## 10 Potato  pH_1  pH_change 33     No         3     3       0   0.076 0.133
## 11 Potato  pH_1  pH_change 33     pAK80      3     3       0   0.064 0.133
## 12 Potato  pH_1  pH_change 33     SK11       3     3       0   0.1   0.133
## 13 Potato  pH_1  pH_change 33     Wg2        3     3       0   0.076 0.133
## 14 Potato  pH_1  pH_change No     pAK80      3     3       9   0.059 0.133
## 15 Potato  pH_1  pH_change No     SK11       3     3       9   0.076 0.133
## 16 Potato  pH_1  pH_change No     Wg2        3     3       9   0.072 0.133
## 17 Potato  pH_1  pH_change pAK80  SK11       3     3       7.5 0.197 0.211
## 18 Potato  pH_1  pH_change pAK80  Wg2        3     3       9   0.059 0.133
## 19 Potato  pH_1  pH_change SK11   Wg2        3     3       9   0.076 0.133
## # i 1 more variable: p.adj.signif <chr>

3.4.3 pH_1 - Strain & Protein - pH_change

This subset analysis focus on the first pH change, but the analysis did not contribute to the main conclusion of the paper. Therefore, this analysis is not presented or discussed in the paper.

# load data 
load("R_objects/stat_test_data_pH_change.RData")

# Create subset
dat.clean <- 
  dat.clean.wide %>% 
  filter(!is.na(pH_change) & Day == "pH_1")

# Set names of variables
PREDICTOR <- c("Protein","Strain")
OUTCOME <- "pH_change"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 20 x 6
##    Strain Protein     variable      n   mean    sd
##    <chr>  <chr>       <fct>     <dbl>  <dbl> <dbl>
##  1 33     Casein      pH_change     3 -2.21  0.015
##  2 No     Casein      pH_change     3 -0.23  0.01 
##  3 SK11   Casein      pH_change     3 -0.24  0    
##  4 Wg2    Casein      pH_change     3 -2.17  0.026
##  5 pAK80  Casein      pH_change     3 -0.207 0.049
##  6 33     Hydrolysate pH_change     2 -0.23  0    
##  7 No     Hydrolysate pH_change     2 -0.22  0.028
##  8 SK11   Hydrolysate pH_change     2 -0.2   0    
##  9 Wg2    Hydrolysate pH_change     2 -0.19  0.014
## 10 pAK80  Hydrolysate pH_change     2 -0.16  0.014
## 11 33     Potato      pH_change     3 -0.457 0.032
## 12 No     Potato      pH_change     3 -0.197 0.023
## 13 SK11   Potato      pH_change     3 -0.257 0.021
## 14 Wg2    Potato      pH_change     3 -0.357 0.006
## 15 pAK80  Potato      pH_change     3 -0.24  0    
## 16 33     Water       pH_change     2 -0.175 0.007
## 17 No     Water       pH_change     2 -0.16  0    
## 18 SK11   Water       pH_change     2 -0.155 0.021
## 19 Wg2    Water       pH_change     2 -0.12  0    
## 20 pAK80  Water       pH_change     2 -0.11  0.014

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

##  [1] Strain     Protein    replicate  Day0       Day1       Day2      
##  [7] Day6       norm_1     norm_2     norm_6     Day        pH_change 
## [13] is.outlier is.extreme
## <0 rækker> (eller 0-længde row.names)

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)

# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 x 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.942  0.0156

# Compute Shapiro-Wilk test of normality for each data group
sets <- dat.clean %>%
  group_by(Day, Protein, Strain) %>% 
  summarise(n_unique = n_distinct(pH_change),
            n = n()) %>%
  filter(n > 2, n_unique > 1) %>% 
  transmute(day_protein_strain = paste(Day, Protein, Strain))

## `summarise()` has grouped output by 'Day', 'Protein'. You can override using
## the `.groups` argument.

dat.clean %>% 
  mutate(day_protein_strain = paste(Day, Protein, Strain)) %>%
  filter(day_protein_strain %in% sets$day_protein_strain) %>%
  group_by(Day, Protein, Strain) %>% 
  shapiro_test(vars = "pH_change") %>% 
  adjust_pvalue(method = "fdr") %>%
  filter(p.adj < 0.05)

## # A tibble: 2 x 7
##   Strain Protein Day   variable  statistic     p p.adj
##   <chr>  <chr>   <chr> <chr>         <dbl> <dbl> <dbl>
## 1 No     Potato  pH_1  pH_change      0.75     0     0
## 2 Wg2    Potato  pH_1  pH_change      0.75     0     0

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 x 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1    19    30     0.708 0.782

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

# for each protein source
res.protein <- dat.clean %>% 
  group_by(Strain, Day) %>%
  kruskal_test(pH_change ~ Protein) %>%
  adjust_pvalue(method = "fdr")
res.protein %>% filter(p < 0.05)

## # A tibble: 3 x 9
##   Strain Day   .y.           n statistic    df      p method          p.adj
##   <chr>  <chr> <chr>     <int>     <dbl> <int>  <dbl> <chr>           <dbl>
## 1 33     pH_1  pH_change    10      8.51     3 0.0366 Kruskal-Wallis 0.0707
## 2 SK11   pH_1  pH_change    10      8.18     3 0.0424 Kruskal-Wallis 0.0707
## 3 Wg2    pH_1  pH_change    10      8.56     3 0.0358 Kruskal-Wallis 0.0707

pwc.protein <- dat.clean %>%
  filter(Protein %in% c("Casein", "Potato")) %>%
  group_by(Strain, Day) %>%
  pairwise_wilcox_test(pH_change ~ Protein) %>%
  adjust_pvalue(method = "fdr")
pwc.protein %>% filter(p < 0.05)

## # A tibble: 0 x 11
## # i 11 variables: Strain <chr>, Day <chr>, .y. <chr>, group1 <chr>,
## #   group2 <chr>, n1 <int>, n2 <int>, statistic <dbl>, p <dbl>, p.adj <dbl>,
## #   p.adj.signif <chr>

# for each strain
res.strain <- dat.clean %>% 
  group_by(Day, Protein) %>%
  kruskal_test(pH_change ~ Strain) %>%
  adjust_pvalue(method = "fdr")
res.strain %>% filter(p < 0.05)

## # A tibble: 2 x 9
##   Protein Day   .y.           n statistic    df       p method          p.adj
##   <chr>   <chr> <chr>     <int>     <dbl> <int>   <dbl> <chr>           <dbl>
## 1 Casein  pH_1  pH_change    15      11.9     4 0.0182  Kruskal-Wallis 0.0364
## 2 Potato  pH_1  pH_change    15      13.4     4 0.00943 Kruskal-Wallis 0.0364

pwc.strain <- dat.clean %>%
  filter(Protein %in% c("Casein", "Potato")) %>%
  group_by(Protein) %>%
  pairwise_wilcox_test(pH_change ~ Strain) %>%
  adjust_pvalue(method = "fdr")
pwc.strain %>% filter(p < 0.05)

## # A tibble: 0 x 10
## # i 10 variables: Protein <chr>, .y. <chr>, group1 <chr>, group2 <chr>,
## #   n1 <int>, n2 <int>, statistic <dbl>, p <dbl>, p.adj <dbl>,
## #   p.adj.signif <chr>

3.5.1 DATA

# load data 
load("R_objects/stat_test_data_growth.RData")

# Create subset
dat.clean <- dat %>% 
  filter(!is.na(CFU) & Day != 0) %>%
  subset(Protein %in% c("Potato", "Casein"))%>%
  mutate(Day = as.factor(Day))

3.5.2 PREPARE VARIABLES

# Set names of variables

PREDICTOR <- c("Protein","Strain","Day")
OUTCOME <- "logN"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 16 x 7
##    Strain Protein Day   variable     n  mean    sd
##    <chr>  <chr>   <fct> <fct>    <dbl> <dbl> <dbl>
##  1 33     Casein  1     logN         3  9.21 0.056
##  2 33     Casein  6     logN         2  4.11 0.038
##  3 SK11   Casein  1     logN         3  7.79 0.077
##  4 SK11   Casein  6     logN         3  6.17 0.785
##  5 Wg2    Casein  1     logN         3  9.13 0.019
##  6 Wg2    Casein  6     logN         2  4.46 0.025
##  7 pAK80  Casein  1     logN         3  7.38 0.075
##  8 pAK80  Casein  6     logN         4  7.30 1.02 
##  9 33     Potato  1     logN         3  8.98 0.091
## 10 33     Potato  6     logN         4  6.82 0.034
## 11 SK11   Potato  1     logN         3  8.20 0.509
## 12 SK11   Potato  6     logN         4  7.15 0.14 
## 13 Wg2    Potato  1     logN         3  8.79 0.051
## 14 Wg2    Potato  6     logN         4  6.75 0.101
## 15 pAK80  Potato  1     logN         3  8.26 0.107
## 16 pAK80  Potato  6     logN         4  6.86 0.176

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

3.5.3 PRELIMINARY TEST

Identify outliers

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 x 9
##   Strain Protein Day   Sample        pH     CFU  logN is.outlier is.extreme
##   <chr>  <chr>   <fct> <chr>      <dbl>   <dbl> <dbl> <lgl>      <lgl>     
## 1 SK11   Potato  6     SK11-P1-D6  5.99 8800000  6.94 TRUE       FALSE     
## 2 Wg2    Potato  6     Wg2-P2-D6   5    4020000  6.60 TRUE       FALSE

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.The Shapiro-Wilk test is performed on all data groups when possible with n > 2 and only with unique values.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)

# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 x 3
##   variable         statistic      p.value
##   <chr>                <dbl>        <dbl>
## 1 residuals(model)     0.724 0.0000000182

# Compute Shapiro-Wilk test of normality for each data group
sets <- dat.clean %>%
  group_by(Day, Protein, Strain) %>% 
  summarise(n_unique = n_distinct(logN),
            n = n()) %>%
  filter(n > 2, n_unique > 1) %>% 
  transmute(day_protein_strain = paste(Day, Protein, Strain))

## `summarise()` has grouped output by 'Day', 'Protein'. You can override using
## the `.groups` argument.

dat.clean %>% 
  mutate(day_protein_strain = paste(Day, Protein, Strain)) %>%
  filter(day_protein_strain %in% sets$day_protein_strain) %>%
  group_by(Day, Protein, Strain) %>% 
  shapiro_test(vars = "logN") %>% 
  adjust_pvalue(method = "fdr") %>%
  filter(p.adj < 0.05)

## # A tibble: 1 x 7
##   Strain Protein Day   variable statistic     p p.adj
##   <chr>  <chr>   <fct> <chr>        <dbl> <dbl> <dbl>
## 1 SK11   Casein  1     logN          0.75     0     0

Test homogneity of variance assumption

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 x 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1    15    35      1.47 0.171

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

3.5.4 PERFORM TEST

The data are not normal distributed. The Shapiro-Wilk normality test found that at least one group was not normal distributed. As a result, for statistical data analysis, the Kruskal-Wallis one-way analysis of variance test is used, followed by a Wilcoxon test for pairwise group comparisons.

# for each protein source
res.protein <- dat.clean %>% 
  group_by(Day, Strain) %>%
  kruskal_test(logN ~ Protein) %>%
  adjust_pvalue(method = "fdr")
res.protein %>% filter(p < 0.05)

## # A tibble: 3 x 9
##   Strain Day   .y.       n statistic    df      p method         p.adj
##   <chr>  <fct> <chr> <int>     <dbl> <int>  <dbl> <chr>          <dbl>
## 1 33     1     logN      6      3.86     1 0.0495 Kruskal-Wallis 0.103
## 2 Wg2    1     logN      6      3.86     1 0.0495 Kruskal-Wallis 0.103
## 3 pAK80  1     logN      6      3.86     1 0.0495 Kruskal-Wallis 0.103

# for each strain
res.strain <- dat.clean %>% 
  group_by(Day, Protein) %>%
  kruskal_test(logN ~ Strain) %>%
  adjust_pvalue(method = "fdr")
res.strain %>% filter(p < 0.05)

## # A tibble: 4 x 9
##   Protein Day   .y.       n statistic    df      p method          p.adj
##   <chr>   <fct> <chr> <int>     <dbl> <int>  <dbl> <chr>           <dbl>
## 1 Casein  1     logN     12     10.2      3 0.0166 Kruskal-Wallis 0.0407
## 2 Potato  1     logN     12      9.46     3 0.0237 Kruskal-Wallis 0.0407
## 3 Casein  6     logN     11      8.21     3 0.0418 Kruskal-Wallis 0.0418
## 4 Potato  6     logN     16      8.91     3 0.0305 Kruskal-Wallis 0.0407

pwc <- dat.clean %>%
  group_by(Day, Protein) %>%
  pairwise_wilcox_test(logN ~ Strain) %>%
  adjust_pvalue(method = "fdr")
pwc %>% filter(p < 0.05)

## # A tibble: 2 x 11
##   Protein Day   .y.   group1 group2    n1    n2 statistic     p p.adj
##   <chr>   <fct> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 Potato  6     logN  33     SK11       4     4         0 0.029 0.200
## 2 Potato  6     logN  SK11   Wg2        4     4        16 0.029 0.200
## # i 1 more variable: p.adj.signif <chr>

3.5.5 VISUALIZE

All CFU data

# load data 
load("R_objects/stat_test_data_growth.RData")

# Create subset
dat.clean <- dat %>% filter(!is.na(CFU) & Day != 0) %>%
  mutate(Day = as.factor(Day))

# Set names of variables

PREDICTOR <- c("Protein","Strain","Day")
OUTCOME <- "logN"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 31 x 7
##    Strain Protein     Day   variable     n  mean    sd
##    <chr>  <chr>       <fct> <fct>    <dbl> <dbl> <dbl>
##  1 33     Casein      1     logN         3  9.21 0.056
##  2 33     Casein      6     logN         2  4.11 0.038
##  3 SK11   Casein      1     logN         3  7.79 0.077
##  4 SK11   Casein      6     logN         3  6.17 0.785
##  5 Wg2    Casein      1     logN         3  9.13 0.019
##  6 Wg2    Casein      6     logN         2  4.46 0.025
##  7 pAK80  Casein      1     logN         3  7.38 0.075
##  8 pAK80  Casein      6     logN         4  7.30 1.02 
##  9 33     Hydrolysate 6     logN         2  2.20 3.12 
## 10 SK11   Hydrolysate 1     logN         2  7.27 0.09 
## # i 21 more rows

#Name order in plot
dat.clean$Day <- as.factor(dat.clean$Day)
dat.clean$Strain <- factor(dat.clean$Strain, levels = c("pAK80", "SK11", "Wg2", "33"))
dat.clean$Protein <- factor(dat.clean$Protein, levels = c("Casein", "Potato", "Hydrolysate", "Water"))

#Plot - boxplot
bxp.logN.sub <- dat.clean %>%
  filter(., Protein  %in% c("Casein", "Potato", "Water")) %>%
  filter (., Day %in% c(1, 2, 6)) %>%
  ggboxplot(.,
            x = "Day",
            y = "logN",
            color = "Strain", 
            facet.by = "Protein")+
  labs(x="Time [days]", y="log10(CFU/mL)")+
  scale_y_continuous(limits = c(0, 10))+
  geom_hline(yintercept=6.8, linetype="dashed", color = "#7d7d7dff")+
  facet_wrap(~Protein, 
             nrow=1, 
             labeller = as_labeller(c("Casein"="CCDM", 
                                      "Potato"="PCDM", 
                                      "Water"="CDM")))+
  scale_color_manual(values = pal_rickandmorty("schwifty")(5)[2:5], 
                    labels=c("MS22418 (PrtP-)", 
                             "MS22421 (PrtP/SK11)", 
                             "MS22425 (PrtP/Wg2)", 
                             "MS22427 (PrtP/MS22333)"))+ 

  theme_light()+
  theme(legend.position="bottom")

bxp.logN.sub 

#Save plot

ggsave(filename = "plots/bxp.logN_all_sub.svg", plot = bxp.logN.sub, width = 16.5, height = 10, units =c("cm"), dpi = 600)

Figure S4

4.1.1 LOAD

dat <- read_excel("./data/Data_collection.xlsx", sheet="Organic_compounds", na = "NA")

4.1.2 LOOK

# Take a glimpse
glimpse(dat)

## Rows: 93
## Columns: 12
## $ `First Injection` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1~
## $ Sample            <chr> "pAK80-P1-D1", "pAK80-P2-D1", "pAK80-P3-D1", "pAK80-~
## $ Strain            <chr> "pAK80", "pAK80", "pAK80", "pAK80", "pAK80", "pAK80"~
## $ Protein           <chr> "Potato", "Potato", "Potato", "Casein", "Casein", "C~
## $ Day               <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
## $ Lactose           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
## $ Glucose           <dbl> 26.2580, 26.7594, 27.0690, 30.0366, 30.3484, 30.0170~
## $ Pyruvate          <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.00~
## $ Lactate           <dbl> 6.7332, 6.8342, 6.8214, 1.6452, 1.8118, 1.7192, 13.7~
## $ Formate           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
## $ Acetate           <dbl> 1.5994, 1.6274, 1.6640, 1.4928, 1.3780, 1.2972, 1.68~
## $ Ethanol           <dbl> 3.8348, 3.1974, 2.7324, 3.0098, 3.3522, 2.8502, 2.93~

# Make a explorative summary 
skimr::skim(dat)

Variable type: character

Variable type: numeric

4.1.3 CLEAN

#Reorganization of data
dat <- pivot_longer(dat, cols=6:12, names_to="compound", values_to = "concentration")

# Remove column
out.var <- c("First Injection")

dat <- dat %>% select(-one_of(out.var))

# Change variables
dat <- dat %>% mutate(Day = as.factor(Day))

# Look at cleaned data
skimr::skim(dat)

Variable type: character

Variable type: factor

Variable type: numeric

4.1.4 SAVE

# Save cleaned data
save(dat, file = "R_objects/stat_test_data_OA.RData")

# clear the environment and release memory
rm(list = ls(all.names = TRUE)[ls(all.names = TRUE) != "params"])
invisible(gc())

4.2.1 CATEGORICAL VARIABLES

Chi squared test of independence.

# load data 
load("R_objects/stat_test_data_OA.RData")

# This converts all character columns to factor (if only specific columns should be converted replace "where(is_character)" with "c(<col1>,<col2>,...)"
dat <- dat %>% mutate(across(where(is_character), as_factor))


# Create vector of categorical variables with more than two categories and less than half of the number of samples.
cat_var <- dat %>% 
  select_if(is.factor) %>% 
  select(where(~n_distinct(.) > 1)) %>% 
  select(where(~n_distinct(.) < (nrow(dat)/2))) %>%
  colnames()

length(cat_var) # If less than 2 variables the following code should not be run

## [1] 5

# Test the variables pairwise 
cat_var %>% 
  combn(2) %>% 
  t() %>% 
  as_tibble() %>% 
  rowwise %>% 
  mutate(chisq_test = list(
    table(dat[[V1]], dat[[V2]]) %>% chisq.test()
    ),
    chisq_pval = chisq_test$p.value
    )

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## i In argument: `chisq_test = list(table(dat[[V1]], dat[[V2]]) %>%
##   chisq.test())`.
## i In row 1.
## Caused by warning in `chisq.test()`:
## ! Chi-squared approximation may be incorrect
## i Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.

## # A tibble: 10 x 4
## # Rowwise: 
##    V1      V2       chisq_test chisq_pval
##    <chr>   <chr>    <list>          <dbl>
##  1 Sample  Strain   <htest>     0        
##  2 Sample  Protein  <htest>     4.88e- 85
##  3 Sample  Day      <htest>     1.76e-167
##  4 Sample  compound <htest>     1   e+  0
##  5 Strain  Protein  <htest>     9.90e-  1
##  6 Strain  Day      <htest>     5.72e-  6
##  7 Strain  compound <htest>     1   e+  0
##  8 Protein Day      <htest>     8.92e-  1
##  9 Protein compound <htest>     1   e+  0
## 10 Day     compound <htest>     1   e+  0

4.2.2 QUANTITATIVE VARIABLES

# create vector with the relevant variables 
CON.VARS <- dat %>% select(where(is.numeric)) %>% colnames()

CON.VARS

## [1] "concentration"

4.3.1 PREPARE VARIABLES

# Set names of variables
PREDICTOR <- c("Protein","Strain", "Day")
OUTCOME <- "concentration"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 28 x 7
##    Strain Protein Day   variable          n   mean    sd
##    <chr>  <chr>   <fct> <fct>         <dbl>  <dbl> <dbl>
##  1 33     Casein  1     concentration     3  0     0    
##  2 33     Casein  2     concentration     3  0.002 0.004
##  3 33     Casein  6     concentration     4  0.005 0.005
##  4 No     Casein  1     concentration     3 30.6   0.13 
##  5 No     Casein  6     concentration     4 29.4   0.592
##  6 SK11   Casein  1     concentration     3 29.3   1.10 
##  7 SK11   Casein  2     concentration     3 27.8   0.212
##  8 SK11   Casein  6     concentration     4 22.3   0.382
##  9 Wg2    Casein  1     concentration     3  0     0    
## 10 Wg2    Casein  2     concentration     3  0.009 0.007
## # i 18 more rows

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

4.3.2 PRELIMINARY TEST

Identify outliers

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 5 x 8
##   Strain Protein Day   Sample      compound concentration is.outlier is.extreme
##   <chr>  <chr>   <fct> <chr>       <chr>            <dbl> <lgl>      <lgl>     
## 1 No     Casein  6     No-C1-D6    Glucose          30.2  TRUE       FALSE     
## 2 Wg2    Casein  6     Wg2-C3-D6   Glucose          29.8  TRUE       FALSE     
## 3 pAK80  Casein  6     pAK80-C4-D6 Glucose          16.7  TRUE       FALSE     
## 4 No     Potato  6     No-P4-D6    Glucose           0    TRUE       FALSE     
## 5 pAK80  Potato  6     pAK80-P4-D6 Glucose           7.92 TRUE       FALSE

Look at outliers
Outliers are derfined with a coefficient (coef=1.5), which specify how far the outlier is from the edge of their boxes.

Extreme outliers are defined with a coef=3.0. However, two outliers appear to be extreme in the boxplot (Samples Wg2-C3-D6 and No-P4-D6) why they are tested with different coefficients. Both points are outliers with high coefficient values that are close to 3. Therefore, we defined outliers as extreme with coef above 2.59. The extreme outliers are removed.

#Test coef for outlier Wg2-C3-D6 
dat.clean %>%
  filter(Strain == "Wg2" & Protein == "Casein" & Day == 6) %>%
  pull(concentration) %>% is_outlier(coef = 2.995)

## [1] FALSE FALSE  TRUE FALSE

#Test coef for outlier No-P4-D6
dat.clean %>%
  filter(Strain == "No" & Protein == "Potato" & Day == 6) %>%
  pull(concentration) %>% is_outlier(coef = 2.88)

## [1] FALSE FALSE FALSE  TRUE

#Test coef for outlier pAK80-C4-D6
dat.clean %>%
  filter(Strain == "pAK80" & Protein == "Casein" & Day == 6) %>%
  pull(concentration) %>% is_outlier(coef = 2.59)

## [1] FALSE FALSE FALSE FALSE

#Remove extreme outliers

extreme_outliers <- c("Wg2-C3-D6", "No-P4-D6")

dat.clean <- subset(dat.clean, !Sample %in% extreme_outliers) 

# Create plot without extreme outliers
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

Check normality

QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.The Shapiro-Wilk test is performed on all data groups when possible with n > 2 and only with unique values.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)

# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 x 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.468 1.50e-16

# Compute Shapiro-Wilk test of normality for each data group
sets <- dat.clean %>%
  group_by(Day, Protein, Strain) %>% 
  summarise(n_unique = n_distinct(concentration),
            n = n()) %>%
  filter(n > 2, n_unique > 1) %>% 
  transmute(day_protein_strain = paste(Day, Protein, Strain))

## `summarise()` has grouped output by 'Day', 'Protein'. You can override using
## the `.groups` argument.

dat.clean %>% 
  mutate(day_protein_strain = paste(Day, Protein, Strain)) %>%
  filter(day_protein_strain %in% sets$day_protein_strain) %>%
  group_by(Day, Protein, Strain) %>% 
  shapiro_test(vars = "concentration") %>% 
  adjust_pvalue(method = "fdr") %>%
  filter(p.adj < 0.05)

## # A tibble: 4 x 7
##   Strain Protein Day   variable      statistic       p  p.adj
##   <chr>  <chr>   <fct> <chr>             <dbl>   <dbl>  <dbl>
## 1 33     Casein  2     concentration     0.75  0       0     
## 2 Wg2    Casein  2     concentration     0.75  0       0     
## 3 No     Casein  6     concentration     0.665 0.00415 0.0259
## 4 Wg2    Casein  6     concentration     0.75  0       0

Test homogneity of variance assumption

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 x 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1    27    63      1.36 0.156

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

4.3.3 PERFORM TEST

The data are not normal distributed. The Shapiro-Wilk normality test found that at least one group was not normal distributed. As a result, for statistical data analysis, the Kruskal-Wallis one-way analysis of variance test is used, followed by a Wilcoxon test for pairwise group comparisons.

# for each protein source
res.protein <- dat.clean %>% 
  group_by(Day, Strain) %>%
  kruskal_test(concentration ~ Protein) %>%
  adjust_pvalue(method = "fdr")
res.protein %>% filter(p < 0.05)

## # A tibble: 11 x 9
##    Strain Day   .y.               n statistic    df      p method         p.adj
##    <chr>  <fct> <chr>         <int>     <dbl> <int>  <dbl> <chr>          <dbl>
##  1 33     1     concentration     6      4.35     1 0.0369 Kruskal-Wallis 0.063
##  2 SK11   1     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.063
##  3 Wg2    1     concentration     6      4.35     1 0.0369 Kruskal-Wallis 0.063
##  4 pAK80  1     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.063
##  5 33     2     concentration     6      3.97     1 0.0463 Kruskal-Wallis 0.063
##  6 Wg2    2     concentration     6      3.97     1 0.0463 Kruskal-Wallis 0.063
##  7 pAK80  2     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.063
##  8 33     6     concentration     8      5.40     1 0.0202 Kruskal-Wallis 0.063
##  9 SK11   6     concentration     8      5.33     1 0.0209 Kruskal-Wallis 0.063
## 10 Wg2    6     concentration     7      4.58     1 0.0323 Kruskal-Wallis 0.063
## 11 pAK80  6     concentration     8      5.33     1 0.0209 Kruskal-Wallis 0.063

# for each strain
res.strain <- dat.clean %>% 
  group_by(Day, Protein) %>%
  kruskal_test(concentration ~ Strain) %>%
  adjust_pvalue(method = "fdr")
res.strain %>% filter(p < 0.05)

## # A tibble: 6 x 9
##   Protein Day   .y.               n statistic    df       p method         p.adj
##   <chr>   <fct> <chr>         <int>     <dbl> <int>   <dbl> <chr>          <dbl>
## 1 Casein  1     concentration    15     13.2      4 0.0104  Kruskal-Wal~ 0.0182 
## 2 Potato  1     concentration    15     12.8      4 0.0121  Kruskal-Wal~ 0.0182 
## 3 Casein  2     concentration    12      9.84     3 0.02    Kruskal-Wal~ 0.024  
## 4 Potato  2     concentration    11      9.41     3 0.0243  Kruskal-Wal~ 0.0243 
## 5 Casein  6     concentration    19     16.1      4 0.00289 Kruskal-Wal~ 0.00867
## 6 Potato  6     concentration    19     16.3      4 0.00265 Kruskal-Wal~ 0.00867

pwc <- dat.clean %>%
  filter(Day == 6) %>%
  group_by(Day, Protein) %>%
  pairwise_wilcox_test(concentration ~ Strain) %>%
  adjust_pvalue(method = "fdr")
pwc %>% filter(p < 0.05)

## # A tibble: 10 x 11
##    Protein Day   .y.           group1 group2    n1    n2 statistic     p p.adj
##    <chr>   <fct> <chr>         <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
##  1 Casein  6     concentration 33     No         4     4         0 0.029 0.058
##  2 Casein  6     concentration 33     pAK80      4     4         0 0.029 0.058
##  3 Casein  6     concentration 33     SK11       4     4         0 0.029 0.058
##  4 Casein  6     concentration No     pAK80      4     4        16 0.029 0.058
##  5 Casein  6     concentration No     SK11       4     4        16 0.029 0.058
##  6 Potato  6     concentration 33     pAK80      4     4         0 0.029 0.058
##  7 Potato  6     concentration 33     SK11       4     4         0 0.029 0.058
##  8 Potato  6     concentration pAK80  SK11       4     4        16 0.029 0.058
##  9 Potato  6     concentration pAK80  Wg2        4     4        16 0.029 0.058
## 10 Potato  6     concentration SK11   Wg2        4     4        16 0.029 0.058
## # i 1 more variable: p.adj.signif <chr>

4.3.4 VISUALIZE

#Name order in plot
dat.clean$Day <- as.factor(dat.clean$Day)
dat.clean$Strain <- factor(dat.clean$Strain, 
                           levels = c("No", "pAK80", "SK11", "Wg2", "33"), 
                           labels =c ("Non-inoculated", 
                                      "MS22418 (PrtP-)", 
                                      "MS22421 (PrtP/SK11)", 
                                      "MS22425 (PrtP/Wg2)", 
                                      "MS22427 (PrtP/MS22333)")
                           )
dat.clean$Protein <- factor(dat.clean$Protein, levels = c("Casein", "Potato"))

#Plot - boxplot
bxp.glucose <- dat.clean %>%
  filter(., Protein  %in% c("Casein", "Potato")) %>%
  filter (., Day %in% c(1, 2, 6)) %>%
  ggboxplot(.,
            x = "Day",
            y = "concentration",
            color = "Strain", 
            facet.by = "Protein")+
  labs(x="Time [days]", y="Concentration [mM]")+
  scale_y_continuous()+
  facet_wrap(~Protein, 
             nrow=1, 
             labeller = as_labeller(c("Casein"="CCDM", 
                                      "Potato"="PCDM")
                                    )
             )+
  scale_color_manual(values = c("Non-inoculated"="#E89242FF", 
                                  "MS22418 (PrtP-)"="#82491EFF", 
                                  "MS22421 (PrtP/SK11)"="#24325FFF", 
                                  "MS22425 (PrtP/Wg2)"="#B7E4F9FF", 
                                  "MS22427 (PrtP/MS22333)"="#FB6467FF")
                           )+
  geom_hline(yintercept=30.0, linetype="dashed", color = "#7d7d7dff")+
  theme_light()+
  theme(legend.position="bottom")

save(bxp.glucose, file = "R_objects/Fig1B.RData")

4.5.1 PRELIMINARY TEST

Identify outliers

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 3 x 8
##   Strain Protein Day   Sample      compound concentration is.outlier is.extreme
##   <chr>  <chr>   <fct> <chr>       <chr>            <dbl> <lgl>      <lgl>     
## 1 Wg2    Casein  6     Wg2-C3-D6   Lactate           0.15 TRUE       FALSE     
## 2 pAK80  Casein  6     pAK80-C4-D6 Lactate          20.9  TRUE       FALSE     
## 3 No     Potato  6     No-P4-D6    Lactate          47.9  TRUE       FALSE

Look at outliers

Outliers are derfined with a coefficient (coef=1.5), which specify how far the outliers are from the edge of their box.

Points are extreme outliers if they have high coefficient values that are close to 3. Therefore, we defined outliers as extreme with coef above 2.59. The extreme outliers are removed.

#Test coef for outlier Wg2-C3-D6 
dat.clean %>%
  filter(Strain == "Wg2" & Protein == "Casein" & Day == 6) %>%
  pull(concentration) %>% is_outlier(coef = 2.59)

## [1] FALSE FALSE  TRUE FALSE

#Test coef for outlier No-P4-D6
dat.clean %>%
  filter(Strain == "No" & Protein == "Potato" & Day == 6) %>%
  pull(concentration) %>% is_outlier(coef = 2.98)

## [1] FALSE FALSE FALSE  TRUE

#Test coef for outlier pAK80-C4-D6
dat.clean %>%
  filter(Strain == "pAK80" & Protein == "Casein" & Day == 6) %>%
  pull(concentration) %>% is_outlier(coef = 2.55)

## [1] FALSE FALSE FALSE FALSE

#Remove extreme outliers

extreme_outliers <- c("Wg2-C3-D6", "No-P4-D6")

dat.clean <- subset(dat.clean, !Sample %in% extreme_outliers) 

# Create plot without extreme outliers
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

Check normality

QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.The Shapiro-Wilk test is performed on all data groups when possible with n > 2 and only with unique values.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)

# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 x 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.706 3.32e-12

# Compute Shapiro-Wilk test of normality for each data group
sets <- dat.clean %>%
  group_by(Day, Protein, Strain) %>% 
  summarise(n_unique = n_distinct(concentration),
            n = n()) %>%
  filter(n > 2, n_unique > 1) %>% 
  transmute(day_protein_strain = paste(Day, Protein, Strain))

## `summarise()` has grouped output by 'Day', 'Protein'. You can override using
## the `.groups` argument.

dat.clean %>% 
  mutate(day_protein_strain = paste(Day, Protein, Strain)) %>%
  filter(day_protein_strain %in% sets$day_protein_strain) %>%
  group_by(Day, Protein, Strain) %>% 
  shapiro_test(vars = "concentration") %>% 
  adjust_pvalue(method = "fdr") %>%
  filter(p.adj < 0.05)

## # A tibble: 0 x 7
## # i 7 variables: Strain <chr>, Protein <chr>, Day <fct>, variable <chr>,
## #   statistic <dbl>, p <dbl>, p.adj <dbl>

Test homogneity of variance assumption

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 x 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1    27    63     0.956 0.537

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

4.5.2 PERFORM TEST

The data are not normal distributed, but normality is found within each group by the Shapiro-Wilk normality test. However, the same data groups of the glucose data were not normal distributed.

The same statistical analysis is applied for the two data sets, including a Kruskal-Wallis test and a Wilcoxon pairwise test.

# for each protein source
res.protein <- dat.clean %>% 
  group_by(Day, Strain) %>%
  kruskal_test(concentration ~ Protein) %>%
  adjust_pvalue(method = "fdr")
res.protein %>% filter(p < 0.05)

## # A tibble: 10 x 9
##    Strain Day   .y.               n statistic    df      p method          p.adj
##    <chr>  <fct> <chr>         <int>     <dbl> <int>  <dbl> <chr>           <dbl>
##  1 33     1     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.0693
##  2 SK11   1     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.0693
##  3 Wg2    1     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.0693
##  4 pAK80  1     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.0693
##  5 Wg2    2     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.0693
##  6 pAK80  2     concentration     6      3.86     1 0.0495 Kruskal-Wallis 0.0693
##  7 33     6     concentration     8      5.33     1 0.0209 Kruskal-Wallis 0.0693
##  8 SK11   6     concentration     8      5.33     1 0.0209 Kruskal-Wallis 0.0693
##  9 Wg2    6     concentration     7      4.5      1 0.0339 Kruskal-Wallis 0.0693
## 10 pAK80  6     concentration     8      5.33     1 0.0209 Kruskal-Wallis 0.0693

# for each strain
res.strain <- dat.clean %>% 
  group_by(Day, Protein) %>%
  kruskal_test(concentration ~ Strain) %>%
  adjust_pvalue(method = "fdr")
res.strain %>% filter(p < 0.05)

## # A tibble: 6 x 9
##   Protein Day   .y.               n statistic    df       p method         p.adj
##   <chr>   <fct> <chr>         <int>     <dbl> <int>   <dbl> <chr>          <dbl>
## 1 Casein  1     concentration    15     12.8      4 0.0121  Kruskal-Wal~ 0.0182 
## 2 Potato  1     concentration    15     12.9      4 0.0118  Kruskal-Wal~ 0.0182 
## 3 Casein  2     concentration    12      9.67     3 0.0216  Kruskal-Wal~ 0.0243 
## 4 Potato  2     concentration    11      9.41     3 0.0243  Kruskal-Wal~ 0.0243 
## 5 Casein  6     concentration    19     16.1      4 0.00293 Kruskal-Wal~ 0.00879
## 6 Potato  6     concentration    19     16.3      4 0.00263 Kruskal-Wal~ 0.00879

pwc <- dat.clean %>%
  filter(Day == 6) %>%
  group_by(Day, Protein) %>%
  pairwise_wilcox_test(concentration ~ Strain) %>%
  adjust_pvalue(method = "fdr")
pwc %>% filter(p < 0.05)

## # A tibble: 10 x 11
##    Protein Day   .y.           group1 group2    n1    n2 statistic     p p.adj
##    <chr>   <fct> <chr>         <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
##  1 Casein  6     concentration 33     No         4     4        16 0.029 0.058
##  2 Casein  6     concentration 33     pAK80      4     4        16 0.029 0.058
##  3 Casein  6     concentration 33     SK11       4     4        16 0.029 0.058
##  4 Casein  6     concentration No     pAK80      4     4         0 0.029 0.058
##  5 Casein  6     concentration No     SK11       4     4         0 0.029 0.058
##  6 Potato  6     concentration 33     pAK80      4     4        16 0.029 0.058
##  7 Potato  6     concentration 33     SK11       4     4        16 0.029 0.058
##  8 Potato  6     concentration pAK80  SK11       4     4         0 0.029 0.058
##  9 Potato  6     concentration pAK80  Wg2        4     4         0 0.029 0.058
## 10 Potato  6     concentration SK11   Wg2        4     4         0 0.029 0.058
## # i 1 more variable: p.adj.signif <chr>

4.5.3 VISUALIZE

#Name order in plot
dat.clean$Day <- as.factor(dat.clean$Day)
dat.clean$Strain <- factor(dat.clean$Strain, 
                           levels = c("No", "pAK80", "SK11", "Wg2", "33"),
                           labels =c ("Non-inoculated", 
                                      "MS22418 (PrtP-)", 
                                      "MS22421 (PrtP/SK11)", 
                                      "MS22425 (PrtP/Wg2)", 
                                      "MS22427 (PrtP/MS22333)")
                           )
dat.clean$Protein <- factor(dat.clean$Protein, levels = c("Casein", "Potato"))

#Boxplot
bxp.lactate <- dat.clean %>%
  filter(., Protein  %in% c("Casein", "Potato")) %>%
  filter (., Day %in% c(1, 2, 6)) %>%
  ggboxplot(.,
            x = "Day",
            y = "concentration",
            color = "Strain", 
            facet.by = "Protein")+
  labs(x="Time [days]", y="Concentration [mM]")+
  scale_y_continuous()+
  facet_wrap(~Protein, 
             nrow=1, 
             labeller = as_labeller(c("Casein"="CCDM", 
                                      "Potato"="PCDM")
                                    )
             )+
  scale_color_manual(values = c("Non-inoculated"="#E89242FF", 
                                  "MS22418 (PrtP-)"="#82491EFF", 
                                  "MS22421 (PrtP/SK11)"="#24325FFF", 
                                  "MS22425 (PrtP/Wg2)"="#B7E4F9FF", 
                                  "MS22427 (PrtP/MS22333)"="#FB6467FF")
                           )+

  theme_light()+
  theme(legend.position="bottom")

save(bxp.lactate, file = "R_objects/Fig1C.RData")

5.1.1 LOAD

dat <- read_excel("./data/Data_collection.xlsx", sheet="AA", na = "NA")

5.1.2 LOOK

# Take a glimpse
glimpse(dat)

## Rows: 94
## Columns: 27
## $ Sample_ID         <chr> "Wg2-C1-D1", "Wg2-C2-D1", "Wg2-C3-D1", "Sk11-C1-D1",~
## $ Sample            <chr> "DC1-1", "DC1-2", "DC1-3", "CC1-1", "CC1-2", "CC1-3"~
## $ `First Injection` <chr> "50Lise", "51Lise", "52Lise", "56Lise", "57Lise", "5~
## $ `Sample no`       <chr> "C1", "C2", "C3", "C1", "C2", "C3", "C1", "C2", "C3"~
## $ Day               <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2~
## $ Protein           <chr> "Casein", "Casein", "Casein", "Casein", "Casein", "C~
## $ Strain            <chr> "Wg2", "Wg2", "Wg2", "SK11", "SK11", "SK11", "pAK80"~
## $ conc              <chr> "uM", "uM", "uM", "uM", "uM", "uM", "uM", "uM", "uM"~
## $ His               <dbl> 18.4537472, 17.1673140, 18.5421378, 0.7142539, 0.794~
## $ Asn               <dbl> 13.546499, 12.328580, 15.030006, 1.541399, 2.227930,~
## $ Ser               <dbl> 6.097704, 6.053979, 6.150434, 4.579292, 5.554308, 5.~
## $ `Gln/Arg`         <dbl> 3.3068133, 2.9704246, 3.1039723, 0.2580126, 0.286869~
## $ Gly               <dbl> 59.60045, 58.21350, 57.65550, 66.11592, 83.48793, 82~
## $ Asp               <dbl> 2.3255601, 2.5187825, 2.4126894, 0.5196090, 0.770297~
## $ Glu               <dbl> 31.6012219, 31.9529499, 29.6124309, 0.7158399, 1.459~
## $ Thr               <dbl> 3.6553179, 3.2991876, 3.7922678, 0.2722407, 0.605377~
## $ Ala               <dbl> 22.832045, 22.878213, 26.319410, 1.558777, 1.925683,~
## $ Pro               <dbl> 318.4714880, 309.0975224, 314.1835222, 4.9019494, 5.~
## $ Cys               <dbl> 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.000000~
## $ Lys               <dbl> 5.526448, 6.202491, 4.736371, 2.165367, 2.188678, 2.~
## $ Tyr               <dbl> 95.12991, 90.49597, 96.35402, 15.21761, 16.79873, 17~
## $ Met               <dbl> 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.~
## $ Val               <dbl> 17.711121, 17.242679, 18.731475, 3.099217, 3.132582,~
## $ Ile               <dbl> 24.26801, 24.62159, 25.36100, 22.05858, 24.04160, 24~
## $ Leu               <dbl> 25.0888213, 23.7935900, 25.8345528, 0.3067792, 0.511~
## $ Phe               <dbl> 51.791264, 48.307623, 51.469286, 3.014791, 3.351971,~
## $ Trp               <dbl> 36.56362, 35.95126, 37.81154, 26.58388, 29.87836, 30~

# Make a explorative summary 
skimr::skim(dat)

Variable type: character

Variable type: numeric

5.1.3 CLEAN DATA

# Reorganization of data
dat <-  pivot_longer(dat, cols=9:27, names_to="aa", values_to = "concentration")

# Change variables
dat <- dat %>% mutate(Day = as.factor(Day))

# Look at cleaned data
skimr::skim(dat)

Variable type: character

Variable type: factor

Variable type: numeric

glimpse(dat)

## Rows: 1,786
## Columns: 10
## $ Sample_ID         <chr> "Wg2-C1-D1", "Wg2-C1-D1", "Wg2-C1-D1", "Wg2-C1-D1", ~
## $ Sample            <chr> "DC1-1", "DC1-1", "DC1-1", "DC1-1", "DC1-1", "DC1-1"~
## $ `First Injection` <chr> "50Lise", "50Lise", "50Lise", "50Lise", "50Lise", "5~
## $ `Sample no`       <chr> "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1"~
## $ Day               <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
## $ Protein           <chr> "Casein", "Casein", "Casein", "Casein", "Casein", "C~
## $ Strain            <chr> "Wg2", "Wg2", "Wg2", "Wg2", "Wg2", "Wg2", "Wg2", "Wg~
## $ conc              <chr> "uM", "uM", "uM", "uM", "uM", "uM", "uM", "uM", "uM"~
## $ aa                <chr> "His", "Asn", "Ser", "Gln/Arg", "Gly", "Asp", "Glu",~
## $ concentration     <dbl> 18.453747, 13.546499, 6.097704, 3.306813, 59.600447,~

5.1.4 SAVE

# Save cleaned data
save(dat, file = "R_objects/stat_test_data_AA.RData")

# clear the environment and release memory
rm(list = ls(all.names = TRUE)[ls(all.names = TRUE) != "params"])
invisible(gc())

5.2.1 CATEGORICAL VARIABLES

Chi squared test of independence

# load data 
load("R_objects/stat_test_data_AA.RData")

# This converts all character columns to factor (if only specific columns should be converted replace "where(is_character)" with "c(<col1>,<col2>,...)"
dat <- dat %>% mutate(across(where(is_character), as_factor))


# Create vector of categorical variables with more than two categories and less than half of the number of samples.
cat_var <- dat %>% 
  select_if(is.factor) %>% 
  select(where(~n_distinct(.) > 1)) %>% 
  select(where(~n_distinct(.) < (nrow(dat)/2))) %>%
  colnames()

length(cat_var) # If less than 2 variables the following code should not be run

## [1] 8

# Test the variables pairwise 
cat_var %>% 
  combn(2) %>% 
  t() %>% 
  as_tibble() %>% 
  rowwise %>% 
  mutate(chisq_test = list(
    table(dat[[V1]], dat[[V2]]) %>% chisq.test()
    ),
    chisq_pval = chisq_test$p.value
    )

## Warning: There were 15 warnings in `mutate()`.
## The first warning was:
## i In argument: `chisq_test = list(table(dat[[V1]], dat[[V2]]) %>%
##   chisq.test())`.
## i In row 1.
## Caused by warning in `chisq.test()`:
## ! Chi-squared approximation may be incorrect
## i Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 14 remaining warnings.

## # A tibble: 28 x 4
## # Rowwise: 
##    V1        V2              chisq_test chisq_pval
##    <chr>     <chr>           <list>          <dbl>
##  1 Sample_ID Sample          <htest>     0        
##  2 Sample_ID First Injection <htest>     0        
##  3 Sample_ID Sample no       <htest>     0        
##  4 Sample_ID Day             <htest>     0        
##  5 Sample_ID Protein         <htest>     3.58e-311
##  6 Sample_ID Strain          <htest>     0        
##  7 Sample_ID aa              <htest>     1   e+  0
##  8 Sample    First Injection <htest>     0        
##  9 Sample    Sample no       <htest>     0        
## 10 Sample    Day             <htest>     0        
## # i 18 more rows

5.2.2 QUANTITATIVE VARIABLES

# create vector with the relevant variables 
CON.VARS <- dat %>% select(where(is.numeric)) %>% colnames()

CON.VARS

## [1] "concentration"

5.3.1 DATA

# load data 
load("R_objects/stat_test_data_AA.RData")

# data subset

dat.clean <- dat %>%
  filter(., Protein == "Casein", !is.na(concentration))

# Set names of variables
PREDICTOR <- c("Strain", "Day", "aa")
OUTCOME <- "concentration"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 266 x 7
##    Day   Strain  aa      variable          n   mean    sd
##    <fct> <chr>   <chr>   <fct>         <dbl>  <dbl> <dbl>
##  1 1     MS22333 Ala     concentration     3 41.9   1.27 
##  2 1     MS22333 Asn     concentration     3 32.7   1.88 
##  3 1     MS22333 Asp     concentration     3  4.69  0.286
##  4 1     MS22333 Cys     concentration     3  0.394 0.161
##  5 1     MS22333 Gln/Arg concentration     3  3.07  0.123
##  6 1     MS22333 Glu     concentration     3 81.1   1.42 
##  7 1     MS22333 Gly     concentration     3 74.3   2.42 
##  8 1     MS22333 His     concentration     3 49.9   1.40 
##  9 1     MS22333 Ile     concentration     3 22.9   0.925
## 10 1     MS22333 Leu     concentration     3 28.9   0.328
## # i 256 more rows

Identify outliers

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 53 x 12
##    Day   Strain  aa      Sample_ID Sample `First Injection` `Sample no` Protein
##    <fct> <chr>   <chr>   <chr>     <chr>  <chr>             <chr>       <chr>  
##  1 6     MS22333 Asn     33-C4-D6  EC6-4  32Lise            C4          Casein 
##  2 6     MS22333 Gln/Arg 33-C3-D6  EC6-3  31Lise            C3          Casein 
##  3 6     MS22333 His     33-C3-D6  EC6-3  31Lise            C3          Casein 
##  4 6     MS22333 Met     33-C4-D6  EC6-4  32Lise            C4          Casein 
##  5 6     MS22333 Phe     33-C3-D6  EC6-3  31Lise            C3          Casein 
##  6 6     MS22333 Thr     33-C3-D6  EC6-3  31Lise            C3          Casein 
##  7 6     MS22333 Tyr     33-C4-D6  EC6-4  32Lise            C4          Casein 
##  8 6     MS22333 Val     33-C4-D6  EC6-4  32Lise            C4          Casein 
##  9 6     No      Ala     No-C3-D6  AC6-3  39Lise            C3          Casein 
## 10 6     No      Asn     No-C2-D6  AC6-2  38Lise            C2          Casein 
## # i 43 more rows
## # i 4 more variables: conc <chr>, concentration <dbl>, is.outlier <lgl>,
## #   is.extreme <lgl>

###VISUALIZE

#Name order in plot
dat.clean$Day <- as.factor(dat.clean$Day)
dat.clean$Strain <- factor(dat.clean$Strain, 
                           levels = c("No", "pAK80", "SK11", "Wg2", "MS22333"), 
                           labels =c ("Non-inoculated", 
                                      "MS22418 (PrtP-)", 
                                      "MS22421 (PrtP/SK11)", 
                                      "MS22425 (PrtP/Wg2)", 
                                      "MS22427 (PrtP/MS22333)"))
dat.clean$Protein <- factor(dat.clean$Protein, levels = c("Casein", "Potato"))

#Boxplots

bxp.AA.CCDM <- dat.clean %>%
  ggboxplot(.,
            x = "Day",
            y = "concentration",
            color = "Strain",
            facet.by = "aa")+
  labs(title="CCDM", x="Time [days]", y="Amino acid concentration [uM]")+
  scale_y_continuous()+
  facet_wrap(~aa, 
             nrow=5,
             scales = "free_y")+
  scale_color_manual(values = c("Non-inoculated"="#E89242FF", 
                                  "MS22418 (PrtP-)"="#82491EFF", 
                                  "MS22421 (PrtP/SK11)"="#24325FFF", 
                                  "MS22425 (PrtP/Wg2)"="#B7E4F9FF", 
                                  "MS22427 (PrtP/MS22333)"="#FB6467FF")
                           )+
  theme_light()+
  theme(legend.position="bottom")

#Save plot
ggsave(filename = "plots/bxp.AA.CCDM.svg", plot = bxp.AA.CCDM, width = 16.5, height = 20, units =c("cm"), dpi = 600)

Figure S5A

5.4.1 DATA

# load data 
load("R_objects/stat_test_data_AA.RData")

# data subset

dat.clean <- dat %>%
  filter(., Protein == "Potato", !is.na(concentration))

# Set names of variables
PREDICTOR <- c("Strain", "Day", "aa")
OUTCOME <- "concentration"
SUBJECT <- NULL

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 266 x 7
##    Day   Strain  aa      variable          n   mean    sd
##    <fct> <chr>   <chr>   <fct>         <dbl>  <dbl> <dbl>
##  1 1     MS22333 Ala     concentration     3 10.3   0.274
##  2 1     MS22333 Asn     concentration     3  3.74  1.77 
##  3 1     MS22333 Asp     concentration     3  7.71  0.054
##  4 1     MS22333 Cys     concentration     3  0     0    
##  5 1     MS22333 Gln/Arg concentration     3  9.75  0.105
##  6 1     MS22333 Glu     concentration     3  9.10  0.297
##  7 1     MS22333 Gly     concentration     3 80.8   1.97 
##  8 1     MS22333 His     concentration     3  0.911 0.11 
##  9 1     MS22333 Ile     concentration     3 23.0   0.481
## 10 1     MS22333 Leu     concentration     3  1.74  0.057
## # i 256 more rows

Identify outliers

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 17 x 12
##    Day   Strain  aa      Sample_ID  Sample `First Injection` `Sample no` Protein
##    <fct> <chr>   <chr>   <chr>      <chr>  <chr>             <chr>       <chr>  
##  1 6     MS22333 Ile     33_P1-D6   EP6-1  Lise25            P1          Potato 
##  2 6     MS22333 Leu     33_P1-D6   EP6-1  Lise25            P1          Potato 
##  3 6     MS22333 Val     33_P1-D6   EP6-1  Lise25            P1          Potato 
##  4 6     SK11    Asn     SK11_P1-D6 CP6-1  Lise17            P1          Potato 
##  5 6     SK11    Glu     SK11_P1-D6 CP6-1  Lise17            P1          Potato 
##  6 6     SK11    Gly     SK11_P1-D6 CP6-1  Lise17            P1          Potato 
##  7 6     SK11    Thr     SK11_P4-D6 CP6-4  Lise20            P4          Potato 
##  8 6     SK11    Val     SK11_P1-D6 CP6-1  Lise17            P1          Potato 
##  9 6     Wg2     Gln/Arg Wg2_P3-D6  DP6-3  Lise11            P3          Potato 
## 10 6     Wg2     Ser     Wg2_P1-D6  DP6-1  Lise9             P1          Potato 
## 11 6     Wg2     Thr     Wg2_P1-D6  DP6-1  Lise9             P1          Potato 
## 12 6     Wg2     Tyr     Wg2_P3-D6  DP6-3  Lise11            P3          Potato 
## 13 6     pAK80   Gly     pAK80_P4-~ BP6-4  Lise4             P4          Potato 
## 14 6     pAK80   Leu     pAK80_P3-~ BP6-3  Lise3             P3          Potato 
## 15 6     pAK80   Lys     pAK80_P3-~ BP6-3  Lise3             P3          Potato 
## 16 6     pAK80   Trp     pAK80_P2-~ BP6-2  Lise2             P2          Potato 
## 17 6     pAK80   Val     pAK80_P3-~ BP6-3  Lise3             P3          Potato 
## # i 4 more variables: conc <chr>, concentration <dbl>, is.outlier <lgl>,
## #   is.extreme <lgl>

###VISUALIZE

#Name order in plot
dat.clean$Day <- as.factor(dat.clean$Day)
dat.clean$Strain <- factor(dat.clean$Strain, 
                           levels = c("No", "pAK80", "SK11", "Wg2", "MS22333"), 
                           labels =c ("Non-inoculated", 
                                      "MS22418 (PrtP-)", 
                                      "MS22421 (PrtP/SK11)", 
                                      "MS22425 (PrtP/Wg2)", 
                                      "MS22427 (PrtP/MS22333)"))
dat.clean$Protein <- factor(dat.clean$Protein, levels = c("Casein", "Potato"))

#Boxplots

bxp.AA.PCDM <- dat.clean %>%
  ggboxplot(.,
            x = "Day",
            y = "concentration",
            color = "Strain",
            facet.by = "aa")+
  labs(title="PCDM", x="Time [days]", y="Amino acid concentration [uM]")+
  scale_y_continuous()+
  facet_wrap(~aa, 
             nrow=5,
             scales = "free_y")+
  scale_color_manual(values = c("Non-inoculated"="#E89242FF", 
                                  "MS22418 (PrtP-)"="#82491EFF", 
                                  "MS22421 (PrtP/SK11)"="#24325FFF", 
                                  "MS22425 (PrtP/Wg2)"="#B7E4F9FF", 
                                  "MS22427 (PrtP/MS22333)"="#FB6467FF")
                           )+
  theme_light()+
  theme(legend.position="bottom")

#Save plot

ggsave(filename = "plots/bxp.AA.PCDM.svg", plot = bxp.AA.PCDM, width = 16.5, height = 20, units =c("cm"), dpi = 600)

Figure S5B

5.5.1 NORMALIZED DATA

dat.norm  <- dat %>%
  group_by(Sample) %>%
  mutate(conc_norm = 100*concentration/sum(concentration))

table(is.na(dat.norm$concentration))

## 
## FALSE  TRUE 
##  1767    19

# Extract metadata and aa table
dat.wide <- dat.norm %>%
  filter(!is.na(concentration)) %>%
  pivot_wider(values_from = conc_norm, 
              names_from = aa, 
              id_cols = c(`Sample no`, Strain, Protein, Day, Sample), 
              values_fill = 0)

meta <- dat.wide %>% 
  ungroup() %>%
  select(`Sample no`:Sample) %>% 
  data.frame
row.names(meta) <- meta$Sample

#aa.num defines the aa data in prercent
aa.num <- dat.wide %>%
  ungroup() %>%
  select(His:Trp) %>%
  data.frame
row.names(aa.num) <- meta$Sample

#aa.norm is the calculated z-score used in the heatmap

aa.norm <- aa.num

for (i in 1:ncol(aa.num)){
  x.no <- mean(aa.num[meta$Strain == "No",i])
  aa.norm[,i] <- (aa.num[,i]-x.no)/sd(aa.num[,i])
}

5.5.2 VISUALIZE: Heatmap

#Annotations
ano_df2 =data.frame("Day"=meta$Day, "Medium" = meta$Protein, "Strain" = meta$Strain)
rownames(ano_df2) = rownames(meta)
ano_df2$Day <- as.factor(ano_df2$Day)
ano_df2$Medium <- factor(ano_df2$Medium, levels = c("Casein", "Potato"), labels=c("CCDM", "PCDM"))
ano_df2$Strain <- factor(ano_df2$Strain, levels = c("No", "pAK80", "SK11", "Wg2", "MS22333"), labels=c("A: Non-inoculated", "B: MS22418 (PrtP-)", "C: MS22421 (PrtP/SK11)", "D: MS22425 (PrtP/Wg2)", "E: MS22427 (PrtP/MS22333)"))

annot_colors=list("Medium"=c("CCDM"="#A6EEE6FF", "PCDM"="#FAE48BFF"), "Day" = c("1"= "#d2d3daff", "2"= "#86889bff", "6"="#282931ff"), "Strain" =c("A: Non-inoculated"="#E89242FF", "B: MS22418 (PrtP-)"="#82491EFF", "C: MS22421 (PrtP/SK11)"="#24325FFF", "D: MS22425 (PrtP/Wg2)"="#B7E4F9FF", "E: MS22427 (PrtP/MS22333)"="#FB6467FF"))

#Plot
heatmap <- pheatmap(aa.norm,
         annotation_row = ano_df2, 
         cutree_cols = 3, cutree_rows = 5, border_color = NA,
         color=colorRampPalette(c("navy", "white", "#d50c00ff"))(50), annotation_colors = annot_colors)

#Save
ggsave(filename = "plots/pmap_aa.svg", plot = heatmap,width = 20, height = 32, units = "cm",dpi = 600)

Figure 3

5.6.1 PERMANOVA

# Calculate distances

ed <- vegdist(aa.num, method = "euclidian")

adonis2(ed ~ Day, data = ano_df2)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = ed ~ Day, data = ano_df2)
##          Df SumOfSqs      R2      F Pr(>F)   
## Day       2     5324 0.08926 4.4103  0.005 **
## Residual 90    54328 0.91074                 
## Total    92    59652 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(ed ~ Strain, data = ano_df2)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = ed ~ Strain, data = ano_df2)
##          Df SumOfSqs     R2      F Pr(>F)    
## Strain    4    29063 0.4872 20.902  0.001 ***
## Residual 88    30589 0.5128                  
## Total    92    59652 1.0000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(ed ~ Medium, data = ano_df2)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = ed ~ Medium, data = ano_df2)
##          Df SumOfSqs      R2      F Pr(>F)   
## Medium    1     4415 0.07401 7.2736  0.006 **
## Residual 91    55237 0.92599                 
## Total    92    59652 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(ed ~ Strain*Day*Medium, data = ano_df2)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = ed ~ Strain * Day * Medium, data = ano_df2)
##                   Df SumOfSqs      R2        F Pr(>F)    
## Strain             4    29063 0.48720 123.0467  0.001 ***
## Day                2     5072 0.08502  42.9452  0.001 ***
## Medium             1     4496 0.07536  76.1354  0.001 ***
## Strain:Day         7     3214 0.05388   7.7763  0.001 ***
## Strain:Medium      4    11093 0.18596  46.9646  0.001 ***
## Day:Medium         2     1298 0.02177  10.9943  0.001 ***
## Strain:Day:Medium  7     1579 0.02647   3.8194  0.001 ***
## Residual          65     3838 0.06434                    
## Total             92    59652 1.00000                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# load data 
load("R_objects/stat_test_data.RData")

# Create subset
dat.clean <- dat
dat.num <- dat.clean %>% select(all_of(ends_with("acid")))

# Calculating Bray-Curtis PCoA with capscale
tmp2 <- capscale(aa.num ~ 1, comm = aa.num, distance = "euclidian", metaMDS = TRUE)

## Square root transformation
## Wisconsin double standardization

tmp2

## Call: capscale(formula = aa.num ~ 1, distance = "euclidian", comm =
## aa.num, metaMDSdist = TRUE)
## 
##               Inertia Rank
## Total           1.257     
## Unconstrained   1.257   18
## Inertia is squared Euclidean distance 
## Species scores projected from 'aa.num' 
## 
## Eigenvalues for unconstrained axes:
##   MDS1   MDS2   MDS3   MDS4   MDS5   MDS6   MDS7   MDS8 
## 0.6189 0.2642 0.1496 0.0648 0.0433 0.0377 0.0204 0.0165 
## (Showing 8 of 18 unconstrained eigenvalues)
## 
## metaMDSdist transformed data: wisconsin(sqrt(aa.num))

# Collect data for plotting
mds.samples <- data.frame(tmp2$CA$u)
mds.scfa <- data.frame(tmp2$CA$v)

# Prepare point zero and labels for arrows
mds.scfa$label1 <- row.names(mds.scfa)
mds.scfa$x <- 0
mds.scfa$y <- 0

# Bind with main data
dat.mds <- cbind(ano_df2, mds.samples)

plots <- list()
sets <- data.frame(Day = c(1,2,6))
for (i in 1:nrow(sets)) {
  tmp <- dat.mds[dat.mds$Day == sets$Day[i],]
  plots[[i]] <- ggplot() +
  geom_point(data = tmp, 
             mapping = aes(x = MDS1, y = MDS2, 
                           color = Strain, 
                           shape = Medium),
             size = 3, 
             alpha = 0.7) +
  theme_pubr(legend = "top") + 
  labs(x = "Axis 1", y = "Axis 2") +
  guides(fill = "none") +
  guides(size = "none") +
    scale_color_manual(values = c("A: Non-inoculated"="#E89242FF", 
                                  "B: MS22418 (PrtP-)"="#82491EFF", 
                                  "C: MS22421 (PrtP/SK11)"="#24325FFF", 
                                  "D: MS22425 (PrtP/Wg2)"="#B7E4F9FF", 
                                  "E: MS22427 (PrtP/MS22333)"="#FB6467FF"
                                  )
                       ) +
    coord_cartesian(xlim = c(min(mds.samples$MDS1)*1.1,
                              max(mds.samples$MDS1)*1.1),
                    ylim = c(min(mds.samples$MDS2)*1.1,
                             max(mds.samples$MDS2)*1.1))
  names(plots)[i] <- paste(sets$Medium[i], "day", sets$Day[i], sep = " ")
}
  plots[[4]] <- ggplot() +
  geom_segment(data = mds.scfa, mapping = aes(x=x, y=y, xend=0.35*MDS1, yend=0.35*MDS2),
               lineend = "butt",
               linejoin = "round",
               linewidth = 0.5,
               arrow = arrow(length = unit(0.3, 'cm'))) +
  geom_label_repel(data = mds.scfa,
             mapping = aes(x = 0.35*MDS1, y = 0.35*MDS2), #
             label = mds.scfa$label1,
             size = 4,
             min.segment.length = 0,
             segment.alpha = 0.8,
             box.padding = 0.3,
             force = 1) +
  theme_pubr(legend = "top") + 
  labs(x = "Axis 1", y = "Axis 2") +
  guides(fill = "none") +
    coord_cartesian(xlim = c(min(mds.samples$MDS1)*1.1, max(mds.samples$MDS1)*1.1), ylim = c(min(mds.samples$MDS2)*1.1,max(mds.samples$MDS2)*1.1))
  names(plots)[4] <- "AA weights"

p.out <- ggarrange(plotlist = plots, widths = 12, heights = 8,
          common.legend = TRUE, 
          legend = "bottom", 
          labels = names(plots)
          )