{r setup, include=FALSE} knitr::opts_chunk_set(echo = TRUE, warning = FALSE, message = FALSE, fig.width = 10, fig.height = 7) library(readxl) library(ggplot2) library(tidyr) if(!require(patchwork)) install.packages(“patchwork”) library(patchwork) if(!require(ggcorrplot)) install.packages(“ggcorrplot”) library(ggcorrplot) if (!requireNamespace(“BiocManager”, quietly = TRUE)) install.packages(“BiocManager”) if (!requireNamespace(“mixOmics”, quietly = TRUE)) BiocManager::install(“mixOmics”) library(mixOmics) if(!require(openxlsx)) install.packages(“openxlsx”) library(openxlsx)

1. Data Import and Quality Control Summary

load the raw datasets and perform a basic check on total missing entries.

Load original dataset

my_data <- read_excel(“PharmaDx_metabolomics.xlsx”)

Basic dataset architecture metrics

n_samples_init <- nrow(my_data) n_features_init <- ncol(my_data) - 6 # Subtracting metadata columns (SampleID, Group, Age, BMI, Gender, Batch) total_nas_init <- sum(is.na(my_data))

Basic inspection

dim(my_data) head(my_data) str(my_data) summary(my_data)

Check group distribution

table(my_data$Group)

Check batch distribution

table(my_data$Batch)

2. Preprocessing: Filtering, Half-Min Imputation, Log2 & Autoscaling

Filter out metabolites with > 20% missing values directly

bad_metabolites <- colnames(my_data)[colMeans(is.na(my_data)) > 0.20] my_data_clean <- my_data[, !(colnames(my_data) %in% bad_metabolites)]

Isolate metabolite matrix (Column 7 onwards)

metabolite_data <- my_data_clean[, 7:ncol(my_data_clean)]

Imputation using Half-Minimum method

for (col in colnames(metabolite_data)) { if (any(is.na(metabolite_data[[col]]))) { min_val <- min(metabolite_data[[col]], na.rm = TRUE) metabolite_data[[col]][is.na(metabolite_data[[col]])] <- min_val / 2 } }

Log2(x + 1) Transformation

metabolite_data_log2 <- log2(metabolite_data + 1)

Autoscaling (Z-score)

metabolite_data_autoscaled <- scale(metabolite_data_log2)

Recombine metadata with treated frames

my_data_log2 <- cbind(my_data_clean[, 1:6], metabolite_data_log2) my_data_autoscaled_final <- cbind(my_data_clean[, 1:6], as.data.frame(metabolite_data_autoscaled))

Verification Boxplots (Glucose) tracking pre- and post-normalization steps

plot_before <- ggplot(my_data_clean, aes(x = Group, y = Glucose, fill = Group)) + geom_boxplot(alpha = 0.7, outlier.shape = NA) + geom_jitter(width = 0.15, alpha = 0.5) + labs(title = “Glucose: Raw Data”, y = “Raw Intensity”, x = ““) + theme_minimal() + theme(legend.position =”none”)

plot_after <- ggplot(my_data_autoscaled_final, aes(x = Group, y = Glucose, fill = Group)) + geom_boxplot(alpha = 0.7, outlier.shape = NA) + geom_jitter(width = 0.15, alpha = 0.5) + labs(title = “Glucose: Autoscaled Data”, y = “Z-Score”, x = ““) + theme_minimal() + theme(legend.position =”none”)

plot_before + plot_after

3. Exploratory Analysis: PCA Unsupervised Assessment

pca_result <- prcomp(metabolite_data_autoscaled, scale. = FALSE) var_explained <- round(100 * (pca_result$sdev^2 / sum(pca_result$sdev^2)), 1)

pca_df <- data.frame( PC1 = pca_result$x[, 1], PC2 = pca_result$x[, 2], Batch = as.factor(my_data_clean$Batch), Group = as.factor(my_data_clean$Group) )

Plot A: QC Check Colored by Batch

pca_batch <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Batch)) + geom_point(size = 3.5, alpha = 0.8) + labs(title = “QC Assessment: PCA by Run Batch”, x = paste0(“PC1 (”, var_explained[1], “%)”), y = paste0(“PC2 (”, var_explained[2], “%)”)) + theme_minimal()

