正常的illumina芯片数据可以使用lumi包的lumiR.batch函数读取
正常的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热图展现。
写在文末
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。
腾讯云开发者
