【原】你要挖的公共数据集作者上传了错误的表达矩阵肿么办（如何让高手心甘情愿的帮你呢？）

健明 2021-07-14

展开全文

事情总不会一直顺风顺水，如果有人带你，学习当然舒服，明天我们会发布学徒招募，敬请期待！

装一点沙子，填你挖的坑

最近收到网友求助，问：

尝试一篇文献的表达差异分析和热图重现，主要参考您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./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')