绝大部分数据集只需要使用我们的标准代码，无论是测序还是芯片的多组学，还是单细胞，都很容易全方面分析和出图，只需要极少的代码知识，所以很多小伙伴自嘲自己是调包侠。
而且并不是所有的数据集都是标准的，所以就必然面临各个包的各个函数的参数的选择，所以大家会戏谑自己是调参侠。
这个时候遇到了非常规数据并且处处碰壁的小伙伴就会比较焦灼，认为是自己对生物信息学底层算法的不理解，问我要数学底层知识点书籍，我当然是可以随手给出一大堆，可是so what？难道真的指望一个医学生或者一个生命科学领域浪费了近10年的小伙伴去看公式吗？
虽然是生信底层算法知识很重要，可是也不能说自己分析数据失败就怪底层算法，99%的生信工程师这辈子都不会接触底层算法，没听说过他们都分析数据失败啊！同样是这个小伙伴，焦虑算法，但是问的问题我觉得完全不是算法本身问题：
这个肉眼一看就是批次效应啊， 数据集在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12644
它已经是有两个关联的文献：
Refining molecular pathways leading to calcific aortic valve stenosis by studying gene expression profile of normal and calcified stenotic human aortic valves. Circ Cardiovasc Genet 2009 Oct;2(5):489-98. PMID: 20031625Increased biglycan in aortic valve stenosis leads to the overexpression of phospholipid transfer protein via Toll-like receptor 2. Am J Pathol 2010 Jun;176(6):2638-45. PMID: 20382708其中一个文献说的很清楚，就是两个独立的数据集：
A large-scale quantitative measurements of gene expression was performed on 5 normal and 5 AS valves using Affymetrix GeneChips. A total of 409 and 306 genes were significantly up- and downregulated in AS valves, respectively.To provide a global biological validation of the whole-genome gene expression analysis, the microarray experiment was repeated in a second set of aortic valves with (n=5) or without (n=5) AS.这两个数据集被存放在同一个gse链接里面了，所以是：
GSM317342 hAorticValve_Ctrl1
GSM317343 hAorticValve_Ctrl2
GSM317344 hAorticValve_Ctrl3
GSM317345 hAorticValve_Ctrl4
GSM317346 hAorticValve_Ctrl5
GSM317347 hAorticValve_Case1
GSM317348 hAorticValve_Case2
GSM317349 hAorticValve_Case3
GSM317350 hAorticValve_Case4
GSM317351 hAorticValve_Case5
GSM377368 hAorticValve_Ctrl6
GSM377369 hAorticValve_Ctrl7
GSM377370 hAorticValve_Ctrl8
GSM377371 hAorticValve_Ctrl9
GSM377372 hAorticValve_Ctrl10
GSM377373 hAorticValve_Case6
GSM377374 hAorticValve_Case7
GSM377375 hAorticValve_Case8
GSM377376 hAorticValve_Case9
GSM377377 hAorticValve_Case10
其实不看文章也可以大致猜测出来，很明显上面的1~5的case和control的10个样品，和下面的6~10的case和control的10个样品，应该是不一样的，不然为什么样品排序如此诡异呢？
读取作者提供的cel文件
因为是 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 芯片，所以作者会提供affy芯片的cel文件：
GSE12644_RAW.tar 96.6 Mb TAR (of CEL)
Supplementary file File size
GSM317342.CEL.gz 5.3 Mb
GSM317343.CEL.gz 5.2 Mb
GSM317344.CEL.gz 5.3 Mb
下载这个压缩包，GSE12644_RAW.tar，解压后就可以读取这个文件夹里面的全部的cel文件，代码如下所示：
rm(list = ls())
library(oligo)
library(ggplot2)
celFiles <- list.celfiles(GSE12644_RAW/,listGzipped = T,
full.name=TRUE)
celFiles
cel_data <- oligo::read.celfiles( celFiles )
dim(exprs(cel_data))
exprs(cel_data)[1:4,1:4]
### 标准化(一步完成背景校正、分位数校正标准化和RMA (robust multichip average) 表达整合)cel_data_rma <- oligo::rma(cel_data )
exp_probe <- Biobase::exprs(cel_data_rma)
exp_probe[1:4,1:4]
save(exp_probe,file = exp_probe.Rdata)
或者直接下载作者提供的表达量矩阵
一般来说，affy芯片的cel文件处理有多种算法，mas5或者rma是最经典的， 但是这个选择不需要了解这两个算法的数学底层，既然两个算法都很出名那么他们肯定是效果大差不差，并不会是导致批次效应的， 要不然这么重要的事情学术界怎么可能放任两个算法起冲突给后面的数据分析带来bug呢？
所以可以相信作者，直接读取作者的表达量矩阵，然后走后面的分析：
library(AnnoProbe)
library(GEOquery)
getOption(timeout)
options(timeout=10000)
gse_number <- GSE12644gset <- geoChina(gse_number)
a=gset[[1]]
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数dim(dat)#看一下dat这个矩阵的维度dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列这两个矩阵很容易看看相关性，可以看到批次效应（2个数据集）仍然是大于不同算法（读取cel文件或者下载作者表达量矩阵）带来的差异：
是否需要默认 normalizeBetweenArrays
一般来说，我们拿到了表达量芯片数据的矩阵，会有两个考虑，是否log以及是否可以 normalizeBetweenArrays：
##～～～查看数据是否需要log～～～boxplot(dat[,1:4],las=2)
#dat=log2(dat)boxplot(dat ,las=2)
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat ,las=2)
前面的是否log很简单的，就是看看表达量芯片数据集里面的数据范围是否在100以下集中分布，如果是在成百上千甚至几十万的数值那就log一下即可，后面的normalizeBetweenArrays其实做不做基本上没有什么影响的。
你可以拿normalizeBetweenArrays前后的矩阵看看样品相关性，注释样品的分组，批次，算法等信息。很明显你会看到，normalizeBetweenArrays前后的矩阵虽然boxplot看起来被拉平了，但是后面的PCA和热图，仍然是可以看到肉眼可见的批次效应：
因为 normalizeBetweenArrays 函数本来就不是矫正两个数据集的差异的。
是否需要sva包的ComBat函数
原则上，前面的批次效应肉眼可见，肯定是需要走sva包的ComBat函数，代码如下所示：
load(file = ../GSE12644/step1-output.Rdata)
dat[1:4,1:4]
boxplot(dat,last=2)
# 每个数据集真实的batch是不一样的，需要自己修改batch=rep(c(b1,b2),each=10)
# boxplot(dat)library(sva)
combat_edata1 = ComBat(dat=as.matrix(dat),
batch= batch , mod=NULL,
par.prior=T, prior.plots=F)
combat_edata1[1:4,1:4]
boxplot(combat_edata1)
dat = combat_edata1
可以很明显的看到运行了sva包的ComBat函数之后的表达量矩阵里面确实是没有了两个数据集的差异：
那么问题来了，这个是让我们肉眼心安理得呢，还是真实有帮助呢？
因为sva包的ComBat函数之后的表达量矩阵都可以独立走差异分析，所以就有各自的上下调基因列表，可以韦恩图，但是我们以前介绍过：初学者最喜欢怀疑自己的分析是否正确，比如差异分析的时候就容易陷入上下调基因数量的对比问题，文章可能是上下调一千附近，但是学员自己复现的时候就数量对不上。其实这个问题并不在于上下调基因数量，应该是看质量，这样的对比才有意义。文献差异分析结果是1000个上调基因500个下调基因，但是自己做出来仅仅是50个和25个，其实仅仅是因为使用的筛选阈值不一样。如果画一个差异变化倍数（logFC）散点图，就可以很直观的给出两次分析结果差异了。详见：两次差异分析结果的比较不要局限于韦恩图
我们这个时候并不能很肯定差异基因数量的差异是否真的影响了生物学意义的挖掘，可以直接比较两次分析的通路的一致性。有意思的是，其实变化并不大哦？
那么问题来了，你真的需要去学习生信的底层算法和公式吗？上面的哪个步骤我提到了算法底层， 当然了，每个步骤都有算法的意义，如何看boxplot，什么是批次，如何看pca图，热图，上下调基因如何富集到功能通路。。。
这个数据集里面的20个样品肯定是有差异，但是差异可以区分大小：
第一个差异就是批次，因为是两次芯片每次10个样品，它们之间的差异是非常恐怖的，这个也是芯片早期被大家诟病的地方，因为它测的就不是表达量本身，是一个相对的概念，换一个芯片换一个国家换一个地区换一个时间段每个独立的数据集完全是不可能跟其它数据集对比，所以才有了各式各样的处理批次效应的算法。第二个差异是作者的分组，绝大部分情况下这个差异肯定是很重要也很巨大，比如癌症样品和正常组织，千差万别，癌症样品里面肉眼都能看到很恶心，所以不可能癌症样品跟正常组织混合的。但是，很多时候作者就瞎搞，比如细胞系培养或者处理前后其实作者骗你它根本就没有培养那么你肉眼看数据就会发现很诡异的现象，处理组和对照组几乎是没有区别。。。第三个差异就是流程和算法层面，这个基本上可以忽略在任何成熟的数据处理流程里面，如果你告诉我limma的voom，edgeR和Deseq2结果有差异，你就相当于推翻了过去的十几万篇转录组相关的科学研究，太可怕了，我不需要看任何软件算法的底层原理我都知道这个事情不可能发生学徒作业
检查这个数据集GSE12644对应的两个文献里面的差异基因和通路，看看如何反应在这20个样品里面。
