一网打尽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
- 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()
转自跳圈联盟
