# ============================================================
# 教程配套脚本 04：差异表达分析（limma）
# 输入：预处理后的表达矩阵 + 分组
# 输出：差异基因表 + 火山图/热图/PCA 数据
# ============================================================

library(limma)
library(tidyverse)

expr_symbol <- read.csv("data/processed/GSE14520_expr_symbol.csv",
                        row.names = 1, check.names = FALSE)
group <- readRDS("data/processed/GSE14520_group.rds")

# ---- 1. 只保留肿瘤与正常两组样本 ----
keep <- group %in% c("Tumor", "Normal")
expr_use <- expr_symbol[, keep]
group_use <- factor(group[keep], levels = c("Normal", "Tumor"))

# ---- 2. 设计矩阵 ----
design <- model.matrix(~ group_use)
colnames(design) <- c("Intercept", "Tumor_vs_Normal")

# ---- 3. limma 三件套 ----
fit <- lmFit(expr_use, design)
fit <- eBayes(fit)

# 提取全部基因的差异结果
results <- topTable(fit, coef = "Tumor_vs_Normal", number = Inf, sort.by = "none")
head(results)

# ---- 4. 筛选差异基因：|log2FC|>1 且 adj.P<0.05 ----
deg <- results %>%
  rownames_to_column("gene") %>%
  filter(abs(logFC) > 1 & adj.P.Val < 0.05)

up <- deg %>% filter(logFC > 0) %>% nrow()
down <- deg %>% filter(logFC < 0) %>% nrow()
cat("差异基因总数:", nrow(deg), " 上调:", up, " 下调:", down, "\n")

# ---- 5. 保存结果 ----
write.csv(results, "results/limma_all_results.csv", row.names = FALSE)
write.csv(deg, "results/DEGs_tumor_vs_normal.csv", row.names = FALSE)

# ---- 6. 与焦亡基因集取交集（示例列表，正式研究用文献完整列表）----
pyroptosis_genes <- c("GSDMD", "GSDME", "GSDMA", "GSDMB", "GSDMC",
                      "CASP1", "CASP3", "CASP4", "CASP5", "CASP8",
                      "IL1B", "IL18", "NLRP3", "AIM2", "PYCARD",
                      "IL1A", "IL6", "TNF", "NFKB1", "NFKB2", "PJVK")
deg_pyro <- deg %>% filter(gene %in% pyroptosis_genes)
cat("差异焦亡基因数:", nrow(deg_pyro), "\n")
print(deg_pyro$gene)

write.csv(deg_pyro, "results/DEGs_pyroptosis.csv", row.names = FALSE)

cat("差异表达分析完成！\n")
