# -*- encoding: utf-8 -*-
# ============================================================================
# bioinfo_rna_seq.R — 某细胞类型 RNA-seq 差异表达 + 富集分析模板（R / DESeq2）
# ----------------------------------------------------------------------------
# 技能定位（自 v2.x 起）：R 用于生信分析；投稿级绘图统一走 Python 引擎
# (prism_theme.py)，本脚本只负责分析产出（DEG 表 / 富集表 / 分析图）。
#
# 用法：
#   1) RStudio 打开本文件，逐段运行（或 source 全文件）；
#   2) 默认 demo=TRUE：用「真实人类基因名 + 合成计数」端到端跑通全流程，
#      不涉及任何真实数据，仅验证管线可用；
#   3) 接入真实数据：demo=FALSE，准备 counts.csv 与 colData.csv（见下方说明），
#      并调整分组/对比方向。
#
# 输出：./bio_out/  （DEG 表 CSV、MA/火山图/PCA、GO/KEGG/GSEA 表）
# 依赖：DESeq2, clusterProfiler, org.Hs.eg.db, ggplot2；缺失会自动提示安装。
# 数据安全：本脚本纯本地运行，无任何联网上传代码；真实数据只在本机被读取。
# ============================================================================

# ---------------------------------------------------------------------------
# 0. 环境自检（依赖守卫）：缺失即打印安装命令并停止，绝不伪造结果
# ---------------------------------------------------------------------------
# 关键区分（v2.9.5）：requireNamespace() 失败有两种完全不同的原因，必须分开报，
# 否则会把人引向错误的修法：
#   ABSENT —— 从未安装          → 装一次即可
#   BROKEN —— 目录在、却加载失败（安装残缺/被中断、R 大版本升级后未重装）
#             → 必须先 remove.packages() 再 force 重装，直接重装同一个包往往还是坏
need <- c("DESeq2", "clusterProfiler", "org.Hs.eg.db", "ggplot2", "ragg")  # ragg 供 agg_png 出图
inst <- rownames(installed.packages())
load_err <- list()
have <- vapply(need, function(p) {
  ok <- suppressWarnings(suppressMessages(requireNamespace(p, quietly = TRUE)))
  if (!ok) {
    # 抓真实加载错误（如 org.Hs.eg.db 的 "DB file '' not found"），供排错用
    load_err[[p]] <<- tryCatch({
      suppressWarnings(loadNamespace(p)); "未知原因（加载失败但未抛出可读错误）"
    }, error = function(e) gsub("[\r\n]+", " ", conditionMessage(e)))
  }
  ok
}, logical(1))

if (!all(have)) {
  miss    <- need[!have]
  absent  <- miss[!miss %in% inst]   # 压根没装
  broken  <- miss[miss %in% inst]    # 装了但坏
  cat("== R 依赖检查未通过（已中止，未运行分析，避免伪结果） ==\n")
  if (length(absent)) {
    cat(sprintf("MISSING_PKG: %s\n", paste(absent, collapse = ", ")))
    cat("  原因：这些包从未安装。\n")
    cat("  修复（在 RStudio 里逐行运行）：\n")
    cat("    install.packages('BiocManager')\n")
    cat(sprintf("    BiocManager::install(c('%s'))\n", paste(absent, collapse = "', '")))
  }
  if (length(broken)) {
    cat(sprintf("BROKEN_PKG: %s\n", paste(broken, collapse = ", ")))
    cat("  原因：包目录存在，但加载时报错——通常是安装残缺/中断，或 R 大版本升级后未重装。\n")
    cat("        **这不是「没装」，重复装同一个包大概率仍然是坏的。**\n")
    for (p in broken) cat(sprintf("    · %s → %s\n", p, load_err[[p]]))
    cat("  修复（在 RStudio 里逐行运行；先删干净再强制重装）：\n")
    cat(sprintf("    remove.packages(c('%s'))\n", paste(broken, collapse = "', '")))
    cat(sprintf("    BiocManager::install(c('%s'), force = TRUE)\n", paste(broken, collapse = "', '")))
    cat("    # 重装后自检：sapply(c('DESeq2','org.Hs.eg.db'), requireNamespace, quietly=TRUE)\n")
  }
  stop("依赖缺失：", paste(miss, collapse = ", "), "；已中止（未运行分析，避免伪结果）")
}
suppressMessages({
  library(DESeq2)
  library(clusterProfiler)
  library(org.Hs.eg.db)
  library(ggplot2)
})

# 输出目录
out_dir <- "bio_out"
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

# 图输出助手：PDF（cairo 矢量，投稿用）+ PNG（ragg 300dpi，预览用）。
# 图内文字一律用英文/ASCII，规避 headless Rscript 下中文字体与编码问题。
save_pdf_png <- function(p, stem, w = 6, h = 5, dpi = 300) {
  grDevices::cairo_pdf(paste0(stem, ".pdf"), width = w, height = h)
  print(p); grDevices::dev.off()
  ragg::agg_png(paste0(stem, ".png"), width = w, height = h, units = "in", res = dpi)
  print(p); grDevices::dev.off()
  cat("  [figure] ", basename(stem), ".pdf/.png\n", sep = "")
}

