GEO数据下载与处理：三种平台注释方法详解

生信学长
2025年10月20日
生信基础知识
GEO数据下载与处理：三种平台注释方法详解

点击卡片关注，一起学习生信分析！



大家晚上好，最近收到一些同学反馈GEO数据下载的一些问题，晚上重新整理一下GEO数据下载的两种方式和平台注释方法。

1. 背景介绍

GEO（Gene Expression Omnibus）是NCBI维护的基因表达数据库，包含了全球科研人员上传的大量芯片和测序数据。与TCGA不同，GEO数据来源更广泛，包括各种疾病、组织和实验条件的数据，是生信分析的重要数据来源。

从GEO下载的芯片数据通常包含探针ID（Probe ID），需要转换成基因名称（Gene Symbol）才能进行后续分析。不同的芯片平台有不同的注释方法，本文将介绍三种常用的平台注释方法。

「应用场景」 ：

  • 下载公共数据进行验证分析
  • 整合多个数据集进行meta分析
  • 探针ID到基因Symbol的转换
  • 处理不同平台的芯片数据

2. 分析原理

2.1 GEO数据结构

GEO数据主要包含以下几个部分：

  • 「GSE（Series）」 ：一个完整的研究，包含多个样本
  • 「GSM（Sample）」 ：单个样本的数据
  • 「GPL（Platform）」 ：芯片平台信息，包含探针注释

2.2 芯片数据注释方法

根据平台信息的获取方式，主要有三种注释方法：

「方法一：使用txt格式的平台注释文件」

  • 适用于：大多数平台都有对应的txt注释文件
  • 优点：信息全面，包含多种注释
  • 缺点：文件较大，需要手动下载

「方法二：使用SOFT格式的平台文件」

  • 适用于：GEO提供的标准格式
  • 优点：可以通过GEOquery自动下载
  • 缺点：下载速度可能较慢

「方法三：使用Bioconductor注释包」

  • 适用于：常用的标准平台（如Affymetrix）
  • 优点：简单快捷，不需要下载文件
  • 缺点：仅限于有注释包的平台

2.3 数据处理流程

  1. 下载表达矩阵（series_matrix文件）
  2. 判断是否需要log转换
  3. 提取临床信息
  4. 获取平台注释信息
  5. 探针ID转基因Symbol
  6. 处理重复基因（保留最大表达值，也可以平均值）
  7. 导出标准化表达矩阵

3. 环境准备

3.1 R包安装 

# 安装字符串处理包
install.packages("stringr")
install.packages("dplyr")
# 安装Bioconductor管理器
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

# 安装GEOquery包（用于下载GEO数据）
BiocManager::install("GEOquery")
# 安装平台注释包（根据需要选择）
# 示例：Affymetrix Human Genome U133 Plus 2.0 Array
BiocManager::install("hgu133plus2.db")

3.2 数据准备

「从GEO下载数据的方式」 ：

  1. 「自动下载」 （推荐）：使用GEOquery包的 getGEO() 函数
  2. 「手动下载」 ：

    • 访问 GEO数据库
    • 搜索GSE编号（如GSE74777）
    • 下载 series_matrix.txt.gz 文件
    • 下载平台文件（GPL文件）

4. 代码实现

方法一：使用txt格式平台注释文件（示例：GPL17586）

4.1 加载R包 

# 加载必要的R包
library(GEOquery)    # 用于下载和处理GEO数据
library(stringr)     # 字符串处理
library(dplyr)       # 数据处理

4.2 设置工作目录 

# 设置工作目录
# 建议为每个数据集创建独立文件夹
setwd("你的工作目录/GSE74777")

# 检查当前工作目录
getwd()

4.3 下载并读取表达矩阵 

# 下载GEO数据集
# 如果目录下已存在对应文件，会自动读取而不重新下载
gset <- getGEO("GSE74777",
               destdir = ".",           # 保存到当前目录
               AnnotGPL = FALSE,        # 不下载注释信息（我们手动处理）
               getGPL = FALSE)          # 不下载平台文件

cat("数据下载完成！\n")
# 提取表达矩阵
# gset是一个列表，通常第一个元素就是我们需要的数据
dat <- exprs(gset[[1]])

# 查看矩阵维度
dim(dat)
cat("表达矩阵维度:", nrow(dat), "个探针 x", ncol(dat), "个样本\n")

