# ============================================================
# 02_differential_analysis.R
# 差异表达分析（limma）
# ============================================================

library(limma)
library(ggplot2)
library(ggrepel)
library(pheatmap)

# ---- 读取数据 ----
expr <- read.csv("data/GSE38974_COPD_expression_clean.csv", 
                 row.names = 1, check.names = FALSE)
sample_info <- read.csv("data/GSE38974_COPD_sample_info.csv")

# 确保分组是因子
sample_info$Group <- factor(sample_info$Group, levels = c("Control", "COPD"))
expr <- expr[, sample_info$Sample]

# ---- limma差异分析 ----
design <- model.matrix(~ 0 + sample_info$Group)
colnames(design) <- levels(sample_info$Group)

contrast.matrix <- makeContrasts(
  COPD_vs_Control = COPD - Control,
  levels = design
)

fit <- lmFit(expr, design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

# 提取结果
DEGs <- topTable(fit2, adjust.method = "fdr", number = Inf, sort.by = "B")

# 添加差异分组
DEGs$group <- case_when(
  DEGs$logFC > 1 & DEGs$adj.P.Val < 0.05 ~ "Up",
  DEGs$logFC < -1 & DEGs$adj.P.Val < 0.05 ~ "Down",
  TRUE ~ "NotSig"
)

cat("差异基因统计:\n")
print(table(DEGs$group))

# 保存
write.csv(DEGs, "output/COPD_DEGs_all.csv")
deg_sig <- DEGs[DEGs$group != "NotSig", ]
write.csv(deg_sig, "output/COPD_DEGs_significant.csv")

# ---- 火山图 ----
plot_data <- DEGs
plot_data$negLogP <- -log10(plot_data$adj.P.Val)

top_up <- head(rownames(plot_data[plot_data$group == "Up", ]), 10)
top_down <- head(rownames(plot_data[plot_data$group == "Down", ]), 10)
plot_data$label <- ifelse(rownames(plot_data) %in% c(top_up, top_down),
                          rownames(plot_data), "")

ggplot(plot_data, aes(x = logFC, y = negLogP, color = group)) +
  geom_point(alpha = 0.6, size = 1.5) +
  scale_color_manual(values = c("Down" = "#2E86AB", 
                                 "NotSig" = "grey70", 
                                 "Up" = "#F24236")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
  geom_text_repel(aes(label = label), size = 3, max.overlaps = 30) +
  labs(x = "log2 Fold Change", y = "-log10(FDR)",
       title = "COPD vs Control Volcano Plot") +
  theme_minimal()

ggsave("output/04_volcano_plot.png", width = 10, height = 8, dpi = 300)

# ---- 热图 ----
top50_genes <- rownames(head(DEGs[order(DEGs$adj.P.Val), ], 50))
expr_top50 <- expr[top50_genes, ]

annotation_col <- data.frame(
  Group = factor(sample_info$Group),
  row.names = sample_info$Sample
)

expr_scaled <- t(scale(t(expr_top50)))

pheatmap(expr_scaled,
         annotation_col = annotation_col,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         show_colnames = FALSE,
         show_rownames = TRUE,
         fontsize_row = 6,
         main = "Top 50 DEGs Heatmap",
         filename = "output/05_deg_heatmap.png",
         width = 10, height = 12)

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