# ============================================================
# 01_data_preprocess.R
# 数据下载与预处理
# ============================================================

library(GEOquery)
library(limma)
library(tidyverse)

# ---- 设置工作目录 ----
# setwd("/path/to/your/project")

# ---- 第1步：下载GEO数据 ----
# 示例数据集
datasets <- list(
  COPD = "GSE38974",
  Diabetes = "GSE25724"
)

# 下载COPD数据
cat("正在下载 COPD 数据集:", datasets$COPD, "...\n")
gse_copd <- getGEO(datasets$COPD, GSEMatrix = TRUE, getGPL = FALSE)
expr_copd_raw <- exprs(gse_copd[[1]])
cat("原始矩阵维度:", dim(expr_copd_raw), "\n")

# 下载糖尿病数据（可选，如果做交集分析）
# cat("正在下载 Diabetes 数据集:", datasets$Diabetes, "...\n")
# gse_dm <- getGEO(datasets$Diabetes, GSEMatrix = TRUE, getGPL = FALSE)
# expr_dm_raw <- exprs(gse_dm[[1]])

# ---- 第2步：提取样本信息 ----
pd_copd <- pData(phenoData(gse_copd[[1]]))
cat("样本信息列名:\n")
print(colnames(pd_copd))

# 根据实际数据集调整分组提取逻辑
# 以下示例针对 GSE38974，其他数据集需要修改
group_copd <- pd_copd$characteristics_ch1.1
group_copd <- ifelse(grepl("COPD", group_copd), "COPD", "Control")
cat("分组情况:\n")
print(table(group_copd))

# ---- 第3步：探针注释 ----
# 获取平台注释
gpl_id <- annotation(gse_copd[[1]])
cat("平台ID:", gpl_id, "\n")

gpl_table <- Table(getGEO(gpl_id))
cat("平台注释列名:\n")
print(colnames(gpl_table))

# 创建探针到基因的映射
# 注意：不同平台列名不同，需要根据实际调整
probe_map <- gpl_table[, c("ID", "Gene Symbol")]
colnames(probe_map) <- c("ProbeID", "GeneSymbol")
probe_map <- probe_map[probe_map$GeneSymbol != "" & !is.na(probe_map$GeneSymbol), ]

# 过滤表达矩阵
expr_copd <- expr_copd_raw[rownames(expr_copd_raw) %in% probe_map$ProbeID, ]
probe_map <- probe_map[match(rownames(expr_copd), probe_map$ProbeID), ]

# ---- 第4步：合并同一基因的多个探针 ----
df_copd <- as.data.frame(expr_copd)
df_copd$GeneSymbol <- probe_map$GeneSymbol

# 去除多基因映射
df_copd <- df_copd[!grepl("///", df_copd$GeneSymbol), ]

# 按基因分组取平均
df_copd_long <- df_copd %>%
  pivot_longer(cols = -GeneSymbol, names_to = "Sample", values_to = "Expression")

expr_copd_gene <- df_copd_long %>%
  group_by(GeneSymbol, Sample) %>%
  summarise(Expression = mean(Expression), .groups = "drop") %>%
  pivot_wider(names_from = Sample, values_from = Expression) %>%
  column_to_rownames("GeneSymbol")

expr_copd_gene <- as.matrix(expr_copd_gene)
cat("注释后基因数:", nrow(expr_copd_gene), "\n")

# ---- 第5步：log2转换 ----
if (max(expr_copd_gene) > 50) {
  expr_copd_gene <- log2(expr_copd_gene + 1)
  cat("已执行log2转换\n")
}

# ---- 第6步：过滤低表达基因 ----
keep <- rowSums(expr_copd_gene > 1) >= ncol(expr_copd_gene) * 0.2
expr_copd_filtered <- expr_copd_gene[keep, ]
cat("过滤后基因数:", nrow(expr_copd_filtered), "\n")

# ---- 第7步：保存 ----
dir.create("data", showWarnings = FALSE)
dir.create("output", showWarnings = FALSE)

write.csv(expr_copd_filtered, "data/GSE38974_COPD_expression_clean.csv")

sample_info <- data.frame(
  Sample = colnames(expr_copd_filtered),
  Group = group_copd[match(colnames(expr_copd_filtered), rownames(pd_copd))]
)
write.csv(sample_info, "data/GSE38974_COPD_sample_info.csv", row.names = FALSE)

cat("预处理完成！最终矩阵维度:", dim(expr_copd_filtered), "\n")