Clean up the environment

rm(list = ls())

Load necessary libraries for reading data, plotting, and statistical analysis

library(readxl)       # For reading Excel files
library(ggpubr)       # For creating publication-ready plots
library(dplyr)        # For data manipulation
library(tidyr)        # For reshaping data
library(ggh4x)        # For advanced facetting functions in ggplot2
library(rstatix)      # For statistical tests and adding significance markers
library(knitr)        # For displaying tibbles in a nice table format
library(kableExtra)   # For enhancing table formatting

Load data

# Load targeted serum, liver, brain, and urine data from the Excel file
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

# Define the list of metabolites for each tissue type
ls.serum.metabolites <- c("Tryptophan", "Kynurenine", "Quinolinate", "Serotonin")
ls.liver.metabolites <- c("Tryptophan", "Kynurenine", "NAD", "NAD+NADP","Serotonin")
ls.brain.metabolites <- c("Tryptophan", "Serotonin")
ls.urine.metabolites <- c("Tryptophan", "Quinolinate", "1-MN", "3-IS", "2PY")
Data.serum <- inner_join(Serum, PigInfo, by = "PigID") %>%
               select(Diet, all_of(ls.serum.metabolites)) %>%
               mutate(Region = "Serum") %>%
               pivot_longer(cols = -c(Diet, Region), names_to = "variable", values_to = "value") %>%
               mutate(Class = "I")


Data.liver <- inner_join(Liver, PigInfo, by = "PigID") %>%
               mutate(`NAD+NADP` = `NAD+` + `NADP+`) %>%
               rename(`NAD` = "NAD+") %>% # fix name
               select(Diet, all_of(ls.liver.metabolites)) %>%
               mutate(Region = "Liver") %>%
               pivot_longer(cols = -c(Diet, Region), names_to = "variable", values_to = "value") %>%
               mutate(Class = "II")

Data.brain <- inner_join(Brain, PigInfo, by = "PigID") %>%
               select(Diet, Region, all_of(ls.brain.metabolites)) %>%
               pivot_longer(cols = -c(Diet, Region), names_to = "variable", values_to = "value") %>%
               mutate(Class = "II")


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

Data.urine <- Urine %>%
                rename(`1-MN` = "1-Methylnicotinamide") %>% # fix name
                rename(`2PY` = "N-Methyl-2-pyridone-5-carboxamide") %>% # fix name
                rename(`3-IS` = "3-Indoxylsulfate") %>% # fix name
                mutate(across(`1-MN`:`β-Pseudouridine`, ~ .x / Effective.Osmolality)) %>% # Osmolality adjustment
                inner_join(PigInfo, by = "PigID") %>%
                select(Diet, all_of(ls.urine.metabolites)) %>%
                mutate(Region = "Urine") %>%
                pivot_longer(cols = -c(Diet, Region), names_to = "variable", values_to = "value") %>%
                mutate(Class = "III")
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")))
Fig <- All.data %>% 
          filter(variable %in% c("Kynurenine", "Quinolinate", "NAD+", "1-MN", "3-IS", "2PY", "NAD", "NAD+NADP")) %>%
          mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
          ggerrorplot(x = "variable", y = "value", color = "Diet", 
              palette = c("turquoise3", "purple", "orange"),
              na.rm = TRUE, 
              xlab = "", 
              ylab = "Concentration", 
              add = "mean_se", 
              size = 0.6,
              position = position_dodge(0.7)) +
              facet_nested(~Region + variable, scales = "free", independent = "y") +  # Free x axis with proportional facet widths using facet_nested
              theme_classic2() +
              theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

Fig

# save the figure in pdf:
ggsave(plot=Fig, height=4, width=13, dpi=300, filename="Figure 5/Fig5A.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") %>%
          select(Kynurenine, Quinolinate, Diet, SowID, Formula.consumed.g, Time.postprandial.min) %>%
          pivot_longer(cols = -c(Diet, SowID, Formula.consumed.g, Time.postprandial.min), names_to = "variable", values_to = "value") %>%
          filter(Diet %in% c("ALAC", "WPI")) %>%
          group_by(variable) %>%
          mutate(glog.value = glog(value*10)) %>%
          anova_test(glog.value ~ Diet + SowID + Formula.consumed.g + Time.postprandial.min, type = 1) %>%
          as_tibble() %>%
          filter(Effect == "Diet") %>%
          kable() %>% # Print the tibble in a nice table format
          kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
