零基础入门转录组分析——数据处理（GEO数据库——芯片数据）

零基础入门转录组分析——数据处理（GEO数据库——芯片数据）

---

---

GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据，也就是说只要是目前已经发表的论文，论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。

并且GEO网站这个网站作为各种高通量实验数据的公共存储库。这些数据包括基于单通道和双通道微阵列的实验，检测mRNA，基因组DNA和蛋白质丰度，以及非阵列技术，如基因表达系列分析（SAGE），质谱蛋白质组学数据和高通量测序数据。可以按照文献数据集编号等众多形式进行检索。但是在这篇教程中仅介绍如何从GEO网站上根据数据集编号下载所需要的GEO数据集，并且下载后在R中对数据集进一步处理成后续分析所要的形式。

---

---

本项目以非小细胞肺癌GSE37745数据集（肺腺癌）作为展示

选用的数据库是GEO。

实验分组：疾病组，对照组。

我这里使用的R版本是4.2.2

要用到的R包：tidyverse，GEOquery

1. 数据集获取

首先进入GEO网站官网（如下图所示），在检索位置输入数据集编号，点击箭头指向的位置进一步运行搜索。

搜索之后会弹出如下界面：首先需要检查物种类型（Homo sapiens），之后查看数据集的类型是否是高通量测序/芯片数据，我这里是芯片数据（Expression profiling by array），页面往下拉。

如下图所示：包含了使用该数据集的一些文献，以及该数据集对应的注释文件，并且还列出来了数据集中包含的样本。在这里我们比较关注的是注释文件GPL570。

可以点击GPL570查看注释文件信息，页面拉到最后，如下图所示，这个注释文件里包含了芯片的ID及其对应的基因symbol，这也是我们后面需要用到的内容。（并且我们可以看到基因symbol中有的基因后面有三条反斜杠，那是他的别名，后续需要删除）

注意：GEO数据集如果是芯片数据则不需要手动下载，如果是高通量测序的数据则需要手动下载表达矩阵。

2. 数据处理（Rstudio）

获取以上信息后即可直接进入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数据库——芯片数据）配套资源

注：配套资源只要改个路径就能运行，本人已检测过可以跑通，请放心食用，食用过程遇到问题，可先自行百度，实在解决不了可以私信

---

• 目录部分跳转链接：零基础入门生信数据分析——导读