一网打尽GEO数据集全目录信息

1 GEO菜单初览

1.1 GDS目录

1.2 GPL目录

1.3 GSE目录

1.4 GSM目录

2 使用R包rvest下载GEO菜单

2.1 以GPL数据的一个页面为例

2.2 批量下载GPL数据目录

2.3 批量下载GSE数据目录

2.4 批量下载GSM数据目录

2.5 保存及更新下载的数据

3 使用R包tidyverse整理GEO菜单

3.1 整理GPL数据目录

3.2 整理GSE数据目录

3.3 整理GSM数据目录

3.4 保存清洁的GEO数据目录

3.5 提取感兴趣的某个GSE的相关数据条目

4 有了GEO菜单数据，能做些什么？

5 sessionInfo

1. GEO菜单初览

GEO数据库存储了海量的功能基因组学数据，来源于全球的研究者们利用各种芯片和测序技术产生的数据。无论是做肿瘤还是非肿瘤的小伙伴，都会觉得这个数据库很香，不是吗。“用别人的数据，发自己的文章”，再或者“用别人的数据，筛自己感兴趣的疾病分子”，会不会做梦都笑出了声。

GEO中有哪些具体的数据，有没有一个菜单？肯定是有的，而且菜单相当厚，给点菜带来了诸多不便。

从GEO网站主页，我们可以看到，GEO数据分为了4个大类，也就是GDS、GSE、GPL、GSM。其中GSM的数据条目最多，超过了350万个。在开始创建我们的菜单前，我们先分别来看看数据条目：

image.png

1.1 GDS目录
GDS的数据量最少，且2016/2/1后就没再更新。GDS目录下包含了GDS编号、GDS名称、物种、GDS对应的GPL平台、GDS来源的GSE数据和该GDS数据的样本数量等信息。

image.png

1.2 GPL目录

GPL目录页面，分页列出了所有GPL的简要信息，包括GPL编号、GPL名称、GPL使用的技术类别、物种、该GPL的注释信息行数、应用该GPL平台的GSM数量、应用该GPL平台的GSE数量、联系方式、发布日期等信息。

注意此处的Page size：20，最大可以设置为500/页；总页数：1043，还有Release date，可以按照日期升序或者降序排列。后续我们使用R爬取数据时可以更改这些设置。

image.png

1.3 GSE目录

GSE目录同上，包含了GSE编号、GSE标题、GSE数据的类型、该GSE的样本数量、以及对应的GDS链接、补充文件的格式、联系方式以及该GSE的发布日期。

image.png

1.4 GSM目录

GSM目录包含了所有GSM数据编号、GSM标题、GSM数据类型、物种、发布日期等信息。

image.png

2. 使用R包rvest下载GEO菜单

因为GDS数据不太常用，而且其在GEO中存储的页面没有固定的链接格式。我们此次下载其他三类数据，即GPL、GSE、GSM数据目录。

2.1 以GPL数据的一个页面为例

2.1.1 加载R包

library(rvest)
library(tidyverse)
library(magrittr)
library(lubridate)
library(R.utils)

2.1.2 获取GPL数据目录url链接

将页面设置为每页显示最大数量条目，500/页；将数据按照时间升序排列，方便后续更新数据时无缝衔接；发现目前共有42页，选择第2页。地址栏链接为：https://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms&sort=date&display=500&page=2

首先下载该页面的数据条目表格，代码如下：

url_test <- "https://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms&sort=date&display=500&page=2"

url_test <- "https://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms&sort=date&display=500&page=2"

test_table <- read_html(url_test) %>%  
# 读取网页链接
  html_table() %>%           
# 解析提取网页中的table
.[[1]]                               
# 提取table

print(test_table[1:6,])

2.2 批量下载GPL数据目录

单个页面的数据可以下载，那批量的话，就是一个for/while循环可以解决的事情了。

首先，我们从GPL数据的第一页，来获取当前GPL数据条目的条目总数、总页数。

url_1st_GPL <- "https://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms&sort=date&display=500&page=1"
total_number_GPL <- read_html(url_1st_GPL) %>% 
 html_node("#count") %>% 
# "#count"来源：Chrome浏览器中在需要下载的数据上右键单击，点击检查，弹出的页面中，找到相应数据，右键点击，选择copy selector
  html_text() %>% 
  str_extract(pattern = 
"\\d+"
) %>% as.numeric()

total_number_GPL

page_number_GPL <- read_html(url_1st_GPL) %>% 
  html_node(
"#page_cnt"
) %>% 
  html_text() %>% 
as.numeric()

page_number_GPL

开始批量下载：

