Clean up the environment

rm(list = ls())

Load required libraries for data reading, plotting, and data manipulation

library(readxl)       # For reading Excel files
library(ggpubr)       # For creating publication-quality plots
library(dplyr)        # For data manipulation
library(rstatix)      # For statistical analysis and adding significance markers
library(knitr)        # For displaying tibbles in a nice table format
library(kableExtra)   # For enhancing table formatting

Load data

# Serum metabolome data at necropsy
Serum <- read_excel("Data/Data.xlsx", sheet = "Nec_serum_targeted")

# Pig info
PigInfo <- read_excel("Data/Data.xlsx", sheet = "PigInfo")
# Metadata related to necropsy, including feed intake information
Nec_metadata <- read_excel("Data/Data.xlsx", sheet = "Nec_metadata")

Data organization and figure generation

# Join serum data with pig information, reorder diet levels, and create an error plot for serum free tryptophan
Fig1 <- Serum %>% inner_join(PigInfo, by = "PigID") %>%
          mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
          ggerrorplot(y = "Free.Tryptophan", x = "Diet", color = "Diet", 
            palette = c("turquoise3", "purple", "orange"),
            na.rm = TRUE, 
            ylab = "Serum free tryptophan (uM)", 
            xlab = "", 
            add = "mean_se", 
            size = 0.6) +
            theme_classic2() +
          theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

# Save the figure in pdf format
ggsave(plot=Fig1, height=4, width=2.6, dpi=300, filename="Figure 3/Fig3A_left.pdf", useDingbats=FALSE)

Fig1

Evaluate group difference

#This function transforms the input values by the generalized #log function.
glog <- function(y) {
  #Using lambda = 1 
  yt <- log(y+sqrt(y^2+1))
  return(yt)
}

Free Trp - between ALAC and WPI (adjust for Sow ID, formula volume and postprandial time)

Serum %>% inner_join(PigInfo, by = "PigID") %>%
          inner_join(Nec_metadata, by = "PigID") %>%
            filter(Diet %in% c("ALAC", "WPI")) %>%
            mutate(glog.Free.Tryptophan = glog(Free.Tryptophan * 10)) %>% # Apply generalized log transformation
            anova_test(glog.Free.Tryptophan ~ Diet + SowID + Formula.consumed.g + Time.postprandial.min, type = 1)
## ANOVA Table (type I tests)
## 
##                  Effect DFn DFd     F     p p<.05      ges
## 1                  Diet   1  17 4.459 0.050     * 0.208000
## 2                 SowID   2  17 2.053 0.159       0.195000
## 3    Formula.consumed.g   1  17 0.005 0.944       0.000304
## 4 Time.postprandial.min   1  17 0.009 0.927       0.000513

Free Trp - between formula-fed groups and SF (adjust for Sow ID)

Serum %>% inner_join(PigInfo, by = "PigID") %>%
            mutate(glog.Free.Tryptophan = glog(Free.Tryptophan*10)) %>%
            tukey_hsd(glog.Free.Tryptophan ~ Diet + SowID, type = 1) %>% # Perform Tukey's HSD test
            filter(term == "Diet") %>%
            add_significance(p.col = "p.adj", cutpoints = c(0, 0.05, 0.1, 1), symbols = c("*", "#", "ns")) %>%
            kable() %>% # Print the tibble in a nice table format
            kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
Diet ALAC SF 0 0.2430263 -1.308361 1.7944136 0.9190 ns
Diet ALAC WPI 0 -1.2273333 -2.503316 0.0486492 0.0610 #
Diet SF WPI 0 -1.4703595 -2.998762 0.0580430 0.0609 #

Estiamte % of free trp

Serum %>% inner_join(PigInfo, by = "PigID") %>%
            mutate(Per_free_trp = Free.Tryptophan / Tryptophan * 100) %>%
            get_summary_stats(Per_free_trp, type = "common") %>% # Get summary statistics for percentage of free tryptophan
            kable() # Print the tibble in a nice table format
