Clean up the environment

rm(list = ls())

Load library

library(readxl)       # For reading Excel files
library(ggpubr)       # For creating publication-ready plots
library(dplyr)        # For data manipulation

Load data

# Load serum, urine, liver and brain metabolome data
Serum <- read_excel("Data/Data.xlsx", sheet = "Nec_serum_targeted")
Urine <- read_excel("Data/Data.xlsx", sheet = "Nec_urine_NMR")
Liver <- read_excel("Data/Data.xlsx", sheet = "Liver_targeted")
Brain <- read_excel("Data/Data.xlsx", sheet = "Brain_targeted")

# Load pig information
PigInfo <- read_excel("Data/Data.xlsx", sheet = "PigInfo")

Data organization and transformation

#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.data <- Serum %>% inner_join(PigInfo, by = "PigID") %>%
                  mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                  mutate(across(`2-Aminoadipate` : `Kynurenine`, ~.x * 10 )) %>%
                  mutate(across(`2-Aminoadipate` : `Kynurenine`, .fns = list(log = glog), .names = "{.col}_{.fn}")) %>%
                  na.omit()

Serum.data.log <- Serum.data %>% select(`2-Aminoadipate_log` : `Kynurenine_log`)
#Calculate effective osmolality
Effective.Osmolality <- Urine$`Urine osmolality` - (Urine$Urea/1000)

Urine.data <- Urine %>% inner_join(PigInfo, by = "PigID") %>%
                  mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                  mutate(across(`1-Methylnicotinamide` : `β-Pseudouridine`, ~.x / Effective.Osmolality * 100 )) %>% # adjust by effective osmolality
                  mutate(across(`1-Methylnicotinamide` : `β-Pseudouridine`, .fns = list(log = glog), .names = "{.col}_{.fn}")) %>%
                  na.omit()

Urine.data.log <- Urine.data %>% select(`1-Methylnicotinamide_log` : `β-Pseudouridine_log`)
Liver.data <- Liver %>% inner_join(PigInfo, by = "PigID") %>%
                          mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                          mutate(across(`2-Aminobutyrate` : `Kynurenine`, ~.x * 10 )) %>%
                          mutate(across(`2-Aminobutyrate` : `Kynurenine`, .fns = list(log = glog), .names = "{.col}_{.fn}")) %>%
                          na.omit()

Liver.data.log <- Liver.data %>% select(`2-Aminobutyrate_log` : `Kynurenine_log`)
Hippocampus.data <- Brain %>% inner_join(PigInfo, by = "PigID") %>%
                                mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                                filter(Region == "Hippocampus") %>%
                                mutate(across(`2-Aminobutyrate` : `Kynurenine`, ~.x * 10 )) %>%
                                mutate(across(`2-Aminobutyrate` : `Kynurenine`, .fns = list(log = glog), .names = "{.col}_{.fn}")) %>%
                                na.omit()

Hippocampus.data.log <- Hippocampus.data %>% select(`2-Aminobutyrate_log` : `Kynurenine_log`)
Hypothalamus.data <- Brain %>% inner_join(PigInfo, by = "PigID") %>%
                       mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                       filter(Region == "Hypothalamus") %>%
                       mutate(across(`2-Aminobutyrate` : `Kynurenine`, ~.x * 10 )) %>%
                       mutate(across(`2-Aminobutyrate` : `Kynurenine`, .fns = list(log = glog), .names = "{.col}_{.fn}")) %>%
                       na.omit()

Hypothalamus.data.log <- Hypothalamus.data %>% select(`2-Aminobutyrate_log` : `Kynurenine_log`)
PFC.data <- Brain %>% inner_join(PigInfo, by = "PigID") %>%
                        mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                        filter(Region == "PFC") %>%
                        mutate(across(`2-Aminobutyrate` : `Kynurenine`, ~.x * 10 )) %>%
                        mutate(across(`2-Aminobutyrate` : `Kynurenine`, .fns = list(log = glog), .names = "{.col}_{.fn}")) %>%
                        na.omit()

PFC.data.log <- PFC.data %>% select(`2-Aminobutyrate_log` : `Kynurenine_log`)
Striatum.data <- Brain %>% inner_join(PigInfo, by = "PigID") %>%
                              mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                              filter(Region == "Striatum") %>%
                              mutate(across(`2-Aminobutyrate` : `Kynurenine`, ~.x * 10 )) %>%
                              mutate(across(`2-Aminobutyrate` : `Kynurenine`, .fns = list(log = glog), .names = "{.col}_{.fn}")) %>%
                              na.omit()

Striatum.data.log <- Striatum.data %>% select(`2-Aminobutyrate_log` : `Kynurenine_log`)

PCA