urls_GPL <- str_c(
"https://www.ncbi.nlm.nih.gov/geo/browse/?view=platforms&sort=date&display=500&page="

1:page_number_GPL)

urls_GPL[page_number_GPL]

创建一个list，把每个页面中爬取的data.frame表格作为一个元素保存在list中。GEO中的数据每天都在更新，为了节约大家的时间，我已经把大部分数据保存为list，大家需要做的就是根据这个代码进行数据更新了。

load("GEO_menu_lists_upto_20200503.Rda")   
# 加载的数据中包含已下载的GPL_list、GSE_list、GPL_list，后续的所有操作都在这个的基础上进行
length(GPL_list) 

length(GSE_list)
## [1] 259
length(GSM_list)

在上述数据的基础上分别更新相应的数据，代码如下：

total_number_GPL
page_number_GPL

i <- length(GPL_list)
while (i < (page_number_GPL + 1)) {
  tryCatch(expr = {withTimeout(expr = { GPL_list[[i]] <- read_html(urls_GPL[i]) %>% html_table() %>% .[[1]]
# 首先更新保存的GPL_list中的最后一页数据，然后下载新的数据
print(paste0("Page ", i, ", completed"))        
i = i + 1}, 
timeout = 15)     
# 超过15s下载不完成则重试，可根据网络情况适当调整}, 
error = function(e) print(paste0("For page ", i, ", try again")))
# 下载失败时，提供报错信息，但不停止运行，接着再试
}

2.3 批量下载GSE数据目录

思路及代码同前，批量下载GSE数据目录，代码如下：

url_1st_GSE <- "https://www.ncbi.nlm.nih.gov/geo/browse/?view=series&sort=date&display=500&page=1"
total_number_GSE <- read_html(url_1st_GSE) %>% 
  html_node("#count") %>%   html_text() %>%   str_extract(pattern = "\\d+") %>% as.numeric()
total_number_GSE

page_number_GSE <- read_html(url_1st_GSE) %>% 
  html_node("#page_cnt") %>%   html_text() %>% as.numeric()
page_number_GSE

urls_GSE <- str_c("https://www.ncbi.nlm.nih.gov/geo/browse/?view=series&sort=date&display=500&page="
1:page_number_GSE)
urls_GSE[page_number_GSE]


i <- length(GSE_list)
while (i < (page_number_GSE + 1)) {
  tryCatch(expr = {
      withTimeout(expr = {
        GSE_list[[i]] <- read_html(urls_GSE[i]) %>% html_table() %>% .[[1]]
# 首先更新保存的GSE_list中的最后一页数据，然后下载新的数据

print(paste0("Page ", i, ", completed"))
        i = i + 1}, timeout = 15)     
# 超过15s下载不完成则重试，可根据网络情况适当调整}, 
error = function(e) print(paste0("For page ", i, ", try again")))
# 下载失败时，提供报错信息，但不停止运行，接着再试}

2.4 批量下载GSM数据目录

思路及代码同前，批量下载GSM数据目录，代码如下：

url_1st_GSM <- "https://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&sort=date&display=500&page=1"
total_number_GSM <- read_html(url_1st_GSM) %>% 
  html_node("#count") %>% 
  html_text() %>% 
  str_extract(pattern = "\\d+") %>% 
as.numeric()

total_number_GSM
## [1] 3579586

page_number_GSM <- read_html(url_1st_GSM) %>%
  html_node("#page_cnt") %>% 
  html_text() %>% as.numeric()

page_number_GSM
## [1] 7160

urls_GSM <- str_c("https://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&sort=date&display=500&page=",1:page_number_GSM)

urls_GSM[page_number_GSM]
## [1] "https://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&sort=date&display=500&page=7160"

i <- length(GSM_list)
i
## [1] 7160

while (i < (page_number_GSM + 1)) {

  tryCatch(expr = {
      withTimeout(expr = {
        GSM_list[[i]] <- read_html(urls_GSM[i]) %>% html_table() %>% .[[1]]
# 首先更新保存的GSM_list中的最后一页数据，然后下载新的数据

print(paste0("Page ", i, ", completed"))
        i = i + 1}, timeout = 15)     
# 超过15s下载不完成则重试，可根据网络情况适当调整}, 
error = function(e) print(paste0("For page ", i, ", try again")))
# 下载失败时，提供报错信息，但不停止运行，接着再试
}
## [1] "Page 7160, completed"

2.5 保存及更新下载的数据

保存数据为Rdata格式，下次更新可在此基础上再次运行以上代码即可。

save(GPL_list, GSE_list, GSM_list, file = 
"GEO_menu_lists_upto_20200503.Rda")
# 保存为Rda文件后，全部数据条目的list未超过50MB。

3. 使用R包tidyverse整理GEO菜单

R中最习惯和常用的数据格式当然还是数据框了。下一步就是把刚刚下载的lists整理为数据框。

3.1 整理GPL数据目录

首先整理GPL数据，代码如下：

GPL_all <- bind_rows(GPL_list) %>%   
# 将list合并为数据框
 arrange(desc(Samples)) %>%
# 按照每个GPL所有的GSM样本数降序排列
mutate(`Release date`= mdy(`Release date`),
# 更改日期的显示格式
across(.fns = ~ str_replace_all(.x, "\\s{2,}", " / "))
# 将数据框中的全部超过2个的连续空格替换为 /
) %>%   rename_with(~ (paste0("GPL_", str_replace(.x, " ", "_"))))
# 给数据框列名统一重命名，加前缀GPL，去掉列名中的空格
nrow(GPL_all) == total_number_GPL
## [1] TRUE

查看GPL中样本数的Top10：

print(GPL_all[1:10, c(1, 2, 6, 7)])

3.2 整理GSE数据目录

代码类似，整理GSE数据框：

GSE_all <- GSE_list %>%
  dplyr::bind_rows() %>% 
  mutate(`Release date`
= mdy(`Release date`), across(.fns = ~ str_replace_all(.x, "\\s{2,}", " / "))) %>% 
  rename_with(~ (paste0("GSE_", str_replace(.x, " ", "_"))))
nrow(GSE_all) == total_number_GSE
## [1] TRUE

3.3 整理GSM数据目录

GSM_all <- GSM_list %>%
  purrr::map(~ mutate(.x, across(!is.character, as.character))) %>%

# 将GSM_list中元素不是字符串类型的全部转换为字符串
  dplyr::bind_rows() %>%
# 合并数据框
  dplyr::mutate(`Release date`= mdy(`Release date`),

# 更改日期格式
across(.fns = ~ str_replace_all(.x, "\\s{2,}", "/"))# 将2个以上的连续空格替换为 /) 
%>% dplyr::rename_with(~ paste0("GSM_", str_replace(.x, " ", "_")))

# 重命名colnames
# 该步骤略耗时，内存小的话可能会很卡
nrow(GSM_all) == total_number_GSM
## [1] TRUE

3.4 保存清洁的GEO数据目录

保存已整理好的数据目录：

save(GPL_all, GSE_all, GSM_all, file ="GEO_menu_dataframes_upto_20200503.Rda")

3.5 提取感兴趣的某个GSE的相关数据条目

创建一个函数GSE_extract，提取任意的GSE数据相关的条目：

GSE_extract <- function(GSE) {
  GSE_list <- list("GSE"= NULL,
"GSM"= NULL,"GPL"= NULL)

GSE_list[[1]] <- GSE_all %>% dplyr::filter(GSE_Accession %in% GSE)

  GSE_list[[2]] <-  GSM_all %>% dplyr::filter(str_detect(GSM_all$GSM_Series, pattern = paste0(GSE, "\\b")))
  GSE_list[[3]] <-  GPL_all %>% dplyr::filter(GPL_Accession %in% unique(GSE_list[[2]]$GSM_Platform))
return(GSE_list)
}

# 示例1
GSE24206_list <- GSE_extract("GSE24206") 
GSE24206_list[[1]] 
GSE24206_list[[2]] 
GSE24206_list[[3]]

# 示例2

GSE98131_list <- GSE_extract("GSE98131") 
GSE98131_list[[1]] 
GSE98131_list[[2]] 
GSE98131_list[[3]]

4. 有了GEO菜单数据，能做些什么？

数据在手，各种GEO数据的统计分析，你值得拥有。比如，可以回答这个问题：目前测序技术很流行，那近5年还有没有流行的表达谱芯片平台呢？最常选择的是？

top_array <- GSM_all %>%
select(6,7,10) %>% 
 filter(str_sub(GSM_Release_date, 1, 4) %in% as.character(2016:2020)) %>% 
  distinct( ) %>% 
  mutate(GPL_Accession = GSM_Platform) %>% 
  group_by(GPL_Accession) %>% 
  summarise(n = n()) %>% 
  arrange(desc(n)) %>% 
  left_join(GPL_all, by = "GPL_Accession") %>% 
  filter(GPL_Technology != "high-throughput sequencing") %>% 
  slice(1:10)
top_array[, 1:4]

Affymetrix的GPL570已经屈居第二位了，第一名被Illumina的GPL10558平台占领
如果你现在需要做个表达谱芯片的话，还会纠结选择哪个平台吗？
数据已在手，可以做的分析还有很多……

5. sessionInfo

sessionInfo()

转自跳圈联盟