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