Plot.PCA.with.centroid<-function(metabolites.log, meta.df, title){
  # metabolites.log = metabolites after log transformation
  # meta.df = dietary group
  pca = mixOmics::pca(metabolites.log, ncomp = 4, scale = F, center = T)
  df = data.frame(pc1 = pca$variates$X[,1], 
                  pc2 = pca$variates$X[,2], 
                  Diet = meta.df)


  centroids.x = aggregate(pc1~Diet, df, mean)
  centroids.y = aggregate(pc2~Diet, df, mean)

  centroids = inner_join(centroids.x, centroids.y, by = "Diet") %>%
              rename(pc1.mean = pc1) %>%
              rename(pc2.mean = pc2)

  df = df %>% inner_join(centroids, by = "Diet")


p.centroid = ggplot(df, aes(x = pc1, y = pc2, colour = Diet)) +
  geom_point(size=2.5) + 
  xlab(paste("PC1 (", round(unlist(pca$prop_expl_var)[1]*100,1), "%)", sep="")) +
  ylab(paste("PC2 (", round(unlist(pca$prop_expl_var)[2]*100,1), "%)", sep="")) +

  geom_point(aes(x=pc1.mean, y=pc2.mean),size=4, shape = 1) +
  geom_segment(aes(x=pc1.mean, y=pc2.mean, xend=pc1, yend=pc2), alpha = 0.8) +

  scale_color_manual(values =  c("turquoise3","purple", "orange")) + 
  ggtitle(title) + 
  stat_ellipse(type = "norm", linetype = 2) +
  
  theme_bw() +

  theme(axis.line = element_line(colour = "black"),
        plot.title = element_text(hjust = 0.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title = element_text(size = 12),
        axis.text.x = element_text(size = 11), 
        axis.text.y = element_text(size = 11), 
        panel.background = element_rect(colour = "black", linewidth=0.5))

return(p.centroid)
}
Fig1 <- Plot.PCA.with.centroid(Serum.data.log, Serum.data$Diet, title = "Serum")
Fig2 <- Plot.PCA.with.centroid(Urine.data.log, Urine.data$Diet, title = "Urine")
Fig3 <- Plot.PCA.with.centroid(Liver.data.log, Liver.data$Diet, title = "Liver")
Fig4 <- Plot.PCA.with.centroid(Hippocampus.data.log, Hippocampus.data$Diet, title = "Hippocampus")
Fig5 <- Plot.PCA.with.centroid(Hypothalamus.data.log, Hypothalamus.data$Diet, title = "Hypothalamus")
Fig6 <- Plot.PCA.with.centroid(PFC.data.log, PFC.data$Diet, title = "Prefrontal cortex")
Fig7 <- Plot.PCA.with.centroid(Striatum.data.log, Striatum.data$Diet, title = "Striatum")
Fig.top <- ggarrange(Fig1, Fig2, Fig3, ncol = 3, nrow = 1, legend = "none")
Fig.bottom <- ggarrange(Fig4, Fig5, Fig6, Fig7, ncol = 4, nrow = 1, legend = "bottom", common.legend = TRUE)
Fig.all <- ggarrange(Fig.top, Fig.bottom, ncol = 1, nrow = 2, legend = "bottom", common.legend = TRUE)

Fig.all

# Save the figure in pdf format:
ggsave(plot=Fig.all, height=8, width=14, dpi=300, filename="SI Figure 6/SI.Fig6.pdf", useDingbats=FALSE)
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] dplyr_1.1.3   ggpubr_0.6.0  ggplot2_3.5.1 readxl_1.4.3 
## 
## loaded via a namespace (and not attached):
##  [1] ggrepel_0.9.4       Rcpp_1.0.11         lattice_0.20-45    
##  [4] tidyr_1.3.0         corpcor_1.6.10      digest_0.6.33      
##  [7] utf8_1.2.4          RSpectra_0.16-1     plyr_1.8.9         
## [10] R6_2.5.1            cellranger_1.1.0    backports_1.4.1    
## [13] ellipse_0.5.0       evaluate_1.0.1      highr_0.10         
## [16] pillar_1.9.0        rlang_1.1.2         rstudioapi_0.15.0  
## [19] car_3.1-2           jquerylib_0.1.4     Matrix_1.5-4.1     
## [22] rmarkdown_2.28      textshaping_0.3.7   labeling_0.4.3     
## [25] rARPACK_0.11-0      BiocParallel_1.32.6 stringr_1.5.1      
## [28] igraph_1.5.1        munsell_0.5.1       broom_1.0.5        
## [31] compiler_4.2.2      xfun_0.40           systemfonts_1.0.5  
## [34] pkgconfig_2.0.3     htmltools_0.5.6.1   tidyselect_1.2.0   
## [37] tibble_3.2.1        gridExtra_2.3       codetools_0.2-19   
## [40] matrixStats_1.0.0   fansi_1.0.6         withr_3.0.1        
## [43] MASS_7.3-60         grid_4.2.2          jsonlite_1.8.8     
## [46] gtable_0.3.5        lifecycle_1.0.4     magrittr_2.0.3     
## [49] scales_1.3.0        stringi_1.7.12      cli_3.6.2          
## [52] cachem_1.0.8        carData_3.0-5       farver_2.1.1       
## [55] ggsignif_0.6.4      reshape2_1.4.4      bslib_0.5.1        
## [58] ragg_1.2.6          generics_0.1.3      vctrs_0.6.5        
## [61] cowplot_1.1.1       RColorBrewer_1.1-3  mixOmics_6.22.0    
## [64] tools_4.2.2         glue_1.6.2          purrr_1.0.2        
## [67] abind_1.4-5         parallel_4.2.2      fastmap_1.1.1      
## [70] yaml_2.3.7          colorspace_2.1-0    rstatix_0.7.2      
## [73] knitr_1.45          sass_0.4.7