example_code2.R

#import library
source("oac_pca_util.R")
library(acPCA)
library(ggplot2)
library(patchwork)
library(dplyr)
library(tidyr)
library(tibble)

file_path = "data\\exampledata2.xlsx"
location_names <- getSheetNames(file_path)

i <- 1
print(i)
name <- location_names[i]

palette <- c( "#F7D9E3","#DDF2FD","#D990B0", "#93BCD5","#BB487E", "#4987AD","#9E004C", "#005285")
split_num <- 48
xdata <- read_data(i ,file_path)$data
lipid_names <- read_data(i ,file_path)$lipidname
lipid_classes <- read_data(i ,file_path)$lipid_classes
xdata <- t(xdata)
xdata <- xdata[c(1:split_num),]


lipid_classes <- lipid_classes[colSums(xdata == 0) == 0]
lipid_names <- lipid_names[colSums(xdata == 0) == 0]
xdata <- xdata[, colSums(xdata == 0) == 0]


colnames(xdata) <- lipid_names

col_means <- colMeans(xdata, na.rm = TRUE)
col_lipid_names <- colnames(xdata)
keep_idx <- unlist(
  tapply(seq_along(col_lipid_names), col_lipid_names, function(i) {
    i[which.max(col_means[i])]
  })
)

xdata <- xdata[, keep_idx, drop = FALSE]
lipid_classes <- lipid_classes[keep_idx]
lipid_names <- lipid_names[keep_idx]

Ydata <- create_Y(xdata)

row_names <- rownames(xdata)
Env <- stringr::str_extract(row_names, "GF|SPF")
Sex <- stringr::str_extract(row_names, "_[MF]_") |> stringr::str_replace_all("_", "")
Age <- stringr::str_extract(row_names, "\\d+[wmM]")
Age_Env <- paste(Age, Env, sep = "_")
Age_Env_Sex <- paste(Age_Env,Sex,sep = '_')
weeks <- sapply(Age, to_weeks)
Age <- ifelse(
  grepl("w$", Age),
  sub("w$", " weeks", Age),
  sub("M$", " months", Age)
)
Age_Env <- Age_Env |>
  sub("^9w_", "2 months ", x = _) |>
  sub("^12M_", "12 months ", x = _) |>
  sub("^18M_", "19 months ", x = _) |>
  sub("^24M_", "24 months ", x = _)
X <- scale(xdata)
Y <- scale(Ydata)
pca_results <- prcomp(X, center=T)
pca_scores <- data.frame(
  PC1 = pca_results$x[, 1],
  PC2 = pca_results$x[, 2],
  Age_Env = Age_Env,
  Sex = Sex,
  Age_Env_Sex = Age_Env_Sex
)
pca_scores$Age_Env <- factor(
  pca_scores$Age_Env,
  levels = c("2 months GF","2 months SPF","12 months GF","12 months SPF","19 months GF","19 months SPF","24 months GF","24 months SPF"))
pca_contribution <- pca_results$sdev^2 / sum(pca_results$sdev^2) * 100
x_name = paste0("PC1 : ", round(pca_contribution[1], 1), "%")
y_name = paste0("PC2 : ", round(pca_contribution[2], 1), "%")
p_pca <- plot_pcascore_scatter(pca_scores,palette = palette,title = 'PCA',x_axis_name=x_name,y_axis_name=y_name)
p_pca_box <- plot_pcascore_box(pca_scores,x = Age_Env,y = PC2,fill = Sex,xlab = '')
p_pca
p_pca_box
###AC-PCA
h=8
result_tune <- acPCAtuneLambda(X=X, Y=Y, nPC=2, lambdas=seq(0, 10, 0.05),anov=F, kernel = "gaussian",
                               bandwidth=h, quiet=T)
ac_result <- acPCA(X=X, Y=Y, lambda=result_tune$best_lambda,kernel="gaussian", bandwidth=h, nPC=2)
acpca_scores <- data.frame(
  PC1 = ac_result$Xv[,1],
  PC2 = ac_result$Xv[,2],
  Age_Env = Age_Env,
  Sex = Sex,
  Age_Env_Sex = Age_Env_Sex
)
acpca_scores$Age_Env <- factor(
  acpca_scores$Age_Env,
  levels = c("2 months GF","2 months SPF","12 months GF","12 months SPF","19 months GF","19 months SPF","24 months GF","24 months SPF"))
