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(reshape2)     # For reshaping data from wide to long format
library(ggh4x)        # For enhanced facet functions in ggplot2
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

# Load targeted metabolomics data from different tissues and biofluids
Serum <- read_excel("Data/Data.xlsx", sheet = "Nec_serum_targeted")
Liver <- read_excel("Data/Data.xlsx", sheet = "Liver_targeted")
Brain <- read_excel("Data/Data.xlsx", sheet = "Brain_targeted")
Urine <- read_excel("Data/Data.xlsx", sheet = "Nec_urine_NMR")

# 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

Data.serum <- inner_join(Serum, PigInfo, by = "PigID") %>%
               select(Diet, Tryptophan) %>%
               mutate(Region = "Serum") %>%
               melt(id = c("Diet", "Region")) %>%
               mutate(Class = "I")


Data.liver <- inner_join(Liver, PigInfo, by = "PigID") %>%
               select(Diet, Tryptophan) %>%
               mutate(Region = "Liver") %>%
               melt(id = c("Diet", "Region")) %>%
               mutate(Class = "II")

Data.brain <- inner_join(Brain, PigInfo, by = "PigID") %>%
               select(Diet, Region, Tryptophan) %>%
               melt(id = c("Diet", "Region")) %>%
               mutate(Class = "II")


#Calculate effective osmolality
Effective.Osmolality <- Urine$`Urine osmolality` - (Urine$Urea/1000)

Data.urine <- Urine %>%
                mutate(across(`1-Methylnicotinamide`:`β-Pseudouridine`, ~ .x / Effective.Osmolality)) %>% # osmolality adjustment
                inner_join(PigInfo, by = "PigID") %>%
                select(Diet, Tryptophan) %>%
                mutate(Region = "Urine") %>%
                melt(id = c("Diet", "Region")) %>%
                mutate(Class = "III")
# Combine all data into a single dataset for analysis
All.data <- bind_rows(Data.serum, Data.liver, Data.brain, Data.urine) %>%
            mutate(Region = factor(Region, levels = c("Serum", "Liver", "Hippocampus", "Hypothalamus", "PFC", "Striatum", "Urine")))
# Create an error plot for tryptophan concentration across different regions and diets
Fig <- All.data %>% filter(variable == "Tryptophan") %>%
                      mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                      ggerrorplot(x = "Region", y = "value", color = "Diet", 
                                 palette = c("turquoise3","purple", "orange"),
                                 na.rm = TRUE, xlab = "", 
                                 ylab ="Tryptophan concentration",  add = "mean_se", position = position_dodge(0.7)) +
                      facet_grid2(.~ Class , scales = "free", independent = "y", space = "free_x") + 
                      theme_classic() +
                      theme(strip.background = element_blank(), strip.text.x = element_blank())
                
Fig

# save the figure to a pdf file:
ggsave(plot=Fig, height=4, width=9, dpi=300, filename="Figure 3/Fig3B.pdf", useDingbats=FALSE)

Evaluate group differences

#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)
}

ALAC vs WPI: Serum and liver data, adjusting for intake volume and postprandial time

