Untitled

# ================================
# Publication-ready Maaslin2 plots
# ================================

# Load packages
library(readr)
library(dplyr)
library(ggplot2)
library(pheatmap)
library(tidyr)
library(forcats)

# ----------------
# Load results
# ----------------
res <- read_tsv("maaslin2_out/all_results.tsv")

# Load your metadata and abundance (adjust filenames if needed)
metadata <- read_tsv("your_metadata.tsv")
abundance <- read_tsv("your_abundance.tsv") %>%
  column_to_rownames("feature")

# Filter significant results (q ≤ 0.05)
sig_res <- res %>%
  filter(qval <= 0.05)

# ----------------
# Panel A: Volcano plot
# ----------------
p_volcano <- ggplot(res, aes(x = coef, y = -log10(qval), color = qval < 0.05)) +
  geom_point(alpha = 0.7, size = 2) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "red") +
  theme_minimal(base_size = 14) +
  scale_color_manual(values = c("grey70", "firebrick")) +
  labs(x = "Effect size (coef)", y = "-log10(q-value)", color = "Significant")

# ----------------
# Panel B: Barplot of top taxa
# ----------------
top_sig <- sig_res %>%
  arrange(qval) %>%
  slice_head(n = 20)

p_bar <- ggplot(top_sig, aes(x = fct_reorder(feature, coef), y = coef, fill = metadata)) +
  geom_col() +
  coord_flip() +
  theme_minimal(base_size = 14) +
  labs(x = "Feature", y = "Effect size (coef)", fill = "Metadata")

# ----------------
# Panel C: Heatmap of significant taxa
# ----------------
sig_features <- unique(sig_res$feature)
mat <- abundance[sig_features, ]

# Make sure samples line up
mat <- mat[, rownames(metadata), drop = FALSE]

pheatmap(mat,
         scale = "row",
         clustering_distance_rows = "euclidean",
         clustering_distance_cols = "euclidean",
         annotation_col = metadata %>% column_to_rownames("sample_id"),
         fontsize = 10,
         main = "Heatmap of significant features")

# ----------------
# Panel D: Boxplots of top taxa
# ----------------
# Take top 5 significant features
top5 <- top_sig %>%
  pull(feature) %>%
  unique() %>%
  head(5)

# Convert abundance to long format
ab_long <- abundance %>%
  rownames_to_column("feature") %>%
  filter(feature %in% top5) %>%
  pivot_longer(-feature, names_to = "sample_id", values_to = "abundance") %>%
  left_join(metadata, by = "sample_id")

p_box <- ggplot(ab_long, aes(x = group, y = abundance, fill = group)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.6, size = 1) +
  facet_wrap(~ feature, scales = "free_y") +
  theme_minimal(base_size = 14) +
  labs(x = "", y = "Abundance", fill = "Group")

# ----------------
# Save plots
# ----------------
ggsave("PanelA_Volcano.png", p_volcano, width = 6, height = 4, dpi = 300)
ggsave("PanelB_Barplot.png", p_bar, width = 6, height = 6, dpi = 300)
ggsave("PanelD_Boxplots.png", p_box, width = 8, height = 6, dpi = 300)
# Heatmap is saved manually by pheatmap if you add filename="PanelC_Heatmap.png"