重复不出来我费九牛二虎之力重复不出来的GSEA文章？

Original 生信技能树生信技能树 2022-06-06

前些天我跟大家进行了一个GSEA分析的合理性的讨论，见：做GSEA分析你的基因到底该如何排序

里面的一个教程：费九牛二虎之力也无法重现的GSEA图很多读者反馈说我的代码无法重复，真的是搞笑，难道说真的是重复不出来我费九牛二虎之力重复不出来的GSEA文章？

所以我亲自加他们好友，怼了他们，结果，尴尬了，原来是电脑系统的问题，所以我找了公司的实习生和首席工程师分别来解决粉丝的报错问题。

首先看实习生

使用getGEO下载一个文件的时候遇到一个报错

rm(list = ls()) ## 魔幻操作，一键清空~ options(stringsAsFactors = F) library(GEOquery) # stromal cells from tumor-draining lymph nodes gset <- getGEO('GSE113250', destdir=".", AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件

排除电脑问题

诡异的是这个代码在mac电脑和Linux上都能运行成功，起初怀疑是电脑的问题，换个windows电脑尝试后还是报同样的错.

排除R版本及R包问题

既然不是电脑的问题，那看一下是不是R版本的问题

> sessionInfo() #查看R运行环境 R version 3.5.2 (2018-12-20) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows >= 8 x64 (build 9200) Matrix products: default locale: [1] LC_COLLATE=Chinese (Simplified)_China.936 LC_CTYPE=Chinese (Simplified)_China.936 [3] LC_MONETARY=Chinese (Simplified)_China.936 LC_NUMERIC=C [5] LC_TIME=Chinese (Simplified)_China.936 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] GEOquery_2.50.5 Biobase_2.42.0 BiocGenerics_0.28.0 loaded via a namespace (and not attached): [1] Rcpp_1.0.0 tidyr_0.8.3 crayon_1.3.4 dplyr_0.8.0.1 assertthat_0.2.0 R6_2.4.0 [7] magrittr_1.5 pillar_1.3.1 stringi_1.3.1 rlang_0.3.1 curl_3.3 rstudioapi_0.9.0 [13] limma_3.38.3 xml2_1.2.0 tools_3.5.2 readr_1.3.1 glue_1.3.0 purrr_0.3.1 [19] hms_0.4.2 yaml_2.2.0 compiler_3.5.2 pkgconfig_2.0.2 tidyselect_0.2.5 tibble_2.0.1

R已是最新版本，查看R包版本

> packageVersion("GEOquery") [1] ‘2.50.5’ #也是最新版本，但我决定还是重装一下看看 remove.packages(GEOquery) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") if (!requireNamespace("BiocManager", quietly = TRUE)) + install.packages("BiocManager") library(BiocManager) BiocManager::install("GEOquery")

更新了好些依赖包之后，再次执行下载那个文件，仍然报错。

