# Load metabolomics data
metabo_data <- [Link]("metabolomics_abundance.csv")
# Make compound names unique and set as row names
metabo_data$Metabolite <- [Link]([Link](metabo_data$Metabolite))
rownames(metabo_data) <- metabo_data$Metabolite
metabo_data <- metabo_data[, -1]
# Remove "X" prefix from sample names if present
colnames(metabo_data) <- sub("^X", "", colnames(metabo_data))
# Load metadata
metadata <- [Link]("[Link]", header = TRUE, [Link] = 1, sep = "\t")
# Clean up sample names
colnames(metabo_data) <- trimws(colnames(metabo_data))
rownames(metadata) <- trimws(rownames(metadata))
# Reorder metadata to match columns of metabo_data (preserve all columns)
metadata <- metadata[colnames(metabo_data), , drop = FALSE] # THAT'S IT, WE
SHOULD Use drop = FALSE to preserve data structure even if metadata is a single
column
# Check alignment
stopifnot(all(colnames(metabo_data) == rownames(metadata)))
# 1. Remove metabolites with too many missing values
missing_threshold <- 0.2
keep_metabolites <- rowSums([Link](metabo_data)) < ncol(metabo_data) *
missing_threshold
metabo_data_filtered <- metabo_data[keep_metabolites, ]
# 2. Keep those detected in at least 3 samples
min_samples <- 3
keep_detected <- rowSums(metabo_data_filtered > 0) >= min_samples
metabo_data_filtered <- metabo_data_filtered[keep_detected, ]
# 3. Keep top 75% high-variance metabolites
metabo_vars <- apply(metabo_data_filtered, 1, var, [Link] = TRUE)
keep_var <- metabo_vars > quantile(metabo_vars, 0.25, [Link] = TRUE)
metabo_data_filtered <- metabo_data_filtered[keep_var, ]
# Log2 transform
metabo_data_log <- log2(metabo_data_filtered + 1)
# Impute remaining missing values with min/2
metabo_data_log[[Link](metabo_data_log)] <- min(metabo_data_log, [Link] = TRUE) / 2
# Transpose and scale (samples as rows)
metabo_data_scaled <- scale(t(metabo_data_log), center = TRUE, scale = TRUE)
# Perform PCA
pca_result <- prcomp(metabo_data_scaled, center = TRUE, scale. = FALSE)
percent_var <- round(100 * pca_result$sdev^2 / sum(pca_result$sdev^2), 1)
# Create PCA data frame with all metadata
pca_data <- [Link](pca_result$x)
pca_data <- cbind(pca_data, metadata[rownames(pca_data), ])
# Ensure categorical variables are clean
pca_data$FrostCondition <- factor(gsub(" ", "_", pca_data$FrostCondition), levels =
c("Before_frost", "After_frost"))
pca_data$TissueType <- factor(gsub(" ", "_", pca_data$TissueType))
pca_data$Rootstock <- factor(pca_data$Rootstock)
# Create PCA plot with jitter and reordered legends
library(ggplot2)
[Link](123)
ggplot(pca_data, aes(x = PC1, y = PC2, color = TissueType, shape = FrostCondition))
+
geom_point(aes(alpha = Rootstock), size = 4, position = position_jitter(width =
0.5, height = 0.5)) + # Add jitter and alpha for Rootstock
scale_alpha_manual(values = c("M26" = 0.3, "B9" = 1), name = "Rootstock") + #
Legend for Rootstock
scale_shape_manual(values = c("Before_frost" = 16, "After_frost" = 17)) + #
Circles and triangles
xlab(paste0("PC1 (", percent_var[1], "% variance)")) +
ylab(paste0("PC2 (", percent_var[2], "% variance)")) +
ggtitle("PCA of Metabolomics Data") +
theme_minimal(base_size = 14) +
scale_color_brewer(palette = "Set1") +
theme(
[Link] = element_rect(fill = "white", color = NA),
[Link] = element_rect(fill = "white", color = NA)
) +
guides(shape = guide_legend(title = "FrostCondition", order = 1), #
FrostCondition first
alpha = guide_legend(title = "Rootstock", order = 2), # Rootstock
second
color = guide_legend(title = "TissueType", order = 3))
ggsave("[Link]", width = 8, height = 6, dpi = 300, bg = "white")
# Save processed data
[Link](metabo_data_log, "[Link]")
[Link](metabo_data_scaled, "[Link]")
[Link](pca_data, "[Link]")