Metabolomics Statistical Analysis

Univariate Analysis

library(tidyverse)

Load normalized data

data <- read.csv('normalized_data.csv', row.names = 1) groups <- factor(read.csv('sample_info.csv')$group)

T-test for each feature

ttest_results <- apply(data, 2, function(x) { test <- t.test(x ~ groups) c(pvalue = test$p.value, fc = mean(x[groups == levels(groups)[2]]) - mean(x[groups == levels(groups)[1]])) }) ttest_results <- as.data.frame(t(ttest_results)) ttest_results$fdr <- p.adjust(ttest_results$pvalue, method = 'BH')

Significant features

sig_features <- ttest_results[ttest_results$fdr < 0.05, ] cat('Significant features (FDR<0.05):', nrow(sig_features), '\n')

Fold Change Calculation

Calculate fold change between groups

calculate_fc <- function(data, groups) { group_means <- aggregate(data, by = list(groups), FUN = mean, na.rm = TRUE) rownames(group_means) <- group_means$Group.1 group_means <- group_means[, -1]

fc <- as.numeric(group_means[2, ]) / as.numeric(group_means[1, ])
log2fc <- log2(fc)

return(data.frame(feature = colnames(data), fold_change = fc, log2fc = log2fc))

}

fc_results <- calculate_fc(data, groups)

Volcano Plot

library(ggplot2)

Combine statistics

results <- merge(ttest_results, fc_results, by.x = 'row.names', by.y = 'feature') results$significant <- results$fdr < 0.05 & abs(results$log2fc) > 1

Plot

ggplot(results, aes(x = log2fc, y = -log10(pvalue), color = significant)) + geom_point(alpha = 0.6) + scale_color_manual(values = c('gray', 'red')) + geom_hline(yintercept = -log10(0.05), linetype = 'dashed') + geom_vline(xintercept = c(-1, 1), linetype = 'dashed') + labs(x = 'Log2 Fold Change', y = '-log10(p-value)', title = 'Metabolomics Volcano Plot') + theme_bw()

ggsave('volcano_plot.png', width = 8, height = 6)

PCA

library(pcaMethods)

PCA

pca_result <- pca(data, nPcs = 5, method = 'ppca')

Scores plot

scores <- as.data.frame(scores(pca_result)) scores$group <- groups

ggplot(scores, aes(x = PC1, y = PC2, color = group)) + geom_point(size = 3) + stat_ellipse(level = 0.95) + labs(x = paste0('PC1 (', round(pca_result@R2[1] * 100, 1), '%)'), y = paste0('PC2 (', round(pca_result@R2[2] * 100, 1), '%)')) + theme_bw()

ggsave('pca_plot.png', width = 8, height = 6)

Loadings

loadings <- as.data.frame(loadings(pca_result)) top_pc1 <- loadings[order(abs(loadings$PC1), decreasing = TRUE)[1:20], ]

PLS-DA

library(mixOmics)

PLS-DA

plsda_result <- plsda(as.matrix(data), groups, ncomp = 3)

Cross-validation

perf_plsda <- perf(plsda_result, validation = 'Mfold', folds = 5, nrepeat = 50) plot(perf_plsda, col = color.mixo(5:7))

Optimal components

ncomp_opt <- perf_plsda$choice.ncomp['BER', 'centroids.dist'] cat('Optimal components:', ncomp_opt, '\n')

Final model

final_plsda <- plsda(as.matrix(data), groups, ncomp = ncomp_opt)

Plot

plotIndiv(final_plsda, group = groups, ellipse = TRUE, legend = TRUE)

VIP scores

vip <- vip(final_plsda) top_vip <- sort(vip[, ncomp_opt], decreasing = TRUE)[1:20] print(top_vip)

sPLS-DA (Sparse)

Tune number of features to select

tune_splsda <- tune.splsda(as.matrix(data), groups, ncomp = 3, validation = 'Mfold', folds = 5, nrepeat = 50, test.keepX = c(5, 10, 20, 50, 100))

optimal_keepX <- tune_splsda$choice.keepX

Final sparse model

splsda_result <- splsda(as.matrix(data), groups, ncomp = ncomp_opt, keepX = optimal_keepX)

Selected features

selected_features <- selectVar(splsda_result, comp = 1)$name cat('Selected features (comp 1):', length(selected_features), '\n')

OPLS-DA (Orthogonal PLS-DA)

library(ropls)

OPLS-DA

oplsda <- opls(data, groups, predI = 1, orthoI = NA)

Summary

print(oplsda)

Scores plot

plot(oplsda, typeVc = 'x-score')

S-plot (loadings vs correlation)

plot(oplsda, typeVc = 'x-loading')

VIP

vip_scores <- getVipVn(oplsda) top_vip <- sort(vip_scores, decreasing = TRUE)[1:20]

Random Forest

library(randomForest)

Random Forest classification

rf_model <- randomForest(x = data, y = groups, importance = TRUE, ntree = 500)

Importance

importance <- importance(rf_model) top_features <- rownames(importance)[order(importance[, 'MeanDecreaseAccuracy'], decreasing = TRUE)[1:20]]

Plot importance

varImpPlot(rf_model, n.var = 20)

ROC Analysis

library(pROC)

ROC for top biomarker

top_feature <- 'feature_123' # Replace with actual feature name roc_result <- roc(groups, data[, top_feature])

Plot

plot(roc_result, main = paste('AUC =', round(auc(roc_result), 3)))

Multiple features

biomarkers <- c('feature_1', 'feature_2', 'feature_3') for (feat in biomarkers) { roc_i <- roc(groups, data[, feat]) cat(feat, ': AUC =', round(auc(roc_i), 3), '\n') }

Heatmap

library(pheatmap)

Top differential features

top_features <- rownames(sig_features)[1:50] data_top <- data[, top_features]

Annotation

annotation_row <- data.frame(Group = groups) rownames(annotation_row) <- rownames(data)

pheatmap(t(data_top), annotation_col = annotation_row, scale = 'row', clustering_method = 'ward.D2', filename = 'heatmap.png', width = 10, height = 12)

Related Skills

• normalization-qc - Data preparation

• pathway-mapping - Functional interpretation

• multi-omics-integration/mixomics-analysis - Advanced multivariate methods