怀疑了不该怀疑的人
1.起因：8.8k的文件报错内存不足
这个问题有两个学生反馈过，但我又重复不出来这个报错，所以先帮他们搞定了需要的文件，然后准备找个网络好的时间远程看看。
但是呢，昨晚十点半，其中一个学生又要了一次数据。曾老板两次看到这个问题，就问我是怎么回事，是不是我写的函数有问题？
??怀疑我??
geo_download是我封装的函数，装了GEOquery和AnnoProbe，图方便，不用每次复制粘贴20多行代码。
我这么靠谱，我怎么会埋坑。干坏事的肯定是另有其人啊！
随后，曾老板拿了原始GEOquery代码给学生跑，才发现原来是GEOquery的锅。
但此时(休息日的晚上十点半，顶着13℃的严寒)我已经从床上弹起来了，还是干点啥吧，比如干翻GEOquery!
2. 5分钟给出替代方案
以前为什么要用GEOquery读取呢，我们是图方便，文件太乱，而GEOquery两三句代码就搞定了。
图方便才用它，它出问题那就不方便了，所以为什么而还要用它？
然后我就敲啊敲，敲出来以下代码，纯R语言基础嘛，在我们的生信入门和数据挖掘直播课里就有，我讲第7年了啊！
首先是从GSE网页上手动下载表达矩阵文件，放在工作目录下：
lines = readLines('GSE7305_series_matrix.txt.gz')
library(stringr)
pd <- str_remove_all(lines, "^!|\"") %>%
str_split("\t",simplify = T)
k = pd[,1] %>% str_starts("Sample_")
pd = pd[k,]
rownames(pd) = pd[,1] %>% str_remove("Sample_")
pd = pd[,-1]
pd = t(pd)
# exp
exp = read.delim('GSE7305_series_matrix.txt.gz',comment.char = '!',row.names = 1)
exp = as.matrix(exp)
谁懂这一刻的含金量啊，主打一个我命由我不由GEOquery，差点害我背锅是吧？那锅给你砸了。
3.其实人人都可以写函数
小白可能要困惑了，人家写了函数，两三句就可以完成。但我们替代方案的代码是不是有点长？
首先整理数据的代码大有可为，学学自己写代码是好事。
其次，我们也可以写函数啊，超简单。
geo_parser <- function(GSE) {
series_matrix_file = paste0(GSE,'_series_matrix.txt.gz')
lines <- readLines(series_matrix_file)
pd <- stringr::str_remove_all(lines, "^!|\"") %>%
str_split("\t", simplify = TRUE)
k <- pd[,1] %>% str_starts("Sample_")
pd <- pd[k, ]
rownames(pd) <- pd[,1] %>%
str_remove("Sample_")
pd <- t(pd[,-1])
exp <- read.delim(series_matrix_file,
comment.char = "!",
row.names = 1) %>%
as.matrix()
return(list(pd = pd, exp= exp))
}
4.写函数就不用每次大段复制粘贴
eset <- geo_parser('GSE7305')
pd = eset$pd
pd[1:4,1:4]
## title geo_accession status
## [1,] "Endometrium/Ovary-Disease 1" "GSM175766" "Public on Apr 09 2007"
## [2,] "Endometrium/Ovary-Disease 2" "GSM175767" "Public on Apr 09 2007"
## [3,] "Endometrium/Ovary-Disease 3" "GSM175768" "Public on Apr 09 2007"
## [4,] "Endometrium/Ovary-Disease 4" "GSM175769" "Public on Apr 09 2007"
## submission_date
## [1,] "Mar 19 2007"
## [2,] "Mar 19 2007"
## [3,] "Mar 19 2007"
## [4,] "Mar 19 2007"
exp = eset$exp
exp[1:4,1:4]
## GSM175766 GSM175767 GSM175768 GSM175769
## 1007_s_at 263.11400 317.26790 323.64440 305.86987
## 1053_at 86.07838 79.46674 70.95246 84.76953
## 117_at 70.43405 54.92562 46.36516 66.78484
## 121_at 212.35556 198.59961 227.35240 197.35152
这不也是三五句搞定了嘛。
注意：除了芯片以外的其他数据，exp都会是0行，包括转录组在内，其他类型的数据要么没有表达矩阵，要么是在GSE页面最后的补充资料里。
5.来自tinyarray的展望
目前GEOquery的问题只发现出现在mac电脑+最新版GEOquery上面，如果你是win电脑，或者早就已经安装过GEOquery，是不受影响的。
我写的tinyarray包里面的geo_download函数就是完成表达矩阵、临床信息和GPL编号提取的，只不过依赖了GEOquery，现在我们有了替代版，tinyarray可以少一个依赖包了，哼。
因为cran的包投稿审核需要时间，使用geo_download遇到问题的同学需要等待下一个版本😄,且代码不用变。
library(tinyarray)
eset = geo_download('GSE7305')
6.人与人之间的信任呢？
遥想当年，以为gtf文件是权威，怀疑了曾老板的代码。结果抓出问题，反而证实了曾老板的代码没有问题。是我不该怀疑他！
所以曾老板给起了推文标题：怀疑了不该怀疑的人，发出去了:
两年前射出的箭，今天正中自己的眉心，他以为GEOquery是权威，怀疑我的函数有问题，所以还是用这个标题吧，呼应上了啊。
本文结论：你可以永远相信小洁老师。
腾讯云开发者
