如果我想批量下载GEO的表达量矩阵，有没有什么好的办法？-洪萨配资

我们前面一篇文章介绍了 NCBI 下载测序数据的方法一篇果蝇多组学（ChIP-seq、RNA-seq、ATAC-seq）经典文章，模式生物，附代码，因此非常容易复现

那如果我们要直接下载 GEO 表达量数据呢？如果样本不多，到 GEO 网站直接下载即可。

比如这个数据集：GSE1563，我们到 GEO 网站（https://www.ncbi.nlm.nih.gov/geo/）上搜索：

可以看到这个数据集有 62 个样本，其名称以 GSM 开头。

我们要知道 GEO 数据库有这几种数据组织结构：

•平台文件（GPL）：描述芯片设计、探针或可检测元件
•样本文件（GSM）：包含个体实验测量数据
•系列文件（GSE）：将相关样本分组，通常代表完整研究
•数据集文件（GDS）：由 GEO 工作人员整理，代表具有生物学和统计学可比性的数据集

我们先下载系列矩阵文件。点击下方的：Series Matrix File(s)，进入 FTP 下载页面：

下载后解压，可以看到这个文件前面是感叹号（!）开头的注释信息，包括系列的主题、GEO 访问号以及样本的临床信息等。

最后 !series_matrix_table_begin 和 !series_matrix_table_end 之间的是基因的表达量矩阵：行为基因，列为样本。

理想情况下，我们可以从这个文件中解析出样本的注释信息以及基因的表达量矩阵。但有些时候，这些信息却分散在 SOFT 格式文件或 Supplementary file 中，自行解析比较繁琐。

比较方便的做法是使用 GEOquery 这样的工具协助下载数据。

安装 GEOquery

# 1. 安装BiocManager（若未安装） if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") # 2. 安装GEOquery BiocManager::install("GEOquery") # 3. 加载GEOquery library(GEOquery)

下载表达矩阵+样本注释

参数GSEMatrix = TRUE（默认值）会优先下载标准化的GSExxxx_series_matrix.txt.gz文件，该文件已整合表达矩阵与核心样本注释，解析速度快，适合绝大多数基础分析。

gse <- getGEO( GEO = "GSE1563", # GEO数据集编号 destdir = ".", # 数据保存路径（"."代表当前工作目录） GSEMatrix = TRUE, # 优先下载series matrix文件 getGPL = TRUE, # 同时下载芯片平台注释（可选，默认FALSE） AnnotGPL = TRUE # 获取详细的平台注释（可选，默认FALSE） ) # getGEO 返回的是一个列表，因为一个 GSE 可能包含多个平台 # 通常我们取第一个 gse_data <- gse[[1]]

• 下载后，当前目录会生成GSE1563_series_matrix.txt.gz文件（若设置getGPL = TRUE，还会生成平台文件如GPL96_family.soft.gz）;
•数据本地缓存：使用 destdir = "." 参数可以将下载的数据保存在当前目录。下次运行时，getGEO 会优先读取本地文件，避免重复下载。

深入理解 ExpressionSet 对象

数据下载完成，它通常以 ExpressionSet 对象的形式存在。我们需要学会从中提取三类核心信息：

提取表达矩阵 (Expression Matrix)

这是进行定量分析的基础（如基因表达量）。

exp_matrix <- exprs(gse_data) # 查看前几行 head(exp_matrix[, 1:5])

显示：

> head(exp_matrix[, 1:5]) GSM26805 GSM26806 GSM26807 GSM26808 GSM26809 1000_at 954.9 1044.7 1112.5 1123.5 1021.0 1001_at 229.2 192.3 342.1 318.1 601.6 1002_f_at 521.3 365.6 300.1 302.3 321.9 1003_s_at 275.8 198.0 220.5 227.3 252.5 1004_at 718.0 785.1 824.7 745.5 587.3 1005_at 509.6 1323.1 957.5 1157.4 644.0 >

提取样本表型信息 (Metadata/Phenodata)

这里包含样本的分组、年龄、性别等临床指标。

sample_info <- pData(gse_data) sample_info

显示：

> sample_info title geo_accession status submission_date last_update_date type channel_count source_name_ch1 organism_ch1 GSM26805 C1PBL GSM26805 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26806 C2PBL GSM26806 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26807 C3PBL GSM26807 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26808 C4PBL GSM26808 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26809 C5PBL GSM26809 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26810 C6PBL GSM26810 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26811 C7PBL GSM26811 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26812 C8PBL GSM26812 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens GSM26813 AR2PBL GSM26813 Public on Jul 14 2004 Jul 14 2004 Mar 16 2009 RNA 1 PBL Homo sapiens

提取平台注释信息 (Feature Data)

将探针 ID 转换为基因 Symbol。

feature_info <- fData(gse_data) head(feature_info)

显示：

> head(feature_info) ID Gene title Gene symbol Gene ID UniGene title UniGene symbol UniGene ID 1000_at 1000_at mitogen-activated protein kinase 3 MAPK3 5595 1001_at 1001_at tyrosine kinase with immunoglobulin like and EGF like domains 1 TIE1 7075 1002_f_at 1002_f_at cytochrome P450 family 2 subfamily C member 19 CYP2C19 1557 1003_s_at 1003_s_at C-X-C motif chemokine receptor 5 CXCR5 643 1004_at 1004_at C-X-C motif chemokine receptor 5 CXCR5 643 1005_at 1005_at dual specificity phosphatase 1 DUSP1 1843