knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
knitr::opts_chunk$set(fig.width = 10, fig.height = 7, dpi = 300)
packages <- c("tidyverse", "pheatmap", "limma", "mixOmics", "ggrepel","ggpubr","randomForest","iml","caret","nnet")
lapply(packages, function(p) if (!requireNamespace(p, quietly = TRUE)) install.packages(p))
lapply(packages, library, character.only = TRUE)
metadata <- read.csv("Metadata.csv")
metagenomics <- read.csv("Metagenomics.csv", row.names = 1)
metatranscriptomics <- read.csv("Metatranscriptomics.csv", row.names = 1)
metaproteomics <- read.csv("Metaproteomics.csv", row.names = 1)
metabolomics <- read.csv("Metabolomics.csv", row.names = 1)
# === STEP 1: Clean sample names ===
clean_colnames <- function(df) {
colnames(df) <- gsub("\\.", " ", colnames(df))
return(df)
}
metagenomics <- clean_colnames(metagenomics)
metatranscriptomics <- clean_colnames(metatranscriptomics)
metaproteomics <- clean_colnames(metaproteomics)
metabolomics <- clean_colnames(metabolomics)
# === STEP 2: Align samples ===
common_ids <- Reduce(intersect, list(
colnames(metagenomics),
colnames(metatranscriptomics),
colnames(metaproteomics),
colnames(metabolomics),
metadata$SampleID
))
metadata <- metadata %>%
filter(SampleID %in% common_ids) %>%
arrange(SampleID)
rownames(metadata) <- metadata$SampleID
omics_list <- list(
Metagenomics = metagenomics[, metadata$SampleID],
Metatranscriptomics = metatranscriptomics[, metadata$SampleID],
Metaproteomics = metaproteomics[, metadata$SampleID],
Metabolomics = metabolomics[, metadata$SampleID]
)
# === STEP 3: Transpose + Fix rownames (PLACE THIS HERE!) ===
omics_list <- lapply(omics_list, function(df) {
df_t <- as.data.frame(t(df))
df_t$SampleID <- rownames(df_t)
rownames(df_t) <- df_t$SampleID
df_t$SampleID <- NULL
return(df_t)
})
plot_heatmap <- function(data, stats, title) {
top_feats <- rownames(stats %>% arrange(FDR) %>% head(20))
mat <- data[, top_feats]
annotation <- metadata[rownames(mat), , drop = FALSE]
pheatmap(mat, annotation_row = annotation["Group"],
main = paste("Top Features -", title))
}
plot_volcano <- function(stats, title) {
stats$Feature <- rownames(stats)
ggplot(stats, aes(x = logFC, y = -log10(pval))) +
geom_point(aes(color = FDR < 0.05), alpha = 0.6) +
geom_text_repel(data = subset(stats, FDR < 0.01),
aes(label = Feature), size = 3, max.overlaps = 10) +
scale_color_manual(values = c("grey", "red")) +
theme_minimal() + labs(title = paste("Volcano Plot -", title))
}
plot_bar <- function(stats, title) {
top <- stats %>% arrange(FDR) %>% head(10)
top$Feature <- factor(rownames(top), levels = rownames(top))
ggplot(top, aes(x = Feature, y = logFC, fill = FDR < 0.05)) +
geom_col() +
coord_flip() +
theme_minimal() +
labs(title = paste("Top Differential Features -", title))
}
plot_pca <- function(data, title) {
# Remove columns with all NA values
data <- data[, colSums(is.na(data)) < nrow(data), drop = FALSE]
# Remove rows with any NA values
data <- data[complete.cases(data), , drop = FALSE]
# Check if there's enough data left to perform PCA
if (nrow(data) < 2 || ncol(data) < 2) {
warning(paste("Skipping PCA for", title, "- insufficient data after NA filtering."))
return(NULL)
}
# Perform PCA
pca <- prcomp(data, scale. = FALSE)
# Get PCA scores and add Group info
scores <- as.data.frame(pca$x)
scores$Group <- metadata[rownames(scores), "Group"]
# Plot the PCA
ggplot(scores, aes(x = PC1, y = PC2, color = Group)) +
geom_point(size = 3, alpha = 0.8) +
theme_minimal() +
labs(title = paste("PCA -", title), x = "PC1", y = "PC2")
}
plot_top_boxplots <- function(data, stats, title, n = 6) {
top_feats <- rownames(stats %>% arrange(FDR) %>% head(n))
long_data <- data[, top_feats] %>%
rownames_to_column("SampleID") %>%
pivot_longer(-SampleID, names_to = "Feature", values_to = "Value") %>%
left_join(metadata, by = "SampleID")
ggplot(long_data, aes(x = Group, y = Value, fill = Group)) +
geom_boxplot(alpha = 0.7) +
stat_compare_means(
method = "t.test",
label = "p.signif",
comparisons = list(c("Case", "Control"))
) +
facet_wrap(~Feature, scales = "free_y") +
theme_minimal() +
labs(title = paste("Top", n, "Features by Boxplot -", title))
}
for (omics_name in names(omics_list)) {
cat("Generating plots for:", omics_name, "\n")
data <- omics_list[[omics_name]]
stats <- stats_results[[omics_name]]
plot_heatmap(data, stats, omics_name)
print(plot_volcano(stats, omics_name))
print(plot_bar(stats, omics_name))
print(plot_pca(data, omics_name))
print(plot_top_boxplots(data, stats, omics_name))
}
