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数据分为了4个大类，也就是GDS、GSE、GPL、GSM。此中GSM的数据条目最多，超过了350万个。在开始创建我们的菜单前，我们先分别来看看数据条目：
1.1 GDS目次
GDS的数据量最少，且2016/2/1后就没再更新。GDS目次下包罗了GDS编号、GDS名称、物种、GDS对应的GPL平台、GDS泉源的GSE数据和该GDS数据的样本数目等信息。
1.2 GPL目次
GPL目次页面，分页列出了所有GPL的扼要信息，包罗GPL编号、GPL名称、GPL利用的技能种别、物种、该GPL的表明信息行数、应用该GPL平台的GSM数目、应用该GPL平台的GSE数目、接洽方式、发布日期等信息。
注意此处的Page size：20，最大可以设置为500/页；总页数：1043，尚有Release date，可以按照日期升序大概降序排列。后续我们利用R爬取数据时可以更改这些设置。
1.3 GSE目次
GSE目次同上，包罗了GSE编号、GSE标题、GSE数据的范例、该GSE的样本数目、以及对应的GDS链接、增补文件的格式、接洽方式以及该GSE的发布日期。
1.4 GSM目次
GSM目次包罗了所有GSM数据编号、GSM标题、GSM数据范例、物种、发布日期等信息。
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]] # 提取tableprint(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_GPLpage_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] 259length(GSM_list)在上述数据的底子上分别更新相应的数据，代码如下：
total_number_GPLpage_number_GPLi <- length(GPL_list)while (i < (page_number_GPL + 1)) { tryCatch(expr = {withTimeout(expr = { GPL_list[] <- read_html(urls_GPL) %>% html_table() %>% .[[1]]# 起首更新生存的GPL_list中的末了一页数据，然后下载新的数据print(paste0(" age ", 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_GSEpage_number_GSE <- read_html(url_1st_GSE) %>% html_node("#page_cnt") %>% html_text() %>% as.numeric()page_number_GSEurls_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[] <- read_html(urls_GSE) %>% html_table() %>% .[[1]]# 起首更新生存的GSE_list中的末了一页数据，然后下载新的数据print(paste0("age ", 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] 3579586page_number_GSM <- read_html(url_1st_GSM) %>% html_node("#page_cnt") %>% html_text() %>% as.numeric()page_number_GSM## [1] 7160urls_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] 7160while (i < (page_number_GSM + 1)) { tryCatch(expr = { withTimeout(expr = { GSM_list[] <- read_html(urls_GSM) %>% html_table() %>% .[[1]]# 起首更新生存的GSM_list中的末了一页数据，然后下载新的数据print(paste0("age ", i, ", completed")) i = i + 1}, timeout = 15) # 超过15s下载不完成则重试，可根据网络环境恰当调解}, error = function(e) print(paste0("For page ", i, ", try again")))# 下载失败时，提供报错信息，但不制止运行，接着再试}## [1] "age 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] TRUE3.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] TRUE3.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)}# 示例1GSE24206_list <- GSE_extract("GSE24206") GSE24206_list[[1]] GSE24206_list[[2]] GSE24206_list[[3]]# 示例2GSE98131_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()转自跳圈同盟 |
