第3章 R环境与GEOquery包：从零搭建你的挖掘基地

好了，咱们正式开始动手。

这一章，说白了就是帮你把“武器”装好。我见过不少新手，分析做到一半发现包装不上，或者版本冲突，心态直接崩了。所以，咱们先把地基打牢。嗯，别急，一步步来。

3.1 R与RStudio：你的核心工作台

R是计算引擎，RStudio是图形界面。我个人习惯用RStudio，调试代码、查看数据都方便得多。

3.2 安装R和RStudio

这一步很简单，但有几个坑要避开。

1. 先装R，再装RStudio。 顺序反了，RStudio会找不到R。
2. 安装路径不要有中文。 我曾经因为用户名是中文，折腾了半天才搞定。R对中文路径支持很差。
3. 建议安装R 4.2以上版本。 有些新包已经放弃支持旧版了。比如Bioconductor 3.18就要求R >= 4.2。

装完以后，打开RStudio，你会看到四个窗格。左下角是控制台，我们就在这里敲命令。

3.3 安装GEOquery包：从Bioconductor下载

GEOquery是Bioconductor家族的包。怎么装？记住这一行就够了：

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GEOquery")

BiocManager会自动帮你处理依赖。第一次跑可能会有点慢，毕竟要下载不少东西。我建议你在网络好的时候做这件事，别在咖啡馆蹭网装。我试过，断断续续真的会疯。

3.4 getGEO函数：核心中的核心

这个函数，是你和GEO数据库之间的桥梁。它的作用就是根据你给的GEO系列号（比如GSE1009），去服务器把数据拽下来。

3.4.1 基本用法

最简单的调用方式是这样：

library(GEOquery)
gse <- getGEO("GSE1009")

它会根据GEO ID自动判断数据类型。返回的对象是一个 list，里面可能包含一个或多个ExpressionSet。

3.4.2 关键参数详解

光会用还不够，你得知道怎么调。我给你列个表：

参数	作用	我常用的值
destdir	下载文件存放目录。如果不指定，默认存到临时目录，下次再用还得重新下。	destdir="geo_data/"
AnnotGPL	是否自动下载GPL平台注释文件。设为TRUE，会省去后续手动找探针-基因对应关系的麻烦。	AnnotGPL=TRUE
getGPL	是否下载GPL平台的原始数据。一般设为FALSE，除非你需要看平台细节。	getGPL=FALSE
GSEMatrix	是否下载表达矩阵。设为TRUE才有表达量数据。	GSEMatrix=TRUE

来看一个实际项目中的写法：

gse <- getGEO(
  "GSE1009", 
  destdir = "geo_data/",
  AnnotGPL = TRUE,
  getGPL = FALSE,
  GSEMatrix = TRUE
)

3.5 下载表达矩阵与元数据

数据拿到了，怎么把它变成我们能用的格式？这要看 gse 对象的内部结构。

3.5.1 提取表达矩阵

如果GEO系列只包含一个芯片，那么：

gse_list <- getGEO("GSE1009", destdir = "geo_data/")
exprs_data <- exprs(gse_list[[1]])

exprs() 函数会从ExpressionSet中挤出表达矩阵。返回的是一个 data.frame，行是探针ID，列是样本名。

3.5.2 提取样本元数据（临床信息）

这个更重要。差异分析的分组信息就藏在里面。

pdata <- pData(gse_list[[1]])
head(pdata[, 1:5])  # 先看前几列，找分组信息

pData() 返回的是一个data.frame，每一行是一个样本，每一列是一个元数据字段。比如 title 列往往写着样本的描述，“正常”还是“肿瘤”。

3.5.3 提取平台注释（探针-基因映射）

探针ID我们看不懂，得转成基因名。

fData <- fData(gse_list[[1]])
head(fData[, 1:3])

fData() 返回的是平台注释表。里面会有一列是 Gene symbol 或 Gene ID。这就是我们做映射的依据。

3.6 一次完整的下载流程（示例）

好了，我们把刚才讲的串起来，写一个完整的脚本：

# 加载包
library(GEOquery)

# 设置下载目录
data_dir <- "GSE_downloads"
if (!dir.exists(data_dir)) dir.create(data_dir)

# 下载数据
gse_id <- "GSE1009"
gse_list <- getGEO(
  gse_id,
  destdir = data_dir,
  AnnotGPL = TRUE,
  getGPL = FALSE,
  GSEMatrix = TRUE
)

# 提取表达矩阵
expr_mat <- exprs(gse_list[[1]])
write.csv(expr_mat, file = paste0(gse_id, "_exprs.csv"))

# 提取元数据
pheno_data <- pData(gse_list[[1]])
write.csv(pheno_data, file = paste0(gse_id, "_pheno.csv"))

# 提取平台注释
annot_data <- fData(gse_list[[1]])
write.csv(annot_data, file = paste0(gse_id, "_annot.csv"))

message("数据下载完成！")

跑完这段代码，你会得到三个CSV文件。这就是后续分析的数据基础。

3.7 常见问题与避坑

我遇到过不少问题，挑两个最典型的说说。

• 下载失败？ 有时候GEO服务器不稳定。我的做法是：等几分钟再试，或者换个时间段。凌晨下载通常快一些。
• 返回的list长度为0？ 要么是GEO ID写错了，要么是这个系列的数据格式不标准。我建议先上GEO官网搜一下这个ID，确认它的数据类型是“Expression profiling by array”。
• AnnotGPL=TRUE报错？ 有可能是该GPL的注释文件在Bioconductor上不存在。那就只能手动去GEO FTP下载了，不过那是另一个话题。

环境搭好了，数据也拉下来了。下一步就该琢磨怎么清洗和差异分析了。别急，慢慢来。