---
title: "Olsen_Analyses"
author: "JA"
date: "2024-01-16"
output: html_document
---


# Load Required Libraries
```{r}
library(readxl)   
library(calibrate) 
library(dplyr)    
library(ggplot2)  

# Load the configuration
config <- yaml::read_yaml("config.yml")

# Set base directory and subdirectories
base_directory <- config$base_directory
data_directory <- file.path(base_directory, config$data_directory)
results_directory <- file.path(base_directory, config$results_directory)


source(paste0(base_directory,"/common_functions/common_functions.R"))  # Custom functions
```


#Import Results from Excel Files
```{r}
# Load data from different sheets
HP <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                               sheet = "HILIC-pos", skip = 4), stringsAsFactors = FALSE)
HN <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                               sheet = "HILIC-neg", skip = 4), stringsAsFactors = FALSE)
CP <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                               sheet = "C8-pos", skip = 4), stringsAsFactors = FALSE)
CN <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                               sheet = "C18-neg", skip = 4), stringsAsFactors = FALSE)
```

#Import metadata
```{r}
# Load metadata from the same Excel file
HP_metadata <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                                        sheet = "HILIC-pos", range = "J1:V5"), stringsAsFactors = FALSE)
HN_metadata <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                                        sheet = "HILIC-neg", range = "J1:V5"), stringsAsFactors = FALSE)
CP_metadata <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                                        sheet = "C8-pos", range = "J1:V5"), stringsAsFactors = FALSE)
CN_metadata <- as.data.frame(read_excel(paste0(data_directory,"/22_0322_alphaSyn_fly_pilot_Classes.xlsx"),
                                        sheet = "C18-neg", range = "J1:V5"), stringsAsFactors = FALSE)

# Extract study metadata from HN sheet
s_m <- as.character(HN_metadata[2, -1])
```


#Combine and Clean Data
```{r}
# Combine all datasets into one dataframe
df <- rbind(HP, HN, CP, CN)

# Remove irrelevant columns by their names
df <- df[, !colnames(df) %in% c("...23", "0.05")]

# Rename specific columns for clarity
colnames(df)[colnames(df) == "p value...25"] <- "p value_Syn_vs_Ctrl"
colnames(df)[colnames(df) == "p value...28"] <- "p value_AKI_KD_Syn_vs_AKI_KD_Ctrl"
```


#Identify and Handle Duplicated Metabolites
Metabolites annotated in more than one 
```{r}
#Duplicated metabolites
dup_metabs <- df$Metabolite[!is.na(df$Metabolite)][duplicated(df$Metabolite[!is.na(df$Metabolite)])]
print(dup_metabs)

# Read decision file for duplicates and remove unwanted duplicates
toRemove <- read.csv(paste0(data_directory,"/dup_metabs_decision.csv"))
table(toRemove$Decision)
dim(df)
for(r in 1:nrow(toRemove)){
  if(toRemove$Decision[r] == ""){
    loc_toRm <- which(toRemove$Method[r] == df$Method & toRemove$Compound_ID[r] == df$Compound_ID)
    print(df[loc_toRm,c("Method","Metabolite")])
    df <- df[!(toRemove$Method[r] == df$Method & toRemove$Compound_ID[r] == df$Compound_ID),]
  }
}
dim(df)

#Check any duplicates left:
df$Metabolite[!is.na(df$Metabolite)][duplicated(df$Metabolite[!is.na(df$Metabolite)])]
```

#Remove Unwanted Data
```{r}
# Remove internal standards and incorrectly identified metabolites
df <- df[!grepl("nternal", df$HMDB_ID), ]
df <- df[!grepl("xypurinol", df$Metabolite), ]

# Remove rows with unknown metabolites
df <- df[!is.na(df$Metabolite), ]

```

#Split Data into Metadata and Values
```{r}
# Separate the metadata and measurement data
m <- df[, c(1:10, 23:26)]
a <- as.matrix(df[, c(11:22)])
```

#Impute values
```{r}
# Drop metabolites with no measurements
toDrop <- which(rowSums(is.na(a)) == ncol(a))
if(length(toDrop) > 0){
  m <- m[-toDrop, ]
  a <- a[-toDrop, ]
}

# Apply custom imputation methods
a <- rm.zeros(a)
a <- fun.impute(a)
a_backup <- a
```


#Impute Missing Values and Remove Zero-Variance Metabolites
```{r}
#Drop metabolites with no measurements
toDrop <- which(apply(a,1,FUN = function(x){return(sum(is.na(x)))}) == ncol(a))
#Empty metabolites
if(length(toDrop) > 0){
  m[toDrop,]
  a <- a[-toDrop,]
  m <- m[-toDrop,]
  
}

a <- rm.zeros(a)
a <- fun.impute(a)
a_backup <- a

# Remove zero-variance metabolites (often internal standards if IS-scaled data)
toKeep <- apply(a, 1, sd, na.rm = TRUE) != 0
a <- a[toKeep, ]
```


#Perform PCA (Principal Component Analysis)

```{r}
# Conduct PCA on the log-transformed data
pcaResults <- prcomp(t(log(a)), scale. = TRUE, center = TRUE)

# Create a dataframe with PCA results and additional metadata
PC <- as.data.frame(pcaResults$x) %>%
      mutate(names = rownames(pcaResults$x), s_m = s_m) %>%
      mutate(median = apply(a, 2, median, na.rm = TRUE),
             mean = apply(a, 2, mean, na.rm = TRUE),
             sum = apply(a, 2, sum, na.rm = TRUE))

# Plot PCA results
ggplot(PC, aes(x = PC1, y = PC2, color = s_m)) + 
  geom_point() + 
  theme(legend.title = element_blank()) + 
  ggtitle("PCA Scores Plot for Known Metabolites (HP, HN, CP, CN)") + 
  scale_x_continuous(name = paste0("[PCA1] ", var.fun(pcaResults)[1], " %")) +
  scale_y_continuous(name = paste0("[PCA2] ", var.fun(pcaResults)[2], " %"))

```



#Write 
```{r}
pdf(paste0(results_directory,"/pca_knowns_HP-HN-CP-CN.pdf"),8,8)
ggplot(PC, aes(x = PC1, y = PC2, color=s_m)) + 
  geom_point() + theme(legend.title = element_blank()) + ggtitle("PCA Scores plot knowns HP and HN, CP, CN") + 
  scale_x_continuous(name=paste0("[PCA1] ", var.fun(pcaResults)[1], " %")) +
  scale_y_continuous(name=paste0("[PCA2] ", var.fun(pcaResults)[2], " %"))
dev.off()
```

