【客户福利】本地读取GEO探针表达数据和芯片注释

其实，从报了培训班之后还没做过GEO数据挖掘。大家注意后面代码的可靠性。

常用的GEO数据下载和读取依赖GEOquery::getGEO ，当我们指定选项GEO编号时，其实质还是从要先下载表达矩阵文件和探针注释（getGPL = T）文件到本地，然后再读取。我个人是不太喜欢getGEO这个函数风格的。

对于R语言下载文件的可靠性，我觉得不咋滴，也可能是我写不好代码。

我个人是封了个脚本，专门用于下载GEO数据库文件，该脚本不公开，仅作为客户福利，当客户付费分析金额达到一定额度后才可免费获取。这个脚本的使用方法之前分享过一次，这里再次分享下。

本文以小鼠(Mus musculus)探针数据GSE106172为例演示本地下载GEO文件，并分享下本地读取表达矩阵，提取临床信息，读取和处理芯片数据以及探针ID转换代码。

下载GSE106172数据

脚本特点

• 基础使用，指定GSE106172 编号即可，下载文件的命令会作为标准输出打印，其它提示内容作为标准误打印。
• 脚本参数有2个，第1个参数为必需参数，是以GSE开头的编号；第2个参数默认为"matrix,miniml,soft,suppl" ，即下载全部内容。也可以只指定部分内容，若只指定补充文件部分suppl 。第2个参数一般默认即可。
• 对matrix,miniml,soft,suppl4个部分的每个部分下载尝试最大次数为10次，对于网络不是十分好的时候比较有用。每一步下载成功后会提示Successfully!
• 脚本中文件的下载是指向了数据库中的每个文件，而不是文件夹，因此要对每步进行测试，并抓取所有文件名。
• 当其中某一个步骤尝试了10次都没有下载成功的时候，会打印报错信息，并且在下一行打印再一次运行的信息。这里其实可以更完善下，就是把脚本的第2个参数写出来，但是一般用不到，先这样，回头有需要再改。

ERROR: Maximum attempts reached, unable to download GEO datasets: $gse,$query\nbash $0 $gse

$gse 为GSE编号

$query 为数据库地址的子文件夹链接

基础使用

$ run_NCBI_download_GEO_v0.1.sh GSE106172 > run_download_GSE106172.sh
==>Trying download: GSE106172 ...
Trying download, 第 1 次...
Successfully!
Trying download: GSE106172 matrix ...
Trying download, 第 1 次...
Successfully!
Trying download: GSE106172 miniml ...
Trying download, 第 1 次...
Successfully!
Trying download: GSE106172 soft ...
Trying download, 第 1 次...
Successfully!
Trying download: GSE106172 suppl ...
Trying download, 第 1 次...
Successfully!

下载

• 文件将下载到当前目录的GSE106172 文件夹内。

$ bash run_download_GSE106172.sh
GSE106172_series_matrix.txt.gz                           100% 2008KB  1.5Mb/s    00:02    
Completed: 2008K bytes transferred in 3 seconds
 (4767K bits/sec), in 1 file.
GSE106172_family.xml.tgz                                 100%   23MB 86.7Mb/s    00:04    
Completed: 23844K bytes transferred in 5 seconds
 (35948K bits/sec), in 1 file.
GSE106172_family.soft.gz                                 100%   23MB 88.6Mb/s    00:04    
Completed: 23854K bytes transferred in 5 seconds
 (34826K bits/sec), in 1 file.
GSE106172_RAW.tar                                        100%   83MB 77.3Mb/s    00:11    
Completed: 85830K bytes transferred in 12 seconds
 (56682K bits/sec), in 1 file.
filelist.txt                                             100%  822               --:--    
Completed: 0K bytes transferred in 0 seconds
 (7K bits/sec), in 1 file.

查看下载的文件

• 看文件大小表达矩阵和芯片注释正常。

