# ============================================================
# 教程配套脚本 03：数据预处理（探针注释、log2 转换、分组）
# 输入：GSE14520 原始表达矩阵 + 表型
# ============================================================

library(GEOquery)
library(tidyverse)

# ---- 1. 读取（若已保存）----
expr <- read.csv("data/raw/GSE14520_expr_raw.csv", row.names = 1, check.names = FALSE)
pheno <- read.csv("data/raw/GSE14520_pheno_raw.csv", row.names = 1, check.names = FALSE)

# ---- 2. 平台注释：探针 ID -> 基因 Symbol ----
# GSE14520 主平台为 GPL3921（Affymetrix HG-U133A 2.0）
gpl <- getGEO("GPL3921")                      # 下载平台注释（首次较慢）
anno <- Table(gpl)
# 注释表中通常有 ID 与 Gene Symbol 列，列名以实际为准：
# table 需包含 "ID" 和 "Gene Symbol"（或 "GeneSymbol"）
probe2gene <- anno[, c("ID", "Gene Symbol")]
colnames(probe2gene) <- c("probe", "symbol")
probe2gene <- probe2gene[probe2gene$symbol != "", ]   # 去掉无注释探针

# ---- 3. 探针映射到基因（同一基因多探针取平均/最大值）----
expr_symbol <- expr %>%
  rownames_to_column("probe") %>%
  inner_join(probe2gene, by = "probe") %>%
  select(-probe) %>%
  group_by(symbol) %>%
  summarise(across(everything(), mean)) %>%   # 多探针取均值
  column_to_rownames("symbol")

dim(expr_symbol)

# ---- 4. log2 转换（若数据未 log 化，检测后转换）----
if (max(expr_symbol, na.rm = TRUE) > 100) {
  expr_symbol <- log2(expr_symbol + 1)
}

# ---- 5. 构建分组信息（示例：以表型中的 tissue 列分组）----
# 实际列名以 pheno 为准，常见为 "characteristics_ch1"（内容形如 "normal liver" / "tumor"）
group <- ifelse(grepl("tumor", pheno$characteristics_ch1, ignore.case = TRUE),
                "Tumor", "Normal")
table(group)

# 保存处理后的表达矩阵与分组
write.csv(expr_symbol, "data/processed/GSE14520_expr_symbol.csv", quote = FALSE)
saveRDS(group, "data/processed/GSE14520_group.rds")

cat("预处理完成！基因数:", nrow(expr_symbol), " 样本数:", ncol(expr_symbol), "\n")
