GSE16561数据集的文章图表复现，小众的illumina表达量芯片

1. 首先是下载数据

gset <- getGEO('GSE16561', destdir=".", AnnotGPL = F,  getGPL = F) 

下载完成后一看数据，

发现好多负值，应该是数据经过背景矫正，log2转换之后又经过scale的数据（z-score的）。也就是小洁老师上课时说的那种不能直接用来做DEG分析的芯片数据。心中顿时飞过一匹草泥马~

2. 登录GEO查看原始数据

没办法了，只能是亲自去GEO界面查看该数据集，果然，这个芯片从来没有听说过，应该是很小众了。GPL是Illumina 的 beadchip，GPL号还不在生信技能树大神整理的注释包列表中（http://www.bio-info-trainee.com/1399.html），不过有rio包，注释信息到方便解决，但是依然隐隐有种前途堪忧的预感。

不管三七二十一，先把能下下来的写着RAW的东东都下下来解压后检查一遍再说。

然而......让人难过的事情再次发生了....

定睛一看这个raw data的结构，和Jimmy大神的五年前的教程里的数据结构相去甚远（《用lumi包来处理illumina的bead系列表达芯片》 http://www.bio-info-trainee.com/1944.html ）

另一个raw包里的芯片信息并没有提及rawdata是如何处理的，似乎对后续的数据分析毫无帮助......

3. 生命不止，折腾不息，查看soft和作者原文

但是，小白我并不想放弃这个GSE，就想看看这个rawdata是怎么处理的，基于生命不止，折腾不息的精神于是我又把soft文件下了下来，定睛一看：

似乎我们离真相又靠近了一点，起码知道这个rawdata大概经历了什么才变成了一开始看到的带了负值的scale数据。此外我还下载了该数据的SCI原文，里面的信息也和我的理解差不多。

4. 大胆假设，小心求证，不服就干，反正错了又不要钱

不过又有几个问题摆在了我的面前：

• 这个GeneSpring软件处理的背景矫正用的都是哪些探针？
• 我们看到的rawdata是经过背景矫正的么？

因为查不到具体的作为背景矫正用的control探针的信息，作为一个初生牛犊不怕虎的生信小白，我做出了一个大胆的决假设：假定我们之前看到的rawdata是经过背景矫正的，反正也是死马当活马医、错了也不要钱，不服就干......

一顿操作猛如虎，经过quantile normalization和log2转换后**，我终于画出了一张漂亮的boxplot......**

为了证实我的操作是否正确，我把这个图发给了Jimmy大神，在收到Jimmy大神的肯定之后，我欢快的进行了后续的DEG分析，当我满心欢喜的查看分析结果时，悲剧再次发生了，我的差异基因结果如下：

然而原作者的结果如下：

差别还是有些大......怀着忐忑的心情，我再次联系上了Jimmy大神，Jimmy大神尽然爽快的答应帮我看一看这个数据集，简直不敢相信大神肯出手相救，在这里请大神收下我的膝盖！

高效的Jimmy大神很快将他的结果和代码发给我了，大神就是大神，复现的结果和原作者基本相同（其中log的底数是e）：

我仔细的拜读了一遍大神的代码，发现原来我只顾着进行quantile normalization和log2转换了，却忽略了小洁老师课里给我们讲的去除重复探针的操作。导致最终的结果与原作者相差甚远。这里附上Jimmy大神的代码（未经大神同意就po出大神的代码，还希望大神不要见怪）：

b=read.table('GSE16561_RAW.txt',header = T)
b[1:4,1:4]
rownames(b)=b[,1]
b=b[,-1]
b=log(b+1)
boxplot(b,las=2)
dat=b
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat,las=2)
pd=pData(a) 
colnames(dat) 
library(stringr)
group_list= str_split(colnames(dat) ,'_',simplify = T)[,2]
table(group_list)

dat[1:4,1:4] 
a 
# 这个 AnnoProbe 包是jimmy大神为了拯救我们小白亲自开发的！
library(AnnoProbe)
ids=idmap('GPL6883',type = 'soft')
head(ids)
colnames(ids)=c('probe_id','symbol')  
ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in%  rownames(dat),]

dat[1:4,1:4]   
dat=dat[ids$probe_id,] 

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  #保留每个基因ID第一次出现的信息(median最大的那个)


save(dat,group_list,file = 'step1-output.Rdata')

经过大神的代码处理之后的rawdata再拿来跑一遍差异分析的流程，后续结果就令人十分满意了！

5. 小结

1. 小白这次大胆的探索真的是歪打正着，在这里写下手账，希望可以给那些被类似数据困扰的小伙伴们提供一定的参考。
2. 做rawdata处理的时候一定不要忘记去重！不要忘记去重！不要忘记去重！。

6. 致谢

本人临床背景，生信零基础，本着对生信的一腔热情，联系上了Jimmy大神，参加了生信技能树的学习班！特别感谢Jimmy大神的启发式引导、还有小洁老师精彩的R语言课程和时不时给我灌下的毒鸡汤！

因为个人的工作时间的调整问题，后半程郭老师和张老师的linux课程我大多数时候无法参加线上直播，但好在有回放可以看，这里要给两位老师说一声抱歉，因为不能参与到直播互动中来而浪费了你们精心备好的课！

总之，很荣幸能够找到生信技能树，向各位老师讨教学习，这真的是一段很棒的经历。

by：Jack Sparrow， 一个生信零基础的小白学员。

2020/08/01