variable n min max median iqr mean sd se ci
Per_free_trp 29 0 16.025 8.489 4.209 8.121 3.796 0.705 1.444 
Fig2 <- Serum %>% inner_join(PigInfo, by = "PigID") %>%
          mutate(Per_free_trp = Free.Tryptophan / Tryptophan * 100) %>%
          mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
          ggerrorplot(y = "Per_free_trp", x = "Diet", color = "Diet", 
            palette = c("turquoise3", "purple", "orange"),
            na.rm = TRUE, 
            ylab = "% of free tryptophan", 
            xlab = "", 
            add = "mean_se", 
            size = 0.6) +
            theme_classic2() +
            theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

Fig2

# save the figure in pdf format:
ggsave(plot=Fig2, height=4, width=2.6, dpi=300, filename="Figure 3/Fig3A_right.pdf", useDingbats=FALSE)

% free Trp - between ALAC and WPI (adjust for Sow ID, formula volume and postprandial time)

Serum %>% inner_join(PigInfo, by = "PigID") %>%
          inner_join(Nec_metadata, by = "PigID") %>%
          mutate(Per_free_trp = Free.Tryptophan / Tryptophan * 100) %>%
          filter(Diet %in% c("ALAC", "WPI")) %>%
          anova_test(Per_free_trp ~ Diet + SowID + Formula.consumed.g + Time.postprandial.min, type = 1)
## ANOVA Table (type I tests)
## 
##                  Effect DFn DFd      F     p p<.05   ges
## 1                  Diet   1  17 12.014 0.003     * 0.414
## 2                 SowID   2  17  1.179 0.332       0.122
## 3    Formula.consumed.g   1  17  0.048 0.830       0.003
## 4 Time.postprandial.min   1  17  0.079 0.782       0.005

% free Trp - between formula-fed groups and SF (adjust for Sow ID)

Serum %>% inner_join(PigInfo, by = "PigID") %>%
            mutate(Per_free_trp = Free.Tryptophan / Tryptophan * 100) %>%
            tukey_hsd(Per_free_trp ~ Diet + SowID, type = 1) %>%
            filter(term == "Diet") %>%
            add_significance(p.col = "p.adj", cutpoints = c(0, 0.05, 0.1, 1), symbols = c("*", "#", "ns")) %>%
            kable() %>% # Print the tibble in a nice table format
            kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
Diet ALAC SF 0 -2.077538 -6.085625 1.930550 0.41200 ns
Diet ALAC WPI 0 -4.987122 -8.283687 -1.690557 0.00256
Diet SF WPI 0 -2.909585 -6.858290 1.039120 0.17800 ns 
sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] kableExtra_1.3.4 knitr_1.45       rstatix_0.7.2    dplyr_1.1.3     
## [5] ggpubr_0.6.0     ggplot2_3.5.1    readxl_1.4.3    
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0  xfun_0.40         bslib_0.5.1       purrr_1.0.2      
##  [5] carData_3.0-5     colorspace_2.1-0  vctrs_0.6.5       generics_0.1.3   
##  [9] viridisLite_0.4.2 htmltools_0.5.6.1 yaml_2.3.7        utf8_1.2.4       
## [13] rlang_1.1.2       jquerylib_0.1.4   pillar_1.9.0      glue_1.6.2       
## [17] withr_3.0.1       lifecycle_1.0.4   stringr_1.5.1     munsell_0.5.1    
## [21] ggsignif_0.6.4    gtable_0.3.5      cellranger_1.1.0  ragg_1.2.6       
## [25] rvest_1.0.3       evaluate_1.0.1    labeling_0.4.3    fastmap_1.1.1    
## [29] fansi_1.0.6       highr_0.10        broom_1.0.5       scales_1.3.0     
## [33] backports_1.4.1   cachem_1.0.8      webshot_0.5.5     jsonlite_1.8.8   
## [37] abind_1.4-5       farver_2.1.1      systemfonts_1.0.5 textshaping_0.3.7
## [41] digest_0.6.33     stringi_1.7.12    grid_4.2.2        cli_3.6.2        
## [45] tools_4.2.2       magrittr_2.0.3    sass_0.4.7        tibble_3.2.1     
## [49] tidyr_1.3.0       car_3.1-2         pkgconfig_2.0.3   xml2_1.3.5       
## [53] rmarkdown_2.28    svglite_2.1.2     httr_1.4.7        rstudioapi_0.15.0
## [57] R6_2.5.1          compiler_4.2.2