# ---------------------------------------------------------------------------
# 1. 参数与数据
# ---------------------------------------------------------------------------
demo <- TRUE   # TRUE=合成演示数据；FALSE=读取真实数据（接入真实课题时改这里）

if (demo) {
  # ---- 合成演示数据：真实人类基因名 + 负二项合成计数（仅验证管线，非真实数据）
  set.seed(42)
  demo_genes <- sample(keys(org.Hs.eg.db, keytype = "SYMBOL"), 5000)
  n_sg <- 3; n_nc <- 3                       # 每组样本数（可改，组内重复>=3）
  base_mean <- 2^runif(5000, 4, 12)          # 基线表达量（跨度大，贴近真实分布）
  # 约 200 个"真差异"基因：log2FC ≈ ±2，用于让 demo 产出真实 DEG
  de_idx <- sample(5000, 200)
  fc <- rep(1, 5000); fc[de_idx] <- ifelse(runif(200) < 0.5, 4, 0.25)
  counts <- t(sapply(seq_len(5000), function(g) {
    mu <- base_mean[g] * c(rep(1, n_nc), rep(fc[g], n_sg))   # Ctrl=1x, Treat=fc x
    rnbinom(n_sg + n_nc, mu = mu, size = 5)   # size 控制离散度
  }))
  rownames(counts) <- demo_genes
  colnames(counts) <- c(paste0("Ctrl", 1:n_nc), paste0("Treat", 1:n_sg))
  coldata <- data.frame(
    row.names = colnames(counts),
    group     = factor(rep(c("Ctrl", "Treat"), c(n_nc, n_sg)), levels = c("Ctrl", "Treat"))
  )
  cat(sprintf("[demo] 合成数据: %d 基因 x %d 样本 (处理组 vs 对照组)\n",
              nrow(counts), ncol(counts)))
} else {
  # ---- 真实数据模式（接入真实课题时使用）----
  # counts.csv  ：行=基因（第一列为基因名，建议用 SYMBOL 或 ENTREZID），列=样本
  # colData.csv ：行名=样本名（与 counts 列名一致），含 group 列（处理组/对照组 等）
  counts  <- as.matrix(read.csv("counts.csv",  row.names = 1, check.names = FALSE))
  coldata <- read.csv("colData.csv", row.names = 1, check.names = FALSE)
  coldata$group <- factor(coldata$group)
  cat(sprintf("[real] 读取数据: %d 基因 x %d 样本\n", nrow(counts), ncol(counts)))
}

# ---------------------------------------------------------------------------
# 2. DESeq2 差异表达（处理组 vs 对照组）
# ---------------------------------------------------------------------------
dds <- DESeqDataSetFromMatrix(
  countData = round(counts),          # DESeq2 要求整数计数
  colData   = coldata,
  design    = ~group
)
dds <- DESeq(dds)

# 对比方向：处理组 vs 对照组（log2FC>0 = 处理组上调）。若课题定义相反，交换两水平即可。
res <- results(dds, contrast = c("group", "Treat", "Ctrl"))
cat("\n== 差异表达汇总（处理组 vs 对照组）==\n")
print(summary(res))

# 缩 shrinkage 版的 log2FC 更适合排序/富集（推荐；默认 type 需额外包，可用
# lfcShrink(dds, contrast=c("group","SG","NC"), type="apeglm")，需安装 apeglm）
# 本模板直接使用 Wald 统计量 stat 做 GSEA 排序，避免额外依赖。

# ---- 2.1 样本聚类：PCA（vst 方差稳定变换）----
vst_dds <- vst(dds, blind = TRUE)
p_pca <- plotPCA(vst_dds, intgroup = "group") +
  ggtitle("PCA (Treat vs Ctrl)") + theme_classic(base_size = 12)
save_pdf_png(p_pca, file.path(out_dir, "pca"), w = 6, h = 5)

# ---- 2.2 MA 图 ----
res_df <- as.data.frame(res)
p_ma <- ggplot(res_df, aes(baseMean, log2FoldChange)) +
  geom_point(alpha = 0.3, size = 0.6,
             colour = ifelse(res_df$padj < 0.05 & !is.na(res_df$padj),
                            ifelse(res_df$log2FoldChange > 0, "red", "blue"), "grey50")) +
  scale_x_log10() + geom_hline(yintercept = 0) +
  labs(x = "Mean of normalized counts", y = "log2 fold change",
       title = "MA plot (Treat vs Ctrl)") + theme_classic(base_size = 12)
save_pdf_png(p_ma, file.path(out_dir, "ma_plot"), w = 6, h = 5)