# 查看前几行
dat[1:5, 1:5]

4.4 判断并执行log转换

芯片数据可能已经log转换过，也可能是原始值。我们需要判断并进行适当的转换。 

# 判断是否需要log转换
ex <- dat
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE))

# 判断逻辑：
# 1. 如果99%分位数 > 100，说明是原始值
# 2. 如果最大值-最小值 > 50 且最小值 > 0，说明是原始值
# 3. 如果数据分布在0-2之间，说明已经log转换
LogC <- (qx[5] > 100) ||
        (qx[6] - qx[1] > 50 && qx[2] > 0) ||
        (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
# 执行log2转换
if (LogC) {
  ex[which(ex <= 0)] <- NA    # 将负值和0替换为NA
  dat <- log2(ex)
  cat("log2转换完成\n")
} else {
  cat("数据已经过log转换，无需再次转换\n")
}

4.5 提取临床信息 

# 从GEO对象中提取样本信息（临床数据）
pd <- pData(gset[[1]])

# 查看临床信息的列名
colnames(pd)

# 查看前几行
head(pd[, 1:5])
# 导出临床信息到CSV文件
write.csv(pd, 'clinical_GSE74777.csv', row.names = TRUE)
cat("临床信息已保存到 clinical_GSE74777.csv\n")

4.6 读取平台注释文件（txt格式） 

# 读取GPL平台注释文件
# 这个文件需要从GEO网站手动下载
# 下载地址：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL17586

gpl <- read.table("GPL17586-45144.txt",
                  header = TRUE,
                  fill = TRUE,
                  sep = "\t",
                  comment.char = "#",      # 跳过注释行
                  stringsAsFactors = FALSE,
                  quote = "")

# 查看注释文件的列名
colnames(gpl)

4.7 提取探针ID和基因Symbol 

# 提取探针ID和基因symbol两列
ids <- gpl[, c("ID", "gene_assignment")]

# 修改列名，便于后续处理
colnames(ids) <- c('probe_id', 'symbol')

# 查看前几行
head(ids)
# 从gene_assignment列提取基因symbol
# 该列格式示例："NM_001276352 // LINC02035 // long intergenic non-protein..."
# 我们需要提取"//"之间的第二部分，即基因名称

ids$symbol <- trimws(str_split(ids$symbol, '//', simplify = TRUE)[, 2])

# 查看提取结果
head(ids)

4.8 清理注释数据 

# 去掉没有基因symbol注释的探针
ids <- ids[ids$symbol != '', ]
ids <- ids[ids$symbol != '---', ]
ids <- ids[ids$symbol != '--- ', ]

cat("过滤后保留", nrow(ids), "个有注释的探针\n")
# 将探针ID转换为字符型
# 有些探针ID全是数字，R会自动识别为数值型，需要转换
rownames(dat) <- as.character(rownames(dat))
ids$probe_id <- as.character(ids$probe_id)

4.9 匹配表达数据和注释信息 

# 只保留在表达矩阵中存在的探针
ids <- ids[ids$probe_id %in% rownames(dat), ]

cat("匹配后保留", nrow(ids), "个探针\n")
# 提取匹配的表达数据
dat <- dat[ids$probe_id, ]

# 检查行名是否完全匹配
all_match <- all(rownames(dat) == ids$probe_id)
cat("探针ID匹配检查:", ifelse(all_match, "通过", "失败"), "\n")

4.10 合并注释信息和表达数据 

# 将基因symbol信息添加到表达矩阵前面
dat <- cbind(ids, dat)

# 查看合并后的数据
dat[1:5, 1:7]

4.11 处理重复基因 

# 对于一个基因有多个探针的情况，保留表达量最大的探针
dat <- aggregate(. ~ symbol, data = dat, max)

cat("去重后保留", nrow(dat), "个基因\n")
# 检查第一行是否需要删除
# 有时第一行可能是空值或无效数据
head(dat, 3)

# 如果需要删除，取消下面的注释
# dat <- dat[-1, ]

4.12 格式化并导出数据 

# 将基因symbol设为行名
rownames(dat) <- dat[, 1]

# 删除前两列（probe_id和symbol）
dat <- dat[, -c(1, 2)]

# 查看最终矩阵
dim(dat)
dat[1:5, 1:5]
# 重新整理为导出格式
final_data <- data.frame(ID = rownames(dat), dat)

# 导出为文本文件
write.table(final_data,
            file = "GSE74777.txt",
            sep = "\t",
            quote = FALSE,
            row.names = FALSE)

cat("数据已成功导出到 GSE74777.txt\n")
cat("最终矩阵维度:", nrow(final_data), "个基因 x", ncol(final_data)-1, "个样本\n")

方法二：使用SOFT格式平台文件（示例：GPL25318）

SOFT格式是GEO的标准格式，可以通过GEOquery直接下载。

4.13 下载SOFT格式平台文件 

# 设置工作目录
setwd("你的工作目录/GSE116918")

# 下载数据集
gset <- getGEO("GSE116918",
               destdir = ".",
               AnnotGPL = FALSE,
               getGPL = FALSE)

# 提取表达矩阵
dat <- exprs(gset[[1]])

4.14 判断log转换（同方法一） 

ex <- dat
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE))