p_acpca <- plot_pcascore_scatter(acpca_scores,palette = palette,title = 'AC-PCA')
###OAC-PCA
###t score
rho <- 0.5
x <- scale(xdata)
oac_results <- oacPCA(X,Y,rho)
oac_pca_scores_t <- data.frame(
  PC1 = oac_results$t[, 1],
  PC2 = oac_results$t[, 2],
  Age_Env = Age_Env,
  Sex = Sex,
  Age_Env_Sex = Age_Env_Sex
)
oac_pca_scores_t$Age_Env <- factor(
  oac_pca_scores_t$Age_Env,
  levels = c("2 months GF","2 months SPF","12 months GF","12 months SPF","19 months GF","19 months SPF","24 months GF","24 months SPF"))
p_oacpca_t <- plot_pcascore_scatter(oac_pca_scores_t,palette = palette,title = 'OAC-PCA (t)')
p_oacpca_t_box_PC1 <- plot_pcascore_box(oac_pca_scores_t,x = Age_Env,y = PC1,fill = Sex)
p_oacpca_t_box_PC2 <- plot_pcascore_box(oac_pca_scores_t,x = Age_Env,y = PC2,fill = Sex)

###s score
oac_pca_scores_s <- data.frame(
  PC1 = oac_results$s[, 1],
  PC2 = oac_results$s[, 2],
  Age_Env = Age_Env,
  Sex = Sex,
  Age_Env_Sex = Age_Env_Sex
)
oac_pca_scores_s$Age_Env <- factor(
  oac_pca_scores_s$Age_Env,
  levels = c("2 months GF","2 months SPF","12 months GF","12 months SPF","19 months GF","19 months SPF","24 months GF","24 months SPF"))
p_oacpca_s <- plot_pcascore_scatter(oac_pca_scores_s,palette = palette,title = 'OAC-PCA (s)')



##loading plot
new_names <- sub(".*\\|", "", lipid_names)
oacloading_1 <- oac_PCA_loading(X,oac_results$s,1)
oacloading_2 <- oac_PCA_loading(X,oac_results$s,2)
loading_df <- data.frame(
  value_1 = oacloading_1$loadings,
  q_val_1 = oacloading_1$q_values,
  value_2 =oacloading_2$loadings,
  q_val_2 = oacloading_2$q_values,
  index = new_names,
  lipid_names = lipid_names,
  lipid_classes = lipid_classes
)

top_bottom_df <- loading_df %>%
  arrange(desc(value_2)) %>%
  dplyr::slice(c(1:11, (n() - 9):n())) %>%
  mutate(
    rank_type = ifelse(row_number() <= 11, "Top", "Bottom"),
    label = index
  )

loading_plot2 <- plot_loading_top(df = top_bottom_df,value = value_2,title_name = 'PC2 loadings')
loading_plot2
top_bottom_df <- loading_df %>%
  arrange(desc(value_1)) %>%
  dplyr::slice(c(1:10, (n() - 9):n())) %>%
  mutate(
    rank_type = ifelse(row_number() <= 10, "Top", "Bottom"),
    label = index
  )

loading_plot1 <- plot_loading_top(df = top_bottom_df,value = value_1,title_name = 'PC1 loadings')

loading_plot1
library(grid)
library(gtable)

g <- ggplot2::ggplotGrob(loading_plot1)

panel_col <- g$layout[g$layout$name == "panel", "l"]

g$widths[panel_col] <- grid::unit(6, "cm")

grid::grid.newpage()
grid::grid.draw(g)


threshold <- 0.05
loading_df$color <- with(loading_df,
                         ifelse(q_val_1 <= threshold & q_val_2 <= threshold, "#E180BC",
                                ifelse(q_val_1 <= threshold, "#80DAE0",
                                       ifelse(q_val_2 <= threshold, "#E0D380", "gray"))))

label_index <- c(colnames(xdata))
loading_df$label_flag <- loading_df$index %in% label_index

loading_plot_pc12 <- plot_pc12(df = loading_df)