# ---- 2.3 火山图 ----
p_vol <- ggplot(res_df, aes(log2FoldChange, -log10(pvalue))) +
  geom_point(alpha = 0.3, size = 0.6,
             colour = ifelse(res_df$padj < 0.05 & !is.na(res_df$padj),
                            ifelse(res_df$log2FoldChange > 1, "red",
                                   ifelse(res_df$log2FoldChange < -1, "blue", "grey50")),
                            "grey50")) +
  geom_vline(xintercept = c(-1, 1), linetype = 2) +
  labs(x = "log2 fold change", y = "-log10 p-value",
       title = "Volcano (Treat vs Ctrl)") + theme_classic(base_size = 12)
save_pdf_png(p_vol, file.path(out_dir, "volcano"), w = 6, h = 5)

# ---------------------------------------------------------------------------
# 3. DEG 筛选与导出（|log2FC|>1 且 padj<0.05）
# ---------------------------------------------------------------------------
res_df$gene <- rownames(res_df)
deg <- res_df[!is.na(res_df$padj) & abs(res_df$log2FoldChange) > 1 &
              res_df$padj < 0.05, ]
up   <- deg[deg$log2FoldChange > 0, ]
down <- deg[deg$log2FoldChange < 0, ]
cat(sprintf("DEG: 上调 %d, 下调 %d（阈值 |log2FC|>1, padj<0.05）\n",
            nrow(up), nrow(down)))

write.csv(deg, file.path(out_dir, "DEG_table.csv"), row.names = FALSE)
write.csv(up,  file.path(out_dir, "DEG_up.csv"),   row.names = FALSE)
write.csv(down,file.path(out_dir, "DEG_down.csv"), row.names = FALSE)

# 全量结果（做 GSEA 排序用）
write.csv(res_df, file.path(out_dir, "all_results.csv"), row.names = FALSE)

# ---------------------------------------------------------------------------
# 4. 富集分析（clusterProfiler；KEGG 需联网，失败自动跳过）
# ---------------------------------------------------------------------------
deg_genes <- deg$gene
if (length(deg_genes) > 0) {
  # ---- 4.1 GO（BP/MF/CC），keyType=SYMBOL 直接使用基因符号 ----
  cat("\n== GO 富集（BP）==\n")
  ego <- enrichGO(deg_genes, OrgDb = org.Hs.eg.db, keyType = "SYMBOL",
                  ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05,
                  qvalueCutoff = 0.2, readable = FALSE)
  if (!is.null(ego) && nrow(as.data.frame(ego)) > 0) {
    write.csv(as.data.frame(ego), file.path(out_dir, "GO_BP.csv"), row.names = FALSE)
    print(head(as.data.frame(ego), 5))
  } else cat("（无显著 GO BP 条目，demo 合成数据属正常）\n")

  # ---- 4.2 KEGG（需联网；离线时 tryCatch 捕获并跳过）----
  cat("\n== KEGG 富集（需联网）==\n")
  id_map <- bitr(deg_genes, fromType = "SYMBOL", toType = "ENTREZID",
                 OrgDb = org.Hs.eg.db)
  kk <- tryCatch(
    enrichKEGG(gene = id_map$ENTREZID, organism = "hsa",
               pvalueCutoff = 0.05, qvalueCutoff = 0.2),
    error = function(e) { cat("KEGG 联网失败，已跳过:", conditionMessage(e), "\n"); NULL }
  )
  if (!is.null(kk) && nrow(as.data.frame(kk)) > 0) {
    write.csv(as.data.frame(kk), file.path(out_dir, "KEGG.csv"), row.names = FALSE)
    print(head(as.data.frame(kk), 5))
  } else cat("（无显著 KEGG 条目或已跳过）\n")

  # ---- 4.3 GSEA（以 Wald stat 排序，gseGO 无需外部基因集文件）----
  cat("\n== GSEA（GO BP，按 stat 排序）==\n")
  ranks <- res_df[!is.na(res_df$stat), c("gene", "stat")]
  ranks <- ranks[order(ranks$stat, decreasing = TRUE), ]
  rvec <- setNames(ranks$stat, ranks$gene)
  gse <- tryCatch(
    gseGO(rvec, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP",
          minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.25, seed = TRUE),
    error = function(e) { cat("GSEA 运行失败，已跳过:", conditionMessage(e), "\n"); NULL }
  )
  if (!is.null(gse) && nrow(as.data.frame(gse)) > 0) {
    write.csv(as.data.frame(gse), file.path(out_dir, "GSEA_BP.csv"), row.names = FALSE)
    print(head(as.data.frame(gse), 5))
  } else cat("（无显著 GSEA 条目，demo 合成数据属正常）\n")
} else {
  cat("无 DEG，跳过富集分析。\n")
}

# ---------------------------------------------------------------------------
# 5. 可复现信息（记录版本，评审/投稿必备）
# ---------------------------------------------------------------------------
sink(file.path(out_dir, "session_info.txt"))
cat("R version / key packages:\n")
print(sessionInfo())
sink()

cat(sprintf("\n完成。输出目录: %s\n", normalizePath(out_dir)))
cat("说明：demo 数据为合成数据，仅验证管线；真实结果以接入真实数据后为准。\n")
