# Load necessary libraries
library(tidyverse) # For data manipulation and visualization
library(RColorBrewer) # For color palettes

# --- Simulate a biological dataset: Protein expression profiles ---
# Imagine 50 protein samples, each with 10 features (e.g., phosphorylation sites, PTMs, abundance levels)
# Introduce some biological groups for demonstration (e.g., 'TreatmentA', 'TreatmentB', 'Control')
set.seed(123) # For reproducibility

num_samples <- 50
num_features <- 10

# Create feature names
feature_names <- paste0("Protein_Feature_", 1:num_features)

# Generate random data
data_matrix <- matrix(rnorm(num_samples * num_features, mean = 100, sd = 15), 
                      nrow = num_samples, ncol = num_features)
colnames(data_matrix) <- feature_names

# Introduce group-specific variations to make PCA interesting
groups <- sample(c("Control", "TreatmentA", "TreatmentB"), num_samples, replace = TRUE, prob = c(0.4, 0.3, 0.3))

# Apply group-specific shifts (e.g., TreatmentA increases certain features, TreatmentB decreases others)
for (i in 1:num_samples) {
  if (groups[i] == "TreatmentA") {
    data_matrix[i, c(1, 3, 5)] <- data_matrix[i, c(1, 3, 5)] + rnorm(3, 20, 5) # Increase specific features
    data_matrix[i, c(2, 4)] <- data_matrix[i, c(2, 4)] + rnorm(2, 5, 2) # Slightly increase others
  } else if (groups[i] == "TreatmentB") {
    data_matrix[i, c(6, 8, 10)] <- data_matrix[i, c(6, 8, 10)] - rnorm(3, 15, 5) # Decrease specific features
    data_matrix[i, c(7, 9)] <- data_matrix[i, c(7, 9)] - rnorm(2, 3, 1) # Slightly decrease others
  }
}

# Convert to a data frame for easier manipulation
protein_data <- as.data.frame(data_matrix)
protein_data$Group <- groups

# Display the first few rows and structure
head(protein_data)
str(protein_data)

# --- Pre-processing steps: Scaling is Crucial for PCA ---
# Isolate numerical features for PCA
features_for_pca <- protein_data %>% 
  select(-Group)

# Activate scaling: Center and scale the data
# PCA is sensitive to the scale of variables. Scaling ensures that features
# with larger variances do not disproportionately influence the principal components.
# `scale()` function centers (mean=0) and scales (sd=1) the data by default.
scaled_data <- scale(features_for_pca)

# Verify scaling (mean should be close to 0, sd close to 1)
# apply(scaled_data, 2, mean) # Should be ~0
# apply(scaled_data, 2, sd)   # Should be ~1

# Convert scaled data back to a data frame if needed for certain packages
scaled_df <- as.data.frame(scaled_data)
head(scaled_df)

# --- Perform PCA using prcomp() --- 
# The `prcomp()` function is recommended over `princomp()` as it uses singular value decomposition (SVD),
# which is generally more robust for numerical stability.
# `scale = TRUE` is often redundant if you've already scaled, but it's good practice to ensure.
# We will use the 'scaled_data' from the previous step.
protein_pca_results <- prcomp(scaled_data, scale. = FALSE) # Set scale. = FALSE as data is already scaled

# Display the basic PCA output
print(protein_pca_results)

# --- Interpret Eigenvalues and Explained Variance ---
# Eigenvalues represent the variance explained by each principal component.
# They are derived from the singular values (sdev) squared.
# `sdev` are the standard deviations of the principal components.

eigenvalues <- (protein_pca_results$sdev)^2

# Calculate the proportion of variance explained by each component
variance_explained <- eigenvalues / sum(eigenvalues)

# Calculate the cumulative proportion of variance explained
cumulative_variance <- cumsum(variance_explained)

# Combine into a data frame for easier viewing
variance_df <- data.frame(
  PC = 1:length(eigenvalues),
  Eigenvalue = eigenvalues,
  Proportion_of_Variance = variance_explained,
  Cumulative_Variance = cumulative_variance
)

head(variance_df)

# --- Visualize the Scree Plot --- 
# A scree plot helps determine the number of principal components to retain.
# It plots the eigenvalues in decreasing order.

