rm(list = ls())
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
# 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")
# 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
#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)
}
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
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 | # |
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)
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
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