✦
医学科研新动向
✦
GEO介绍
在生物信息学分析中，公共数据库如GEO（Gene Expression Omnibus）提供了丰富且多样化的基因表达数据，涵盖了从基因到表观遗传的多个层面。GEO数据库不仅包括基因表达数据，还收录了miRNA表达、DNA甲基化、拷贝数变异（CNV）、SNP（单核苷酸多态性）、ChIP-Seq、ATAC-Seq、sc-RNA、RNA-Seq、全基因组关联研究（GWAS）及蛋白质表达等多种类型的数据。这些数据提供了全面的资源，帮助探讨不同条件下的基因表达、调控机制以及与疾病的关联。
----------------------------------
GEO的Entry Type主要分为四类：
Datasets（GDS）：人工整理的关于某个类别的样本集合（多个数据集）。
Series（GSE）：一个实验项目中所有样本的芯片实验数据集合，代表多个样本的实验设计和结果。
Samples（GSM）：单个样本的实验数据。
Platforms（GPL）：描述用于生成数据的芯片平台，例如Affymetrix、Agilent等。
下面整理一个使用R语言从GEO数据库下载和处理数据的完整代码及流程：
/ TO /
/ TO /
01
安装与加载必要的R包
首先，需要确保已经安装了GEOquery和stringr包，这些包将帮助我们从GEO下载数据和进行字符串处理。
rm(list = ls()) # 清除环境options(stringsAsFactors = F)# 安装包install.packages("GEOquery")install.packages("stringr")# 加载包library(GEOquery)library(stringr)# 设置工作路径setwd("C:\\Users\\admin\\Desktop\\01.download")# 获取当前工作路径getwd()
02
获取GEO数据集
在这一步，使用getGEO函数从GEO数据库下载编号为GSE93798的数据集，该数据集包含肾脏样本，分别来自健康对照组和疾病组。该函数从指定的数据集获取基因表达数据，并将其保存到指定目录（通过destdir参数设置）。
AnnotGPL = F：表示不下载平台注释文件，节省下载时间和空间。如果需要详细的基因注释信息，可以将该参数设为TRUE。getGPL = F：表示不下载平台文件（GPL文件）。如果需要探针平台的详细信息，可以设置为TRUE。destdir="."：表示将下载结果保存到当前文件夹。
Sys.setenv("VROOM_CONNECTION_SIZE"=1131072) # 设置环境变量避免读取大文件时报错gset <- getGEO("GSE93798", destdir=".", AnnotGPL = F, getGPL = F) # 下载数据集a <- gset[[1]] # 提取数据集的第一个表达矩阵dat <- exprs(a) # 获取表达矩阵dim(dat) # 查看矩阵的维度
03
检查数据是否需要对数转换
检查数据的范围，判断是否需要进行对数转换。一般来说，如果数值的范围较大，则可能需要进行log变换。如果数值之间差异不大，则不需要对数转换。
# 查看前四行四列的表达数据dat[1:4, 1:4]# 输出结果可能如下所示：# GSM2462533 GSM2462534 GSM2462535 GSM2462536# 100009676_at 4.871449 5.006905 4.914415 4.918524# 10000_at 7.441907 7.712510 7.873540 7.734035# 10001_at 5.788362 5.736114 5.356750 5.889867# 10002_at 3.618703 3.865409 3.861215 3.879499
04
提取样本的临床信息并进行分组
接下来，提取每个样本的临床信息，并根据样本标题中是否包含“control”来进行分组，将样本分为Normal（正常组）和Disease（疾病组）。
# 提取临床信息pd <- pData(a)# 查看临床信息中的前几列head(pd)[, 1:3]# 输出临床信息部分内容如下：# title geo_accession status# GSM2462533 Glomerular compartment from control human kidney H11_F_27_1_Con_L GSM2462533 Public on Jul 03 2017# GSM2462534 Glomerular compartment from control human kidney H12_M_67_1_Con_L GSM2462534 Public on Jul 03 2017# 根据标题中的"control"进行分组group_list <- ifelse(grepl('control', pd$title), 'Normal', 'Disease')table(group_list) # 检查分组结# group_list# Disease Normal# 20 22
05
下载平台探针数据并匹配基因symbol
为了将表达矩阵中的探针ID转换为gene symbol，下载GPL22945平台文件，并将探针ID与gene symbol进行匹配。
# 下载探针平台数据gpl <- getGEO('GPL22945', destdir = ".")id_pre <- gpl@dataTable@table# 提取探针ID和基因符号列ids2 <- id_pre[, c("ID", "Symbol")]# 重命名列名colnames(ids2) <- c("probe_id", "symbol")head(ids2) # 查看前几# 输出如下# probe_id symbol# 1 1000_at CDH2# 2 10000_at AKT3# 3 100009676_at ZBTB11-AS1# 4 10001_at MED6# 5 10004_at NAALADL1# 6 10005_at ACOT8
接下来将探针ID与表达矩阵中的行名进行匹配，并移除重复的gene symbol。
# 将探针ID与表达矩阵行名匹配ids <- as.data.frame(ids2)ids <- ids[match(rownames(dat), ids$probe_id), ] # 匹配探针IDids <- ids[!duplicated(ids$symbol), ] # 删除重复的基因符号# 更新表达矩阵dat <- dat[ids$probe_id, ]rownames(dat) <- ids$symbol # 使用基因符号替代探针ID
06
保存处理后的数据
最后，将处理好的数据保存为R数据文件和文本文件，便于后续分析使用。
save(dat, group_list, file = "GSE93798_download.Rdata")write.table(dat, file = "symbol.txt", sep = "\t", row.names = TRUE, col.names = TRUE, quote = FALSE)
-END-
文字丨本人编写，如有补充，请随时告诉我
