GEO芯片数据下载与处理：从GSE65682看R语言实战流程

核心洞察

先说结论：这篇教程把GEO下载数据到整理成分析格式的流程讲得很清楚。最值得看的是那个探针ID转基因符号的步骤，很多人会在这卡住。整理临床信息和表达矩阵匹配那部分也实用，后续做差异分析、生存分析基本都能直接套用这套模板。

---

找到感兴趣的GEO数据集后，怎么从GEO网站上根据编号把数据下载下来？下载到本地后，又怎么在R里处理成后续分析能用的格式？这里用数据集GSE65682举例，把完整的R脚本操作流程走一遍。

1. 获取数据集

先进GEO官网。在搜索框里输入数据集编号，点旁边的搜索按钮跑一下。

搜索完会跳出来一个页面。重点看两样东西：物种类型（是不是人类样本），还有数据集类型。这次用的数据集是芯片数据，类型标记为“Expression profiling by array芯片表达谱分析，一种测量基因表达水平的技术。”。

接着往下翻页面。能看到这个数据集对应的注释文件，还有数据集里包含的所有样本。注释文件GPL13667是后面要反复用到的。

点进GPL13667查看注释文件信息。把页面滚动到“Data table header descriptions”这块，快速扫一眼这个注释文件里都有啥。这版注释文件包含了芯片的探针ID和它对应的基因符号（gene symbol基因的官方符号（如TP53），用于标识基因名称。），后面处理数据时会用到。

注意一点：有些GEO数据集只给了探针ID和对应的ENTREZID。那种情况就得先把探针ID转成ENTREZID，再把ENTREZID转成gene symbol。

2. 下载数据并处理

拿到上面的信息后就可以进R了，用代码自动下载数据。

先装好需要的R包：

library(BiocManager)
install("GEOquery")

加载这些R包：

library(GEOquery)
library(limma)
library(affy)
library(data.table)
library(dplyr)

连接GEO，在线下载数据集和注释文件。探针ID转symbol的工作后面做：

gset <- getGEO('GSE65682', destdir=".",
               AnnotGPL = TRUE,
               getGPL = TRUE,
               GSEMatrix = TRUE)

提取表达矩阵：

exp <- exprs(gset[[1]])

这时拿到的表达矩阵里，行名是芯片的探针ID，列名是样本ID。

提取样本的临床信息和注释文件：

cli <- pData(gset[[1]])
GPL <- fData(gset[[1]])

从平台信息里提取探针ID和基因符号两列：

gpl <- GPL[, c("ID", "Gene Symbol")]

清洗基因符号这一列。有的探针对应多个基因符号，用"/// "分隔的，只取第一个：

gpl$"Gene Symbol" <- data.frame(sapply(gpl$"Gene Symbol", function(x) unlist(strsplit(x, "/// "))[1]), stringsAsFactors = F)[, 1]

把基因符号前后的空格去掉：

gpl$"Gene Symbol" <- trimws(gpl$"Gene Symbol")

表达矩阵转成数据框格式，加一行探针ID列，然后跟平台注释信息合并：

exp <- as.data.frame(exp)
exp$ID <- rownames(exp)
exp_symbol <- merge(exp, gpl, by = "ID")

移除包含NA值的行。检查基因符号有没有重复的，对重复的基因符号取平均值来去重。最后移除基因符号为"---"的行：

exp_symbol <- na.omit(exp_symbol)
table(duplicated(exp_symbol$"Gene Symbol"))
exp_unique <- avereps(exp_symbol[, -c(1, ncol(exp_symbol))], ID = exp_symbol$"Gene Symbol")
exp_unique <- exp_unique[row.names(exp_unique) != "---", ]
write.csv(exp_unique,"GSE65682_exp_unique.csv")

从临床信息里提取样本ID、28天死亡事件和生存时间这几列。根据28天死亡事件信息创建分组变量：没事件的标为健康组，事件是1的标为死亡，其他标为存活：

group_info <- as.data.frame(cli[, c(1, 52, 55)])
group_info <- group_info %>%
  mutate(group = ifelse(`mortality_event_28days:ch1` == "NA", "Healthy",
                        ifelse(`mortality_event_28days:ch1` == "1", "Dead","Alive")))
group_info <- group_info %>% rename("sample" = "...1",
                                    "status" = "mortality_event_28days:ch1",
                                    "time" = "time_to_event_28days:ch1")
write.csv(group_info,"GSE65682_group_info.csv",row.names = FALSE)

读取之前保存的表达矩阵数据和分组信息。把表达矩阵的行名设成第一列的基因名。提取跟分组信息匹配的样本列：

expr_data <- read.csv("GSE65682_exp_unique.csv")
group_info <- read.csv("GSE65682_group_info.csv")
sample_ids <- unique(group_info$sample)
rownames(expr_data) <- expr_data$X
expr_data_subset <- expr_data[, colnames(expr_data) %in% sample_ids]

从分组信息里筛出非健康样本（就是败血症样本）。提取疾病组样本的表达数据：

group_sepsis <- group_info %>%
  filter(group != "Healthy")
sample_id <- unique(group_sepsis$sample)
exp_sepsis <- expr_data_subset[, colnames(expr_data_subset) %in% sample_id]

创建只包含样本和分组信息的简化数据框。转置表达矩阵，让样本变成行、基因变成列：

group <- group_info[c("sample", "group")]
sample_exp <- t(expr_data_subset)
sample_exp <- as.data.frame(sample_exp)
sample_exp <- tibble::rownames_to_column(sample_exp, var = "sample")
joined_df <- group %>% right_join(sample_exp, by = "sample")
write.csv(joined_df,"GSE65682_exp_group.csv",row.names = F)

拿到处理好的表达矩阵和临床信息后，就可以跑后续的个性化分析了。差异分析、生存分析、风险模型构建这些下游分析都基于这第一步的工作。所以这最初的一步确实是最关键的。

希望这个流程对你有帮助。从GEO数据库下载转录组芯片数据到处理完数据，上面已经走完了完整流程。

核心结论

1. 教程以数据集GSE65682为例（芯片类型为"Expression profiling by array"），其注释文件GPL13667中包含探针ID与对应的基因符号（Gene Symbol）；若数据集仅提供ENTREZID，则需先转为gene symbol。
2. 处理探针ID转基因符号时，对于多个基因符号以"/// "分隔的，只取第一个；去除前后空格；移除NA值及基因符号为"---"的行；对重复的基因符号取平均值进行去重。
3. 临床信息提取自pData中的第1、52、55列，分别对应样本ID、28天死亡事件（mortality_event_28days:ch1）和生存时间（time_to_event_28days:ch1），并根据事件值（NA/1/其他）创建分组为Healthy、Dead、Alive。
4. 最终输出三个CSV文件：GSE65682_exp_unique.csv（探针转基因符号后的表达矩阵）、GSE65682_group_info.csv（样本分组信息）、GSE65682_exp_group.csv（合并分组与转置后的表达数据，样本为行、基因为列）。

常见问题（FAQ）

探针ID转基因符号时遇到多个基因符号怎么办？

使用sapply函数按"/// "分割，只取第一个基因符号，再用trimws去除前后空格。

有些GEO数据集只给了ENTREZID，如何转成gene symbol？

先利用适当包（如biomaRt）将ENTREZID转换为gene symbol，或通过对应关系映射。

处理后的表达矩阵中有很多重复基因符号怎么办？

使用averepslimma包中的函数，用于对重复基因表达值取平均值，实现去重。函数对重复基因符号取平均值去重，并移除值为"---"的行。