# ============================================================
# 教程配套脚本 05：可视化（火山图、热图、PCA）
# 输入：limma 结果 + 表达矩阵
# 输出：figures/ 下的 PNG 图
# ============================================================

library(ggplot2)
library(pheatmap)
library(tidyverse)

results <- read.csv("results/limma_all_results.csv")

# ---- 1. 火山图 ----
results <- results %>%
  mutate(sig = case_when(
    adj.P.Val < 0.05 & logFC > 1  ~ "up",
    adj.P.Val < 0.05 & logFC < -1 ~ "down",
    TRUE                          ~ "ns"
  ))

# 标注几个感兴趣的基因
label_genes <- c("GSDMD", "GSDME", "CASP1", "IL1B", "NLRP3")

p_volcano <- ggplot(results, aes(x = logFC, y = -log10(adj.P.Val), color = sig)) +
  geom_point(size = 0.6, alpha = 0.6) +
  scale_color_manual(values = c(up = "#C44E52", down = "#4C72B0", ns = "#BBBBBB"),
                     labels = c(up = "上调", down = "下调", ns = "无显著差异")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
  ggrepel::geom_text_repel(
    data = filter(results, gene %in% label_genes),
    aes(label = gene), size = 3, color = "black", max.overlaps = 20) +
  labs(x = "log2(Fold Change)", y = "-log10(adj.P.Val)",
       title = "Volcano plot: Tumor vs Normal") +
  theme_minimal() +
  theme(legend.title = element_blank())
ggsave("figures/R_volcano.png", p_volcano, width = 7, height = 5.5, dpi = 300)

# ---- 2. 热图（差异焦亡基因 × 部分样本）----
expr_symbol <- read.csv("data/processed/GSE14520_expr_symbol.csv",
                        row.names = 1, check.names = FALSE)
deg_pyro <- read.csv("results/DEGs_pyroptosis.csv")

# 提取差异焦亡基因表达，取前 40 个样本演示
mat <- expr_symbol[rownames(expr_symbol) %in% deg_pyro$gene, ]
mat <- mat[, 1:min(40, ncol(mat))]

# 行 z-score 标准化（热图常见预处理）
mat_z <- t(scale(t(mat)))

# 分组注释
annotation_col <- data.frame(Group = group[colnames(mat)])
rownames(annotation_col) <- colnames(mat)

pheatmap(mat_z,
         cluster_cols = TRUE, cluster_rows = TRUE,
         annotation_col = annotation_col,
         show_colnames = FALSE,
         color = colorRampPalette(c("#4C72B0", "white", "#C44E52"))(100),
         main = "Pyroptosis-related DEGs (z-score)",
         filename = "figures/R_heatmap.png", width = 8, height = 6)

# ---- 3. PCA 图 ----
library(ggrepel)
# 用全部基因或高变基因做 PCA（此处用全部表达矩阵示例）
pca_res <- prcomp(t(expr_symbol[, keep]), scale. = TRUE)
pca_df <- data.frame(PC1 = pca_res$x[, 1], PC2 = pca_res$x[, 2],
                     Group = group_use)

p_pca <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Group)) +
  geom_point(size = 2.5, alpha = 0.8) +
  stat_ellipse(level = 0.95) +
  labs(x = paste0("PC1 (", round(summary(pca_res)$importance[2, 1] * 100, 1), "%)"),
       y = paste0("PC2 (", round(summary(pca_res)$importance[2, 2] * 100, 1), "%)"),
       title = "PCA: Tumor vs Normal") +
  theme_minimal()
ggsave("figures/R_pca.png", p_pca, width = 6.5, height = 5.5, dpi = 300)

cat("可视化完成！图已保存到 figures/ \n")