# Using ggplot2 for a cleaner scree plot
scree_plot <- ggplot(variance_df, aes(x = PC, y = Proportion_of_Variance)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_line(aes(y = Cumulative_Variance), color = "red", group = 1) +
  geom_point(aes(y = Cumulative_Variance), color = "red") +
  labs(title = "Scree Plot: Variance Explained by Principal Components",
       x = "Principal Component",
       y = "Proportion of Variance Explained") +
  theme_minimal()
print(scree_plot)

# --- Interpret Loadings (Variable Contributions) ---
# Loadings represent the correlation between original variables and principal components.
# High absolute loading values indicate that the original variable contributes strongly to that PC.
# These are stored in `protein_pca_results$rotation`.

loadings_matrix <- protein_pca_results$rotation

# View loadings for the first few principal components
# Sort by absolute value to easily identify the strongest contributors

# PC1 loadings
sorted_loadings_PC1 <- sort(abs(loadings_matrix[, "PC1"]), decreasing = TRUE)
print("Top contributors to PC1:")
print(loadings_matrix[names(sorted_loadings_PC1)[1:5], "PC1"])

# PC2 loadings
sorted_loadings_PC2 <- sort(abs(loadings_matrix[, "PC2"]), decreasing = TRUE)
print("Top contributors to PC2:")
print(loadings_matrix[names(sorted_loadings_PC2)[1:5], "PC2"])

# Often useful to visualize loadings using a biplot or separate loading plot
# (Will be covered in the next section for better visualization)

# Load `factoextra` for enhanced PCA visualization
# install.packages("factoextra") # Uncomment and run if not installed
library(factoextra)
library(ggplot2)

# --- Prepare data for visualization ---
# Add the group information back to the PCA results for coloring
# `fviz_pca_ind` plots the samples (individuals)
# `fviz_pca_var` plots the variables (features)

# Combine sample group information with PCA scores
# Scores are the coordinates of the samples in the new PC space
protein_pca_scores <- as.data.frame(protein_pca_results$x)
protein_pca_scores$Group <- protein_data$Group

# --- Visualize Score Plot: Sample Distribution ---
# This plot shows how samples cluster based on their principal component scores.
# We will color samples by their biological groups to reveal potential separations.
score_plot <- fviz_pca_ind(protein_pca_results, 
                           col.ind = "Group", # Color by group
                           palette = brewer.pal(n = 3, name = "Set1"), # Use a color palette
                           addEllipses = TRUE, # Add confidence ellipses around groups
                           ellipse.type = "convex", # Type of ellipse
                           repel = TRUE, # Avoid text overlapping
                           legend.title = "Biological Group",
                           title = "PCA Score Plot: Protein Samples by Group") +
  theme_minimal()
print(score_plot)

# --- Visualize Biplot: Samples and Variables Combined ---
# A biplot simultaneously displays both sample scores and variable loadings.
# It allows us to see which features contribute to the separation of samples.

# fviz_pca_biplot combines individuals (samples) and variables (features)
biplot <- fviz_pca_biplot(protein_pca_results, 
                          col.ind = protein_data$Group, # Color individuals by group
                          palette = brewer.pal(n = 3, name = "Set1"), 
                          addEllipses = TRUE, 
                          ellipse.type = "convex",
                          col.var = "black", # Color variables (arrows) black
                          repel = TRUE, 
                          legend.title = "Biological Group",
                          title = "PCA Biplot: Samples and Protein Features") +
  theme_minimal()
print(biplot)

# --- Advanced Biplot Interpretation Tip ---
# If the biplot becomes too cluttered, you can selectively display only the top contributing variables.
# Example: Only show top 10 variables by their contribution to PC1+PC2
# Get the contribution of variables to PC1 and PC2
# var_contrib <- get_pca_var(protein_pca_results)$contrib
# top_vars <- names(sort(rowSums(var_contrib[,1:2]), decreasing = TRUE)[1:5]) # Top 5 for clarity

# This is a bit more complex with fviz_pca_biplot, often easier to manually plot or adjust parameters.
# For simpler biplots, base R's `biplot()` function works, but `factoextra` offers more customization.
# For demonstration, we rely on the `repel = TRUE` for clarity, but be mindful of variable overlap.

# --- Best Practices: Validating PCA Results ---
# 1. Reproducibility: Ensure your scripts are reproducible with `set.seed()`.
# 2. Robustness Check: Consider running PCA on subsets of data or using robust PCA methods
#    (e.g., from `rrcov` package) for datasets with outliers.

# --- Advanced Visualization with factoextra for specific contributions ---
# Often, we want to visualize the contributions of variables specifically to each PC.

# Plot the contributions of variables to PC1
fviz_contrib(protein_pca_results, choice = "var", axes = 1, top = 10, 
             title = "Top 10 Protein Feature Contributions to PC1",
             fill = "steelblue", color = "steelblue") +
  theme_minimal()

# Plot the contributions of variables to PC2
fviz_contrib(protein_pca_results, choice = "var", axes = 2, top = 10,
             title = "Top 10 Protein Feature Contributions to PC2",
             fill = "darkred", color = "darkred") +
  theme_minimal()

# --- Interpreting Clusters with Metadata ---
# Once clusters are identified in the score plot, cross-reference them with your original metadata.
# For example, if a cluster emerges, confirm if it corresponds to a specific biological condition,
# time point, or genetic background.

# Example: After observing clusters in the score plot,
# we can examine which features primarily define these groups.
# This involves looking at the loadings of PCs that separate the clusters.

# --- Handling Outliers ---
# Outliers can heavily influence PCA results. Identify them through score plots (samples far from others)
# and consider their biological context. Sometimes, outliers represent novel biology; other times, they are errors.
# Robust PCA methods can mitigate the effect of extreme outliers.

# --- Integrating PCA with Downstream Analysis ---
# The principal component scores themselves can be used as input for other statistical models,
# such as classification (e.g., SVM, Random Forest) or regression.
# This is a powerful way to reduce dimensionality while retaining most variance.

# Example: Using PC scores for a simple clustering algorithm (e.g., K-means)
# (Not executable without selecting a K, but demonstrates the concept)
# kmeans_result <- kmeans(protein_pca_scores[, c("PC1", "PC2")], centers = 3)
# protein_pca_scores$kmeans_cluster <- as.factor(kmeans_result$cluster)

# fviz_pca_ind(protein_pca_results, 
#              col.ind = protein_pca_scores$kmeans_cluster, 
#              palette = "jco", 
#              addEllipses = TRUE, 
#              legend.title = "K-means Cluster") + 
#   theme_minimal()

# --- Final Interpretation Summary ---
# 1. Which PCs capture the most variance? (Scree plot)
# 2. Which original features drive these PCs? (Loadings/Biplot/Contribution plots)
# 3. How do samples group/separate on these PCs, and does it align with biological hypotheses? (Score plot)
# 4. Are there any unexpected outliers or patterns?

print("Interpretation complete. Proceed to generate hypotheses and validate them experimentally.")