$ ll GSE106172/*/* |cut -d ' ' -f5-
2.0M Jun 12 20:41 GSE106172/matrix/GSE106172_series_matrix.txt.gz
 24M Jun 12 20:41 GSE106172/miniml/GSE106172_family.xml.tgz
 24M Jun 12 20:42 GSE106172/soft/GSE106172_family.soft.gz
 84M Jun 12 20:42 GSE106172/suppl/GSE106172_RAW.tar
 822 Jun 12 20:42 GSE106172/suppl/filelist.txt

脚本帮助文档

$ bash /home/zheng/scripts/run_NCBI_download_GEO_v0.1.sh
Function: Download GEO datasets according to accession number of GEO.
Usage:	bash run_NCBI_download_GEO_v0.1.sh  GSE206426 
	bash run_NCBI_download_GEO_v0.1.sh  GSE206426 suppl
	bash run_NCBI_download_GEO_v0.1.sh  GSE206426 "matrix,miniml,soft,suppl"

读取并处理GSE106172 数据

• 第一次自己写代码处理GEO数据，也是第一次这么处理GEO数据。

• 表达矩阵和临床信息文件./GSE106172/matrix/GSE106172_series_matrix.txt.gz

• 芯片注释文件./GSE106172/soft/GSE106172_family.soft.gz

• 指定3个参数gse_number,f_matrix 和f_soft

gse_number 这里的指定只是为了后期导出文件时带上编号GSE106172

f_matrix 表达矩阵文件

f_soft 芯片注释文件

• 利用getGEO指定filename 读取表达矩阵，其中设置getGPL=F ，是防止自动联网获取探针注释。
• 正常提取和处理临床信息
• 利用parseGEO 指定fname 读取注释文件。
• 理解stringr::str_extract 提取symbol部分。
• 合并表达数据exp_ori 和注释文件ids，得到表达矩阵exp
• 后续导出表达矩阵，临床信息和芯片注释信息.这里表达矩阵并没有根据SYMBOL进行去重。
• 最后模仿AnnoProbe::idmap 形式，利用assign 和save 报错芯片注释信息。这保存好，慢慢积攒起来就是自己的小数据库。

library(GEOquery)
library(magrittr)
library(dplyr)
library(stringr)
# soft 注释
# matrix 表达矩阵和临床信息

gse_number <- "GSE106172" # 配合输出保存文件
f_matrix <- "./GSE106172/matrix/GSE106172_series_matrix.txt.gz"
f_soft <- "./GSE106172/soft/GSE106172_family.soft.gz"

# 读取表达矩阵
#eSet <- GEOquery::getGEO(filename = f_matrix,getGPL = F)
eSet <- getGEO(filename = f_matrix, GSElimits = NULL,  AnnotGPL = F, 
                getGPL = F) # getGPL 为F，否则会联网获取GPL16570.soft.gz文件

# exp_ori <- eSet@assayData$exprs
exp_ori0 <- exprs(eSet)

# 提取临床信息
pd <- pData(eSet)

# 表达矩阵样本排序
identical(colnames(exp_ori0),rownames(pd))
if( ! identical(rownames(pd),colnames(exp_ori0)) ) {
  exp_ori0 = exp_ori0[,match(rownames(pd),colnames(exp_ori0))] 
}
# 整理表达矩阵
exp_ori <- data.frame(probe_id= rownames(exp_ori0),exp_ori0)
stopifnot(identical(exp_ori$probe_id, exp_ori$probe_id)) # 多一步检查，也可不用

# 读入探针注释
gpl <- parseGEO(fname = f_soft, GSElimits = NULL,AnnotGPL = F, 
                 getGPL = F)
gpl_ori <- gpl@gpls[[1]]@dataTable@table # 获取平台注释表格
# gpl_ori <- gpl@gpls[[1]]@dataTable %>% GEOquery::Table()
# gpl_ori <- gpl@gpls$GPL16570@dataTable %>% GEOquery::Table()
gpl_number <- gpl@header$platform_id # 获取平台编号

if(F){
  txt = data.table::fread(f_soft, sep = "")[[1]]
  txt = txt[txt != ""]
  # return(.parseGPLTxt(txt))
  gpl_ori <- GEOquery:::.parseGPLTxt(txt) %>% GEOquery::Table()
}
# 进一步处理探针注释
gpl_anno <- gpl_ori %>% 
  dplyr::select(ID,gene_assignment) %>% 
  dplyr::filter(gene_assignment != "---") #%>% 
  # tidyr::separate(gene_assignment,c("drop","symbol"),sep="//") %>% 
  # dplyr::select(-drop)

gpl_anno["symbol"] <- stringr::str_extract(gpl_anno$gene_assignment ,pattern = "//.*?//") %>% 
  stringr::str_remove_all(pattern = "//| ")

ids <- gpl_anno %>% dplyr::select("ID","symbol")
colnames(ids) <- c("probe_id","symbol")

exp <- merge(x = ids, y = exp_ori, by = "probe_id") %>% dplyr::select(-"probe_id")

write.table(exp, file = paste0("./",gse_number,"_exp.xls"),
            sep = "\t",na = "",row.names = F, col.names = T, quote = F)
write.table(ids, file = paste0("./",gpl_number,"_anno_probeid2symbol.xls"),
            sep = "\t",na = "",row.names = F, col.names = T, quote = F)
write.table(gpl_ori, file = paste0("./",gpl_number,"_anno_origin.xls"),
            sep = "\t",na = "",row.names = F, col.names = T, quote = F)
write.table(tibble::rownames_to_column(pd,var = "SampleID"), file = paste0("./",gpl_number,"_clin_infor.xls"),
            sep = "\t",na = "",row.names = F, col.names = T, quote = F)

# 保存Rdata #模仿AnnoProbe::idmap存储形式
assign(x = paste0(gpl_number,"_bioc"), ids)
save(list = paste0(gpl_number,"_bioc"), file = paste0("./",gpl_number,"_bioc.rda"))

结果文件

$ ll G*[rx][dl][as] |cut -d ' ' -f5-
3.3M Jun 13 00:02 GSE106172_exp.xls
492K Jun 13 00:02 GPL16570_anno_probeid2symbol.xls
239M Jun 13 00:02 GPL16570_anno_origin.xls
9.3K Jun 13 00:02 GPL16570_clin_infor.xls
202K Jun 13 00:02 GPL16570_bioc.rda