你要挖的公共数据集作者上传了错误的表达矩阵肿么办（如何让高手心甘情愿的帮你呢？）
事情总不会一直顺风顺水，如果有人带你，学习当然舒服，明天我们会发布学徒招募，敬请期待！
装一点沙子，填你挖的坑
最近收到网友求助，问：
尝试一篇文献的表达差异分析和热图重现，主要参考您Github中GEO-master/GSE42872_main的代码，但我跑出的差异分析列表logFC与文献给出的列表数据不符，尝试了很多次，不清楚是什么原因?
本来我一般是不理会这样的求助的， 毕竟代码都给了，还不会用，总不能怪我了，巧的是我鬼使神差的回复了：
你的问题在哪里，我就没得空去帮你检查，你要是真想我回答，两个办法。
第一个是把你这个文献写一个PPT，介绍这方面背景知识点给我，我学习到了新知识，作为交换，我就帮你修改代码
第二个是，你直接付费我来帮你检查代码
有趣的是，对方马上甩来了一个详细的PPT，让我也学到了知识，所以就投桃报李，帮忙检查代码，结果发现了很有趣的事情，就是这个数据集的作者，居然上传了错误的表达矩阵。
错误的表达矩阵
[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 这个芯片平台怎么可能只有不到五千个探针！
所以肯定是作者上传错了，因为文件大小就不对。
下载CEL文件
这个时候必须要下载原始数据了。
拿到芯片的原始数据，cel文件的下载地址：ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE84nnn/GSE84571/suppl/GSE84571_RAW.tar
下载这个文件很坎坷，因为我们伟大的墙（不过找了还在墙外的学徒帮忙解决了），解压后就处理咯。
处理CEL文件
CEL文件怎么处理呢？
如果你是孤军奋斗，那么恭喜你咯，至少需要四五天的摸索才有可能搞定下面短短3行代码
这个时候我翻的三年前自己整理的代码，参考我的GEO教程，在https://github.com/jmzeng1314/GEO 上面
# BiocManager::install(c( 'oligo' ),ask = F,update = F)
library(oligo)
# BiocManager::install(c( 'pd.hg.u133.plus.2' ),ask = F,update = F)
library(pd.hg.u133.plus.2)
dir='~/Downloads/GSE84571_RAW/'
od=getwd()
setwd(dir)
celFiles <- list.celfiles(listGzipped = T)
celFiles
affyRaw <- read.celfiles( celFiles )
setwd(od)
eset <- rma(affyRaw)
eset
# http://math.usu.edu/jrstevens/stat5570/1.4.Preprocess_4up.pdf
save(eset,celFiles,file = f)
# write.exprs(eset,file="data.txt")
得到的eset这个对象，与我们之前一直讲解的GEOquery包下载是一样的， 所以后续代码不需要变化。
得到表达矩阵和表型信息
a=eset
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
# [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
#dat=log(dat)
dat[1:4,1:4]
boxplot(dat[,1:4])
colnames(exprs)
group_list = paste(substring(celFiles,12,12),
substring(celFiles,24,24),sep = '-')
table(group_list) #统计频率
save(dat,group_list,file = 'step1-output.Rdata')
有了表达矩阵和表型信息后续分析就完全参考我https://github.com/jmzeng1314/GEO代码即可。
差异分析结果跟作者的附件公布的数量相当，基因也基本一致，问题就解决啦。（感兴趣的读者，可以验证一下我所说的！）
全国巡讲约你
第1-11站北上广深杭，西安，郑州， 吉林，武汉，成都，港珠澳（全部结束）
一年一度的生信技能树单细胞线下培训班（已结束）