Serum %>% inner_join(PigInfo, by = "PigID") %>%
          inner_join(Nec_metadata, by = "PigID") %>%
            filter(Diet %in% c("ALAC", "WPI")) %>%
            mutate(glog.Tryptophan = glog(Tryptophan * 10)) %>% # Apply glog transformation
            anova_test(glog.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 0.048 0.830       0.003
## 2                 SowID   2  17 9.721 0.002     * 0.534
## 3    Formula.consumed.g   1  17 0.497 0.490       0.028
## 4 Time.postprandial.min   1  17 0.084 0.775       0.005
Liver %>% inner_join(PigInfo, by = "PigID") %>%
          inner_join(Nec_metadata, by = "PigID") %>%
            filter(Diet %in% c("ALAC", "WPI")) %>%
            mutate(glog.Tryptophan = glog(Tryptophan * 10)) %>% # Apply glog transformation
            anova_test(glog.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 3.526 0.078       0.172
## 2                 SowID   2  17 3.456 0.055       0.289
## 3    Formula.consumed.g   1  17 1.601 0.223       0.086
## 4 Time.postprandial.min   1  17 0.525 0.479       0.030

Between formula-fed groups and SF, adjusting for Sow ID

Serum %>% inner_join(PigInfo, by = "PigID") %>%
            mutate(glog.Tryptophan = glog(Tryptophan * 10)) %>% # Apply glog transformation
            tukey_hsd(glog.Tryptophan ~ Diet + SowID) %>% # Perform Tukey's HSD test
            filter(term == "Diet") %>%
            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.4425405 0.1949892 0.6900919 0.000460 ***
Diet ALAC WPI 0 0.0145481 -0.1890575 0.2181537 0.983000 ns
Diet SF WPI 0 -0.4279924 -0.6718761 -0.1841087 0.000566 ***
Liver %>% inner_join(PigInfo, by = "PigID") %>%
            mutate(glog.Tryptophan = glog(Tryptophan * 10)) %>% # Apply glog transformation
            tukey_hsd(glog.Tryptophan ~ Diet + SowID) %>%
            filter(term == "Diet") %>%
            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.2404133 -0.1660506 0.6468773 0.3190 ns
Diet ALAC WPI 0 -0.2386350 -0.5729428 0.0956728 0.1970 ns
Diet SF WPI 0 -0.4790483 -0.8794902 -0.0786064 0.0169

Liver, brain and urine data, adjusting for Sow ID

Brain %>% inner_join(PigInfo, by = "PigID") %>%
                 select(Diet, Region, Tryptophan, SowID) %>%
                 melt(id = c("Diet", "Region", "SowID")) %>%
                 mutate(glog.value = glog(value * 10)) %>% # Apply glog transformation
                 group_by(Region) %>%
                 tukey_hsd(glog.value ~ Diet + SowID) %>% # 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")) %>%
                 filter(p.adj.signif != "ns") %>%
                 select(Region, group1, group2, p.adj, p.adj.signif) %>%
                 kable() %>% # Print the tibble in a nice table format
                 kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
Region group1 group2 p.adj p.adj.signif
Hippocampus SF WPI 0.037600
Hypothalamus ALAC WPI 0.010400
Hypothalamus SF WPI 0.004670
PFC ALAC SF 0.049400
PFC ALAC WPI 0.016000
PFC SF WPI 0.000104
Striatum ALAC SF 0.001310
Striatum SF WPI 0.000159
Urine %>% mutate(across(`1-Methylnicotinamide`:`β-Pseudouridine`, ~ .x / Effective.Osmolality)) %>% 
          inner_join(PigInfo, by = "PigID") %>%
          select(Diet, Tryptophan, SowID) %>%
          melt(id = c("Diet","SowID")) %>%
          mutate(glog.value = glog(value * 10)) %>% # Apply glog transformation
          tukey_hsd(glog.value ~ Diet + SowID) %>% # 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")) %>%
          filter(p.adj.signif != "ns") %>%
          select(group1, group2, p.adj, p.adj.signif) %>%
          kable() %>% # Print the tibble in a nice table format
          kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
group1 group2 p.adj p.adj.signif
SF WPI 0.0426
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    ggh4x_0.2.6     
## [5] reshape2_1.4.4   dplyr_1.1.3      ggpubr_0.6.0     ggplot2_3.5.1   
## [9] 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   plyr_1.8.9        stringr_1.5.1    
## [21] munsell_0.5.1     ggsignif_0.6.4    gtable_0.3.5      cellranger_1.1.0 
## [25] ragg_1.2.6        rvest_1.0.3       evaluate_1.0.1    labeling_0.4.3   
## [29] fastmap_1.1.1     fansi_1.0.6       highr_0.10        broom_1.0.5      
## [33] Rcpp_1.0.11       scales_1.3.0      backports_1.4.1   cachem_1.0.8     
## [37] webshot_0.5.5     jsonlite_1.8.8    abind_1.4-5       farver_2.1.1     
## [41] systemfonts_1.0.5 textshaping_0.3.7 digest_0.6.33     stringi_1.7.12   
## [45] grid_4.2.2        cli_3.6.2         tools_4.2.2       magrittr_2.0.3   
## [49] sass_0.4.7        tibble_3.2.1      tidyr_1.3.0       car_3.1-2        
## [53] pkgconfig_2.0.3   xml2_1.3.5        svglite_2.1.2     httr_1.4.7       
## [57] rmarkdown_2.28    rstudioapi_0.15.0 R6_2.5.1          compiler_4.2.2