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
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