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