LogC <- (qx[5] > 100) ||
        (qx[6] - qx[1] > 50 && qx[2] > 0) ||
        (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)

if (LogC) {
  ex[which(ex <= 0)] <- NA
  dat <- log2(ex)
  cat("log2转换完成\n")
} else {
  cat("数据已log转换\n")
}

4.15 下载并解析SOFT文件 

# 下载SOFT格式的平台文件
# 如果本地已存在，会直接读取
GPL <- getGEO(filename = 'GPL25318_family.soft.gz',
              destdir = ".")

# 从SOFT文件中提取注释表格
gpl <- GPL@dataTable@table

# 查看列名
colnames(gpl)

4.16 提取基因Symbol 

# 提取探针ID和基因symbol
ids <- gpl[, c("ID", "Gene Symbol")]
colnames(ids) <- c('probe_id', 'symbol')

# SOFT文件中基因symbol格式可能是："GENE1 /// GENE2 /// GENE3"
# 我们取第一个基因名
ids$symbol <- str_split(ids$symbol, '///', simplify = TRUE)[, 1]

# 清理空值
ids <- ids[ids$symbol != '', ]
ids <- ids[ids$symbol != '---', ]

4.17 后续步骤（同方法一） 

# 匹配、合并、去重、导出的步骤与方法一相同
rownames(dat) <- as.character(rownames(dat))
ids$probe_id <- as.character(ids$probe_id)

ids <- ids[ids$probe_id %in% rownames(dat), ]
dat <- dat[ids$probe_id, ]

dat <- cbind(ids, dat)
dat <- aggregate(. ~ symbol, data = dat, max)

rownames(dat) <- dat[, 1]
dat <- dat[, -c(1, 2)]

final_data <- data.frame(ID = rownames(dat), dat)
write.table(final_data, file = "GSE116918.txt",
            sep = "\t", quote = FALSE, row.names = FALSE)

cat("GSE116918数据处理完成\n")

方法三：使用Bioconductor注释包（示例：GPL570）

对于常用的Affymetrix等标准平台，可以直接使用Bioconductor的注释包。

4.18 使用注释包获取基因Symbol 

# 设置工作目录
setwd("你的工作目录/GSE30219")

# 下载数据
gset <- getGEO("GSE30219",
               destdir = ".",
               AnnotGPL = FALSE,
               getGPL = FALSE)

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

4.19 执行log转换（同上） 

ex <- dat
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE))

LogC <- (qx[5] > 100) ||
        (qx[6] - qx[1] > 50 && qx[2] > 0) ||
        (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)

if (LogC) {
  ex[which(ex <= 0)] <- NA
  dat <- log2(ex)
  cat("log2转换完成\n")
} else {
  cat("数据已log转换\n")
}

4.20 使用注释包获取基因信息 

# 加载平台注释包
# GPL570对应的是hgu133plus2芯片
library(hgu133plus2.db)
library(dplyr)

# 查看注释包中可用的数据类型
ls("package:hgu133plus2.db")
# 获取探针ID与基因symbol的对应关系
ids <- toTable(hgu133plus2SYMBOL)

# 查看数据结构
head(ids)

# 修改列名
colnames(ids) <- c("probe_id", "symbol")

4.21 后续处理（同方法一） 

