Clean up the environment
rm(list = ls())
Load libraries for data reading, plotting, 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 faceting functions in ggplot2
library(rstatix) # For statistical tests and adding significance markers
Load data
# Load hormone data
PP_hormones <- read_excel("Data/Data.xlsx", sheet = "PP_hormones")
Nec_hormones <- read_excel("Data/Data.xlsx", sheet = "Nec_hormones")
# Pig info
PigInfo <- read_excel("Data/Data.xlsx", sheet = "PigInfo")
PP.data <- PP_hormones %>% inner_join(PigInfo, by = "PigID") %>%
mutate(`Acylated ghrelin` = Ghrelin) %>%
mutate(`Desacylated ghrelin` = desacylGhrelin) %>%
mutate(GLP1 = GLP1.736 + GLP1.737) %>% # GLP-1 is the sum of GLP-1 736 and 737
mutate(Day = "Day 11 or 12") %>%
select(Diet, Day, TimePoint, `Acylated ghrelin`, `Desacylated ghrelin`, GLP1, Insulin)
Nec.data <- Nec_hormones %>% inner_join(PigInfo, by = "PigID") %>%
mutate(`Acylated ghrelin` = Ghrelin) %>%
mutate(`Desacylated ghrelin` = desacylGhrelin) %>%
mutate(GLP1 = GLP1.736 + GLP1.737) %>% # GLP-1 is the sum of GLP-1 736 and 737
mutate(Day = "Day 16") %>%
mutate(TimePoint = "") %>%
select(Diet, Day, TimePoint, `Acylated ghrelin`, `Desacylated ghrelin`, GLP1, Insulin)
# Combine postprandial and necropsy data
All.data <- bind_rows(PP.data, Nec.data) %>%
mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
mutate(TimePoint = factor(TimePoint, levels = c("baseline", "30 min", "60 min", "120 min")))
Plotting hormone levels by diet and time point
Fig1 <- All.data %>% filter(!is.na(`Acylated ghrelin`)) %>% # Remove NA values
ggerrorplot(x = "TimePoint", y = "Acylated ghrelin", color = "Diet",
palette = c("turquoise3","purple", "orange"),
na.rm = TRUE, xlab = "",
ylab ="Acylated ghrelin (ug/L)", add = "mean_se", position = position_dodge(0.6)) +
facet_nested(.~ Day , scales = "free", space = "free_x", independent = "y") +
theme_classic()
Fig2 <- All.data %>% filter(!is.na(`Desacylated ghrelin`)) %>% # Remove NA values
ggerrorplot(x = "TimePoint", y = "Desacylated ghrelin", color = "Diet",
palette = c("turquoise3","purple", "orange"),
na.rm = TRUE, xlab = "",
ylab ="Desacylated ghrelin (ug/L)", add = "mean_se", position = position_dodge(0.6)) +
facet_nested(.~ Day , scales = "free", space = "free_x", independent = "y") +
theme_classic()
Fig3 <- All.data %>% filter(!is.na(GLP1)) %>% # Remove NA values
ggerrorplot(x = "TimePoint", y = "GLP1", color = "Diet",
palette = c("turquoise3","purple", "orange"),
na.rm = TRUE, xlab = "",
ylab ="GLP-1 (ug/L)", add = "mean_se", position = position_dodge(0.6)) +
facet_nested(.~ Day , scales = "free", space = "free_x", independent = "y") +
theme_classic()
Fig4 <- All.data %>% filter(!is.na(Insulin)) %>% # Remove NA values
ggerrorplot(x = "TimePoint", y = "Insulin", color = "Diet",
palette = c("turquoise3","purple", "orange"),
na.rm = TRUE, xlab = "",
ylab ="Insulin (ug/L)", add = "mean_se", position = position_dodge(0.6)) +
facet_nested(.~ Day , scales = "free", space = "free_x", independent = "y") +
theme_classic()
# Combine all hormone plots into a single figure
Fig.all <- ggarrange(Fig1, Fig2, Fig3, Fig4, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")
Fig.all
# Save the figure in pdf format:
ggsave(plot=Fig.all, height=6.5, width=10, dpi=300, filename="SI Figure 3/SI.Fig3.pdf", useDingbats=FALSE)
Statistical analysis via Mann-Whitney U test
All.data %>% select(-Day) %>%
pivot_longer(cols = c(`Acylated ghrelin`, `Desacylated ghrelin`, GLP1, Insulin), names_to = "variable", values_to = "value") %>% # Reshape data to long format
group_by(variable, TimePoint) %>%
wilcox_test(value ~ Diet, paired = FALSE, p.adjust.method = "none") %>% # Perform Mann-Whitney U test (Wilcox test)
select(-p.adj) %>% # Remove adjusted p-values as they are not used
add_significance(cutpoints = c(0, 0.05, 0.1, 1), symbols = c("*", "#", "ns")) %>%
filter(p.signif != "ns") %>% # Filter out non-significant results
select(TimePoint, variable, group1, group2, p, p.signif)
## # A tibble: 5 × 6
## TimePoint variable group1 group2 p p.signif
## <fct> <chr> <chr> <chr> <dbl> <chr>
## 1 <NA> Acylated ghrelin ALAC WPI 0.069 #
## 2 baseline GLP1 ALAC WPI 0.075 #
## 3 120 min GLP1 ALAC WPI 0.088 #
## 4 <NA> GLP1 ALAC SF 0.062 #
## 5 60 min Insulin ALAC WPI 0.089 #
# Ghrelin pre and post- meal comparison
All.data %>% filter(Day == "Day 11 or 12") %>%
pivot_longer(cols = c(`Acylated ghrelin`, `Desacylated ghrelin`), names_to = "variable", values_to = "value") %>% # Reshape data to long format
group_by(variable) %>%
wilcox_test(value ~ TimePoint, paired = FALSE, p.adjust.method = "none") %>%
select(-p.adj) %>% # Remove adjusted p-values as they are not used
add_significance(cutpoints = c(0, 0.05, 0.1, 1), symbols = c("*", "#", "ns")) %>%
filter(p.signif != "ns") %>% # Filter out non-significant results
select(variable, group1, group2, p, p.signif)
## # A tibble: 7 × 5
## variable group1 group2 p p.signif
## <chr> <chr> <chr> <dbl> <chr>
## 1 Acylated ghrelin baseline 30 min 0.014 *
## 2 Acylated ghrelin baseline 120 min 0.093 #
## 3 Desacylated ghrelin baseline 30 min 0.0000299 *
## 4 Desacylated ghrelin baseline 60 min 0.0000025 *
## 5 Desacylated ghrelin baseline 120 min 0.00000553 *
## 6 Desacylated ghrelin 30 min 120 min 0.026 *
## 7 Desacylated ghrelin 60 min 120 min 0.033 *
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] rstatix_0.7.2 ggh4x_0.2.6 tidyr_1.3.0 dplyr_1.1.3 ggpubr_0.6.0
## [6] 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] htmltools_0.5.6.1 yaml_2.3.7 utf8_1.2.4 rlang_1.1.2
## [13] jquerylib_0.1.4 pillar_1.9.0 glue_1.6.2 withr_3.0.1
## [17] lifecycle_1.0.4 munsell_0.5.1 ggsignif_0.6.4 gtable_0.3.5
## [21] cellranger_1.1.0 ragg_1.2.6 evaluate_1.0.1 labeling_0.4.3
## [25] knitr_1.45 fastmap_1.1.1 fansi_1.0.6 highr_0.10
## [29] broom_1.0.5 scales_1.3.0 backports_1.4.1 cachem_1.0.8
## [33] jsonlite_1.8.8 abind_1.4-5 farver_2.1.1 systemfonts_1.0.5
## [37] textshaping_0.3.7 gridExtra_2.3 digest_0.6.33 grid_4.2.2
## [41] cowplot_1.1.1 cli_3.6.2 tools_4.2.2 magrittr_2.0.3
## [45] sass_0.4.7 tibble_3.2.1 car_3.1-2 pkgconfig_2.0.3
## [49] rmarkdown_2.28 rstudioapi_0.15.0 R6_2.5.1 compiler_4.2.2