gset <- getGEO('GSE76275', destdir=".", AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件

使用getGEO下载别的文件（GSE76275）就没有问题说明不是R包版本的问题，问题在这个文件本身。

尝试修改一些参数：

gset <- getGEO('GSE113249', destdir=".", AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件 gset <- getGEO('GSE113249', destdir=".", fill=TRUE, AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件 ######不成功 gset <- getGEO('GSE113249', destdir=".", sep="\t", AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件 #####不成功 gset <- getGEO('GSE113249', destdir=".", #filename =NULL, #GSElimits = NULL, #GSEMatrix = TRUE, #blank.lines.skip=F, #sep="\t", parseCharacteristics = FALSE, #逐个尝试改变这些参数，均不能读取成功 AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件

排除文件matrix问题

查看文件看有什么问题？

在Linux

grep -v '^!' GSE113249_series_matrix.txt |tr "\t" "#"|less -SN

在windows电脑

在notepad++中查看制表符分割

notepad++ 显示全字符,即显示换行符号CRLF (注：这是在windows下)，空格显示的点‘.’ 而UNIX操作系统中以LF为行结束符，并且自动加上那个CR，以便节省一个字节的长度（这样做在当时存储容量不大的计算机上确实是有意义的，而为了不断兼容，到现在这个规定也没有改变）

可以发现文件中的表达矩阵martix分割并没有问题。

注释部分探究

换个方式来读，可以读进来，说明使用getGEO可以把文件下载下来，只是解析文件的时候出错

gset=read.table('GSE113249_series_matrix.txt.gz',sep = '\t',comment.char = '!',header = T)

既然矩阵没有问题，这个文件由两部分，一个是矩阵，一个是注释信息。同时也可以发现使用read.table来读入的时候，省略了文件的注释信息。有理由怀疑这个文件读入不成功可能是注释部分的原因。

gset1=read.table('GSE113249_series_matrix.txt.gz',sep = '\t',comment.char = '!',header = F) gset2=read.table('GSE113249_series_matrix.txt.gz',sep = '\t') gset2=read.table('GSE113249_series_matrix.txt.gz') gset2=read.table('GSE113249_series_matrix.txt.gz',comment.char = '!')

发现报错，而且和下载时报错类似，说明的确是注释文件的问题，接下来就是解决怎么正确读入注释文件就行了。

read.table函数的fill=T 参数能强制把文件注释信息读进来，但列数完全乱了。

gset3=read.table('GSE113249_series_matrix.txt.gz',sep='\t',fill=T) #能读进来，但是读取的列数有误

仔细观察注释信息和两次报错信息

实习生已经非常努力了，可以说是一顿操作猛如虎，可惜没有解决问题，但是他已经吸收了不少我的逻辑思维能力，仍然是值得表扬的，下面我们看看公司首席工程师的出场秀！

首席工程师

可以发现，注释信息有两部分，“!Series”都是两列，“!Sample”都是三列，getGEO报错是19行没有3列，read.table报错是27行没有2列。

为什么读入注释信息报错不一样？

接下来排查getGEO读入时是注释信息哪部分解析出错？

在Linux

把注释部分拆成两个文件

grep ^\!Series_ GSE113249_series_matrix.txt |less -SN grep ^\!Series_ GSE113249_series_matrix.txt >tmp-series.txt grep ^\!Sample_ GSE113249_series_matrix.txt |less -SN grep ^\!Sample_ GSE113249_series_matrix.txt > tmp-sample.txt

在R中读入看是否报错

tmpseries <- read.table('tmp-series.txt',sep='\t') tmpsample <- read.table('tmp-sample.txt',sep='\t')

说明getGEO报错来自解析注释文件“!Sample”部分。

接下来再探究下为什么报错，然后尝试解决

tmpsample <- read.table('tmp-sample.txt',sep='\t',fill=T) #强制读入查看

可以看的第19行只用两列信息

在Linux中查看这部分，有三列信息在notepad++中查看也是三列信息（信息太长，无法截图）

在R中查看

发现19行在R中解析的时候有乱码。

解决乱码

Sys.setlocale("LC_ALL","English") #设置成英文 tmpsample <- read.table('tmp-sample.txt',sep='\t') gset <- getGEO('GSE113250', destdir=".", AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件

查看默认字符编码

修改默认字符编码

读取成功

是不是感觉很神奇，问题就这样迎刃而解，这就是技术，这就是实力！

如果你完全看不懂上面的操作，那你可能需要下面的课程：

生信技能树（爆款入门培训课）全国巡讲约你

生信技能树（爆款入门培训课）巡讲第一站-重庆（已结束）

生物信息学全国巡讲之粤港澳大湾区专场（已结束）

生信技能树（爆款入门培训课）巡讲第二站-济南 （3.16-3.18正在报名）

生信技能树（爆款入门培训课）巡讲-千呼万唤进北京（连续两场正在报名）