# 数据类型转换
rownames(dat) <- as.character(rownames(dat))
ids$probe_id <- as.character(ids$probe_id)

# 匹配
ids <- ids[ids$probe_id %in% rownames(dat), ]
dat <- dat[ids$probe_id, ]

# 检查匹配
all(rownames(dat) == ids$probe_id)

# 合并
dat <- cbind(ids, dat)

# 去重
dat <- aggregate(. ~ symbol, data = dat, max)

# 导出
rownames(dat) <- dat[, 1]
dat <- dat[, -c(1, 2)]

final_data <- data.frame(ID = rownames(dat), dat)
write.table(final_data,
            file = "GSE30219.txt",
            sep = "\t",
            quote = FALSE,
            row.names = FALSE,
            col.names = TRUE)

cat("GSE30219数据处理完成\n")

5. 结果解读

5.1 输出文件说明

每个数据集处理后会生成两个文件：

「表达矩阵文件」 （如GSE74777.txt）：

  • 第一列：基因Symbol
  • 其余列：样本ID
  • 数值：log2转换后的表达量

「临床信息文件」 （如clinical_GSE74777.csv）：

  • 包含样本的详细信息（分组、处理条件等）
  • 用于后续的差异分析和分组比较

5.2 三种方法的比较

方法
适用平台
优点
缺点
txt注释文件
所有平台
信息最全面
需要手动下载，文件大
SOFT文件
大多数平台
可自动下载
下载速度慢
Bioconductor包
标准平台
最简单快捷
仅限有包的平台

5.3 数据质量检查 

# 检查表达矩阵
summary(as.vector(as.matrix(dat)))

# 检查缺失值
sum(is.na(dat))

# 检查基因数量
nrow(dat)

# 检查样本数量
ncol(dat)

5.4 常见平台对应的注释包

平台编号
芯片名称
Bioconductor包
GPL570
Affymetrix Human Genome U133 Plus 2.0
hgu133plus2.db
GPL96
Affymetrix Human Genome U133A
hgu133a.db
GPL571
Affymetrix Human Genome U133A 2.0
hgu133a2.db
GPL6244
Affymetrix Human Gene 1.0 ST
hugene10sttranscriptcluster.db
GPL10558
Illumina HumanHT-12 V4.0
illuminaHumanv4.db

6. 常见问题

6.1 getGEO()下载失败

「问题」 ：网络连接超时，无法下载数据

「解决方法」 ： 

# 方法1：增加超时时间
options(timeout = 300)

# 方法2：手动下载后读取
gset <- getGEO(filename = "GSE74777_series_matrix.txt.gz")

6.2 探针ID与表达矩阵不匹配

「问题」 ：注释文件的探针ID数量与表达矩阵不一致

「原因」 ：

  • 平台文件包含所有可能的探针
  • 表达矩阵只包含实际使用的探针

「解决方法」 ： 

# 取交集
ids <- ids[ids$probe_id %in% rownames(dat), ]
dat <- dat[ids$probe_id, ]

6.3 基因Symbol格式问题

不同平台的symbol格式可能不同：

  • 单个基因： "TP53"
  • 多个基因： "TP53 /// MDM2"
  • 带描述： "NM_001276352 // LINC02035 // long intergenic..."

需要根据实际情况调整提取代码。

6.4 重复基因处理策略 

# 策略1：保留最大值（默认）
dat <- aggregate(. ~ symbol, data = dat, max)

# 策略2：保留平均值
dat <- aggregate(. ~ symbol, data = dat, mean)

# 策略3：保留中位数
dat <- aggregate(. ~ symbol, data = dat, median)

6.5 log转换判断不准确

如果自动判断失败，可以手动检查： 

# 查看数据分布
summary(as.vector(dat))
hist(as.vector(dat[, 1]), breaks = 100)

# 手动log转换
dat <- log2(dat + 1)

如果需要文章代码，联系老师免费获取。

( 请发送文章链接告知老师获取哪个代码)





送你一份科研小礼物，开启科研生涯！

点击即可领取：

》〉戳我，领取一份科研小礼物



上一篇
单细胞RNA测序鉴定MMP11+癌相关成纤维细胞作为血管生成和膀胱癌进展的驱动因素
下一篇
多数据集合并与批次校正：整合TCGA和GEO数据
返回文章列表