Plot B: Separation Colored by Biological Status

pca_group <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Group)) + geom_point(size = 3.5, alpha = 0.8) + scale_color_manual(values = c(“Healthy” = “#2ca02c”, “PreDiabetic” = “#ff7f0e”, “T2D” = “#d62728”)) + labs(title = “Biological Clustering: PCA by Clinical Group”, x = paste0(“PC1 (”, var_explained[1], “%)”), y = paste0(“PC2 (”, var_explained[2], “%)”)) + theme_minimal()

pca_batch / pca_group

4. Differential Expression: One-Way ANOVA with BH Correction

metabolite_cols <- colnames(metabolite_data) anova_p_values <- numeric(length(metabolite_cols)) names(anova_p_values) <- metabolite_cols

for (met in metabolite_cols) { model <- lm(metabolite_data_autoscaled[, met] ~ my_data_clean$Group) anova_res <- summary(aov(model)) anova_p_values[met] <- anova_res[[1]]["my_data_clean$Group”, “Pr(>F)”] }

bh_adjusted_p <- p.adjust(anova_p_values, method = “BH”) anova_summary <- data.frame(Raw_P_Value = anova_p_values, BH_adjusted_p = bh_adjusted_p)

Faceted Distributions (Top 10 Significant Biomarkers)

sorted_anova <- anova_summary[order(anova_summary[, 2]), ] # Sort by the second column top_10_names <- head(rownames(sorted_anova), 10)

top_10_data <- my_data_autoscaled_final[, c(“Group”, top_10_names)] top_10_long <- pivot_longer(top_10_data, cols = -Group, names_to = “Metabolite”, values_to = “Z_Score”) top_10_long$Metabolite <- factor(top_10_long$Metabolite, levels = top_10_names)

ggplot(top_10_long, aes(x = Group, y = Z_Score, fill = Group)) + geom_boxplot(alpha = 0.7, outlier.shape = NA) + geom_jitter(width = 0.15, alpha = 0.3, size = 1) + facet_wrap(~Metabolite, scales = “free_y”, ncol = 5) +
scale_fill_manual(values = c(“Healthy” = “#2ca02c”, “PreDiabetic” = “#ff7f0e”, “T2D” = “#d62728”)) + labs(title = “Top 10 Most Significant Metabolites Across Clinical States”, y = “Z-Score”, x = ““) + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position =”none”)

5. Biomarker Metrics: Fold Changes and Volcano Modeling

means_healthy <- colMeans(my_data_log2[my_data_log2$Group == "Healthy", metabolite_cols], na.rm = TRUE) means_prediabetic <- colMeans(my_data_log2[my_data_log2$Group == “PreDiabetic”, metabolite_cols], na.rm = TRUE) means_t2d <- colMeans(my_data_log2[my_data_log2$Group == “T2D”, metabolite_cols], na.rm = TRUE)

log2fc_prediabetic_vs_healthy <- means_prediabetic - means_healthy log2fc_t2d_vs_healthy <- means_t2d - means_healthy

fc_summary <- data.frame( Log2FC_PreDiabetic_vs_Healthy = log2fc_prediabetic_vs_healthy, Log2FC_T2D_vs_Healthy = log2fc_t2d_vs_healthy )

Pairwise individual t-tests for T2D vs Healthy volcano display

Create subsets for just the T2D and Healthy groups (using Log2 data)

t2d_sub <- my_data_log2[my_data_log2$Group == "T2D", metabolite_cols] healthy_sub <- my_data_log2[my_data_log2$Group == “Healthy”, metabolite_cols]

Calculate individual t-test p-values for each metabolite

t_test_p <- numeric(length(metabolite_cols)) names(t_test_p) <- metabolite_cols for (met in metabolite_cols) { t_test_p[met] <- t.test(t2d_sub[[met]], healthy_sub[[met]])$p.value }

Apply Benjamini-Hochberg multiple testing correction

t2d_bh_p <- p.adjust(t_test_p, method = “BH”)

Build the Volcano Plot data frame

volcano_df <- data.frame( Log2FC = fc_summary$Log2FC_T2D_vs_Healthy, P_Value = t2d_bh_p, NegLogP = -log10(t2d_bh_p) )

Categorize metabolites into groups for coloring

volcano_df$Significance <- "Not Significant" volcano_df$Significance[volcano_df$Log2FC > 0.58 & volcano_df$P_Value < 0.05] <- “Up-regulated in T2D” volcano_df$Significance[volcano_df$Log2FC < -0.58 & volcano_df$P_Value < 0.05] <- “Down-regulated in T2D”

Volcano Plot: T2D vs Healthy

ggplot(volcano_df, aes(x = Log2FC, y = NegLogP, color = Significance)) + geom_point(alpha = 0.8, size = 3) + geom_hline(yintercept = -log10(0.05), linetype = “dashed”, color = “gray40”) + geom_vline(xintercept = c(-0.58, 0.58), linetype = “dashed”, color = “gray40”) + scale_color_manual(values = c(“Not Significant” = “gray70”, “Up-regulated in T2D” = “#d62728”, “Down-regulated in T2D” = “#1f77b4”)) + labs(title = “Volcano Mapping: T2D vs Healthy Cohorts”, x = “Log2 Fold Change”, y = “-log10(BH-Adjusted P-Value)”) + theme_minimal() + theme(legend.position = “bottom”)

6. Supervised Modeling: PLS-DA and VIP Analysis

Extract your clean metabolite and biological group outcome factor

X <- metabolite_data_autoscaled Y <- as.factor(my_data_clean$Group)

Run the PLS-DA model with 2 components

plsda_model <- plsda(X, Y, ncomp = 2)

Generate PLS-DA score plot

plotIndiv(plsda_model, comp = c(1, 2), group = Y, ind.names = FALSE, ellipse = TRUE, legend = TRUE, title = “PLS-DA Separation Grid”, col.per.group = c(“#2ca02c”, “#ff7f0e”, “#d62728”))

Calculate VIP scores from your mixOmics PLS-DA model

vip_matrix <- vip(plsda_model)

Put the results into a clean dataframe

vip_summary1 <- data.frame(Metabolite = rownames(vip_matrix), VIP_Comp1 = vip_matrix[, 1])

Sort the metabolites from highest importance to lowest importance

vip_summary2 <- vip_summary1[order(-vip_summary1$VIP_Comp1), ]

7. Secondary Association Networks: Co-regulation Heatmap

Calculate the Pearson correlation matrix between columns

corr_matrix <- cor(metabolite_data_autoscaled, method = “pearson”, use = “complete.obs”)

Generate the correlation heatmap

ggcorrplot( corr_matrix, hc.order = TRUE, # Hierarchical clustering (groups similar metabolites together) type = “full”, # Shows the full square matrix outline.col = “white”, # Clean white borders between blocks colors = c(“#1f77b4”, “white”, “#d62728”), # Blue (negative), White (none), Red (positive) title = “Global Covariation Feature Matrix”, lab = FALSE # Set to TRUE only if you want numbers inside squares ) + theme( plot.title = element_text(face = “bold”, size = 14, hjust = 0.5), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 8), # Vertical labels axis.text.y = element_text(size = 8) )

📂 Step 8: Save Your Master Sheet into Excel

Align variables by row names cleanly to prevent alignment shift errors

final_master_sheet <- data.frame( Metabolite = rownames(anova_summary), Raw_ANOVA_P = anova_summary$Raw_P_Value, BH_Adjusted_P = anova_summary$BH_adjusted_p, Log2FC_PreD_vs_H = fc_summary[rownames(anova_summary), “Log2FC_PreDiabetic_vs_Healthy”], Log2FC_T2D_vs_H = fc_summary[rownames(anova_summary), “Log2FC_T2D_vs_Healthy”], VIP_Score_Comp1 = vip_summary2[rownames(anova_summary), “VIP_Comp1”] )

Export straight into Task-1 folder

write.xlsx(final_master_sheet, “PharmaDx_Metabolomics_Final_Statistics.xlsx”, rowNames = FALSE)

Task1.2

Augustine

2026-07-26