赞
踩
GEO数据库全称GENE EXPRESSION OMNIBUS,是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年,收录了世界各国研究机构提交的高通量基因表达数据,也就是说只要是目前已经发表的论文,论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。
并且GEO网站这个网站作为各种高通量实验数据的公共存储库。这些数据包括基于单通道和双通道微阵列的实验,检测mRNA,基因组DNA和蛋白质丰度,以及非阵列技术,如基因表达系列分析(SAGE),质谱蛋白质组学数据和高通量测序数据。可以按照文献数据集编号等众多形式进行检索。但是在这篇教程中仅介绍如何从GEO网站上根据数据集编号下载所需要的GEO数据集,并且下载后在R中对数据集进一步处理成后续分析所要的形式。
本项目以非小细胞肺癌GSE37745数据集(肺腺癌)作为展示
选用的数据库是GEO。
实验分组:疾病组,对照组。
我这里使用的R版本是4.2.2
要用到的R包:tidyverse,GEOquery
首先进入GEO网站官网(如下图所示),在检索位置输入数据集编号,点击箭头指向的位置进一步运行搜索。
搜索之后会弹出如下界面:首先需要检查物种类型(Homo sapiens),之后查看数据集的类型是否是高通量测序/芯片数据,我这里是芯片数据(Expression profiling by array),页面往下拉。
如下图所示:包含了使用该数据集的一些文献,以及该数据集对应的注释文件,并且还列出来了数据集中包含的样本。在这里我们比较关注的是注释文件GPL570。
可以点击GPL570查看注释文件信息,页面拉到最后,如下图所示,这个注释文件里包含了芯片的ID及其对应的基因symbol,这也是我们后面需要用到的内容。(并且我们可以看到基因symbol中有的基因后面有三条反斜杠,那是他的别名,后续需要删除)
注意:GEO数据集如果是芯片数据则不需要手动下载,如果是高通量测序的数据则需要手动下载表达矩阵。
获取以上信息后即可直接进入R用代码自动下载
rm(list = ls()) # 删除工作空间中所有的对象
setwd('/XX/XX/XX') # 设置工作路径
if(!dir.exists('./00_rawdata')){
dir.create('./00_rawdata')
} # 判断该工作路径下是否存在名为00_rawdata的文件夹,如果不存在则创建,如果存在则pass
setwd('./00_rawdata/') # 设置路径到刚才新建的00_rawdata下
加载包:
library(GEOquery)
library(tidyverse)
标注一下需要下载的数据集编号,及前面查看的注释文件编号,并且在当前00_rawdata文件夹下创建一个名为GSE37745的文件夹(这是为了方便管理,如果不单独创建文件夹,数据集很多的话,就会显得很乱)
GEO_data <- 'GSE37745'
gene_annotation <- 'GPL570'
if(!dir.exists(paste0(GEO_data))){
dir.create(paste0(GEO_data))
}
setwd(paste0(GEO_data))
获取数据集:
gset <- getGEO(GEO_data , # 前面创建GEO_data对象
destdir = '.',
GSEMatrix = T,
getGPL = F) # getGEO函数自动从官网中获取对应的数据集
expr <- as.data.frame(exprs(gset[[1]])) # 将获取的数据集的表达矩阵提取出来,如果没有表达矩阵,说明是高通量数据需要到网站自己下载
首先来看一下获取到的表达矩阵长啥样:行名为芯片的ID(后面需要比对注释文件改成对应的基因symbol),列名为样本ID。
获取注释文件GPL570:
gpl <- getGEO(gene_annotation, destdir = '.')
gpl <- Table(gpl)
colnames(gpl)
[1] "ID" "GB_ACC"
[3] "SPOT_ID" "Species Scientific Name"
[5] "Annotation Date" "Sequence Type"
[7] "Sequence Source" "Target Description"
[9] "Representative Public ID" "Gene Title"
[11] "Gene Symbol" "ENTREZ_GENE_ID"
[13] "RefSeq Transcript ID" "Gene Ontology Biological Process"
[15] "Gene Ontology Cellular Component" "Gene Ontology Molecular Function"
这里面我们需要用到的就是ID和Gene Symbol,进一步处理ID和基因ID,前面我们也说了注释文件的基因symbol里有的有别称,需要剔除,用到的就是separate函数,可以按照sep = 'XXX’的形式分割。
probe2symbol <- dplyr::select(gpl, 'ID', 'Gene Symbol')
probe2symbol <- filter(probe2symbol, `Gene Symbol` != '')
probe2symbol <- separate(probe2symbol, `Gene Symbol`, into = c('symbol', 'drop'), sep = '///')
probe2symbol <- dplyr::select(probe2symbol, -drop)
names(probe2symbol) <- c('ID', 'symbol')
可以看一下处理后的probe2symbol长啥样,如下图所示,第一列就是芯片ID,第二列就是基因symbol,做到这大家也就可以想象到后面的步骤了,无非就是做个匹配,将表达矩阵中的ID对应成相应的基因symbol。
如下代码所示进一步处理表达矩阵:
dat <- expr
dat$ID <- rownames(dat)
dat$ID <- as.character(dat$ID)
probe2symbol$ID <- as.character(probe2symbol$ID)
dat <- dat %>%
merge(probe2symbol, by='ID')%>% ## 合并注释文件和表达矩阵
dplyr::select(-ID)%>% ## 去除多余信息——芯片ID
dplyr::select(symbol, everything())%>% ## 重新排列
mutate(rowMean = rowMeans(.[grep('GSM', names(.))]))%>% ## 求出平均数
arrange(desc(rowMean))%>% ## 把表达量的平均值从大到小排序
distinct(symbol, .keep_all = T)%>% ## symbol留下第一个
dplyr::select(-rowMean)%>% ## 反向选择去除rowMean这一列
tibble::column_to_rownames(colnames(.)[1]) ## 把第一列变成行名并删除
处理后的dat文件,如下图所示:行名为基因symbol,列为样本名。到了这一步表达矩阵已经处理好了,接下来就是新的问题:样本是怎么分组的?
前面的getGEO函数不仅能获取数据集的表达矩阵,还可以获取杨信息,只需要通过pData函数即可展示样本信息
a <- gset[[1]]
pd <- pData(a) #
如下图所示为pd对象的样本信息:有样本名,类型等等(可自行查看)
这里我们需要用到的是geo_accession和characteristics_ch1.2这两列(前一列就是样本ID,后一列就是分组信息),当然该数据集中有三种类型的疾病,分别是adeno(肺腺癌),large(大细胞),squamous(肺鳞癌)
table(pd$characteristics_ch1.2)
# histology: adeno histology: large histology: squamous
# 106 24 66
group <- data.frame(sample = pd$geo_accession, group = pd$characteristics_ch1.2)
table(group$group)
group <- group[group$group != 'histology: large', ] # 筛选掉large(大细胞)
table(group$group)
group$group <- ifelse(group$group == 'histology: adeno', 'LUAD', 'LUSC')
table(group$group)
# LUAD LUSC
# 106 66
做到这一步基本可以说是大功告成了,目前我们有了基因表达矩阵,还有了样本分组信息,接下来将我们的分组样本信息和表达矩阵匹配一下,保存文件即可。
dat <- dat[, group$sample]
dat <- na.omit(dat)
write.csv(dat, file = paste0('dat.', GEO_data, '.csv'))
write.csv(group, file = paste0('group.', GEO_data, '.csv'), row.names = F)
注:做到这一步,可以说是生信分析已经完成一半了,后续很多的分析都会基于前面处理的这两个文件,因此这最初的第一步也是最重要的一步,只有基础打好,后面分析才会可靠一些。
结语:
以上就是GEO数据下载及处理的所有过程,如果有什么需要补充或不懂的地方,大家可以私聊我或者在下方评论。
如果觉得本教程对你有所帮助,点赞关注不迷路!!!
与教程配套的原始数据+代码+处理好的数据见零基础入门转录组分析——数据处理(GEO数据库——芯片数据)配套资源
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。