variable Effect DFn DFd F p p<.05 ges
Kynurenine Diet 1 17 2.946 0.104 0.148
Quinolinate Diet 1 17 6.247 0.023
0.269
Liver %>% mutate(`NAD+NADP` = `NAD+` + `NADP+`) %>%
          inner_join(PigInfo, by = "PigID") %>%
          inner_join(Nec_metadata, by = "PigID") %>%
          select(Kynurenine, `NAD+`, `NAD+NADP`, Diet, SowID, Formula.consumed.g, Time.postprandial.min) %>%
          pivot_longer(cols = -c(Diet, SowID, Formula.consumed.g, Time.postprandial.min), names_to = "variable", values_to = "value") %>%
          filter(Diet %in% c("ALAC", "WPI")) %>%
          group_by(variable) %>%
          mutate(glog.value = glog(value*10)) %>%
          anova_test(glog.value ~ Diet + SowID + Formula.consumed.g + Time.postprandial.min, type = 1) %>%
          as_tibble() %>%
          filter(Effect == "Diet") %>%
          kable() %>% # Print the tibble in a nice table format
          kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
variable Effect DFn DFd F p p<.05 ges
Kynurenine Diet 1 17 1.843 0.192 0.098
NAD+ Diet 1 17 5.673 0.029
0.250
NAD+NADP Diet 1 17 5.020 0.039
0.228

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

Serum %>% inner_join(PigInfo, by = "PigID") %>%
          inner_join(Nec_metadata, by = "PigID") %>%
          select(Kynurenine, Quinolinate, Diet, SowID) %>%
          pivot_longer(cols = -c(Diet, SowID), names_to = "variable", values_to = "value") %>%
          group_by(variable) %>%
          mutate(glog.value = glog(value*10)) %>%
          tukey_hsd(glog.value ~ Diet + SowID) %>%
          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") %>%
          kable() %>% # Print the tibble in a nice table format
          kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
variable term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
Kynurenine Diet SF WPI 0 -0.0320569 -0.0636105 -0.0005032 0.0460
Quinolinate Diet ALAC WPI 0 -1.0660824 -2.2532542 0.1210894 0.0842 #
Liver %>% mutate(`NAD+NADP` = `NAD+` + `NADP+`) %>%
          inner_join(PigInfo, by = "PigID") %>%
          inner_join(Nec_metadata, by = "PigID") %>%
          select(Kynurenine, `NAD+`, `NAD+NADP`, Diet, SowID) %>%
          pivot_longer(cols = -c(Diet, SowID), names_to = "variable", values_to = "value") %>%
          group_by(variable) %>%
          mutate(glog.value = glog(value*10)) %>%
          tukey_hsd(glog.value ~ Diet + SowID) %>%
          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") %>%
          kable() %>% # Print the tibble in a nice table format
          kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
variable term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
NAD+ Diet ALAC WPI 0 -0.2358007 -0.4688345 -0.0027670 0.0470
NAD+NADP Diet ALAC WPI 0 -0.1953640 -0.4032249 0.0124968 0.0682 #

Urine data, adjusting for sow ID

Urine %>% mutate(across(`1-Methylnicotinamide`:`β-Pseudouridine`, ~ .x / Effective.Osmolality)) %>% # osmolality adjustment
                inner_join(PigInfo, by = "PigID") %>%
                select(Quinolinate, `1-Methylnicotinamide`, `3-Indoxylsulfate`, `N-Methyl-2-pyridone-5-carboxamide`, Diet, SowID) %>%
                pivot_longer(cols = -c(Diet, SowID), names_to = "variable", values_to = "value") %>%
                group_by(variable) %>%
                mutate(glog.value = glog(value*10)) %>%
                tukey_hsd(glog.value ~ Diet + SowID) %>%
                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") %>%
                kable() %>% # Print the tibble in a nice table format
                kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting
variable term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
1-Methylnicotinamide Diet ALAC SF 0 -1.3392238 -2.153914 -0.5245339 0.00135
1-Methylnicotinamide Diet ALAC WPI 0 -0.9084302 -1.618660 -0.1982001 0.01100
3-Indoxylsulfate Diet ALAC SF 0 -1.1858386 -2.206152 -0.1655252 0.02110
3-Indoxylsulfate Diet ALAC WPI 0 -0.7875184 -1.677007 0.1019702 0.08860 #
N-Methyl-2-pyridone-5-carboxamide Diet ALAC SF 0 -1.1293781 -1.704696 -0.5540603 0.00021
N-Methyl-2-pyridone-5-carboxamide Diet ALAC WPI 0 -0.8046146 -1.306165 -0.3030642 0.00170
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] tidyr_1.3.0      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   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] car_3.1-2         pkgconfig_2.0.3   xml2_1.3.5        rmarkdown_2.28   
## [53] svglite_2.1.2     httr_1.4.7        rstudioapi_0.15.0 R6_2.5.1         
## [57] compiler_4.2.2