正常的illumina芯片数据可以使用lumi包的lumiR.batch函数读取

表达量芯片数据处理，大家应该是非常熟悉了，我们有一个系列推文，在：

它基本上可以应付主流的芯片数据，主要是 affymetrix和illumina以及agilent，当然最简单的就是affymetrix的芯片，但是最近很多小伙伴问illumina芯片数据，主要是因为一些数据产出的作者自己不熟悉，所以 它们并没有按照规则来上传数据，导致大家没办法使用标准代码处理它。

比如数据，在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE58539

可以看到在该页面有两个不同形式的文件，初次接触的小伙伴可能会犹豫下载哪个 ：

File type/resource
GSE58539_Non-normalized_data.txt.gz 4.8 Mb (ftp)(http) TXT
GSE58539_RAW.tar 26.2 Mb (http)(custom) TAR

其实就是 GSE58539_Non-normalized_data.txt.gz 这个 4.8 Mb文件就可以啦，自行下载哦 。

正常的读取该表达量矩阵文件的代码如下所示：

library(GEOquery)
library(limma) 
library(annotate)  
library(lumi) 
studyID='GSE58539' 
fileName <- 'GSE58539_Non-normalized_data.txt' # Not Run
x.lumi <- lumiR.batch(fileName) ##, sampleInfoFile='sampleInfo.txt')
pData(phenoData(x.lumi))   
## Do all the default preprocessing in one step
#lumi.N.Q <- lumiExpresso(x.lumi, normalize = F) 
lumi.N.Q <- lumiExpresso(x.lumi)
### retrieve normalized data
dat <- exprs(lumi.N.Q)
dim(dat)#看一下dat这个矩阵的维度
# GPL13667
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat[ ,1:4] ,las=2) 

save(dat,file = 'dat_from_lumiR.batch.Rdata')

可以看到，这个时候就已经是一个非常正常的芯片表达量矩阵：

             5786057055_A 5786057055_B 5786057055_C 5786057055_D
ILMN_1343291    14.257948    14.120031    14.162906    14.188226
ILMN_1343295    12.500787    12.786823    13.043421    12.972077
ILMN_1651199     6.576717     6.544162     6.593763     6.503210
ILMN_1651209     6.684041     6.733616     6.713588     6.805869

为什么我们不使用标准的geo数据库的gse芯片数据集处理代码呢？让我们看看：

rm(list = ls())   
options(stringsAsFactors = F)  
f='GSE58539_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE58539
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE58539', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE58539_eSet.Rdata')  ## 载入数据
class(gset)  #查看数据类型
length(gset)  #
class(gset[[1]])
gset 

# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
# GPL13667
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat[ ,1:4] ,las=2) 
colMeans(dat[ ,1:4])
save(dat,file = 'dat_from_GEOquery.Rdata')

可以看到，这个时候的表达量矩阵其实是被zscore了，如下所示  ：

> dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
              GSM1413151   GSM1413152 GSM1413153  GSM1413154
ILMN_1343291  0.09844685 -0.043561935 0.00000000  0.02590084
ILMN_1343295 -0.11312294  0.176450730 0.43721294  0.36480618
ILMN_1651199  0.03707123  0.004797936 0.05399847 -0.03390265
ILMN_1651209 -0.02612686  0.022301674 0.00283289  0.09231090

比较一下两个矩阵，代码如下所示：

rm(list = ls())   
options(stringsAsFactors = F)  

library(GEOquery)
library(limma) 
library(annotate)  
library(lumi)

load(file = 'dat_from_GEOquery.Rdata')
dat_from_GEOquery = dat
load(file = 'dat_from_lumiR.batch.Rdata')
dat_from_lumiR.batch = dat

colnames(dat_from_lumiR.batch)
colnames(dat_from_GEOquery)

library(reshape2)
library(ggpubr)
library(patchwork)
df=melt(dat_from_lumiR.batch)
head(df)
ggboxplot(melt(dat_from_lumiR.batch), 
          x = "Var2", y = "value", width = 0.8) +
ggboxplot(melt(dat_from_GEOquery), 
          x = "Var2", y = "value", width = 0.8)
 

出图如下所示：

比较两次表达量矩阵表达量分布

这个时候也可以很容易看出来，如果是标准的geo数据库的gse芯片数据集处理代码拿到表达量矩阵是被zscore的，所以一般来说不建议后续差异分析富集分析等等。但是因为作者给出来了的 GSE58539_Non-normalized_data.txt.gz 这个 4.8 Mb文件，是正常的illumina芯片数据可以使用lumi包的lumiR.batch函数读取后，获得能够进行下游分析的表达量矩阵。

提一个问题：是不是所有的illumina芯片都应该是去下载_Non-normalized_data.txt.gz后缀的文件走我上面给大家的代码呢？

学徒作业

针对这两个表达量矩阵，各自继续后续差异分析富集分析，比较两次后续差异分析富集分析结果的差异。

两次差异分析的结果，以散点图和韦恩图进行展现。

两次富集分析结果，以gsea热图展现。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 [e-mail]

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。