在很多实验的时候都会遇到不同批次的数据整合的情况，那么今天就给大家介绍一个测序数据的批次数据分析的R包sva。首先我们看下包的安装，以及内置数据的安装：
BiocManager::install('bladderbatch')
BiocManager::install('sva')
接下来我们看个体外概念虚拟变量，是人为设定的用于将分类变量引入回归模型中的方法，也就是说分类变量的类型数量就是虚拟变量的变量数。通常情况下，回归分析，逐步回归，分层回归，Logistic回归，PLS回归等这类影响关系研究的方法时，才可能涉及到虚拟变量设置。其它分析方法并不会涉及。公式如下：
接下来我们看下实例：
library(bladderbatch)
data(bladderdata)
dat <- bladderEset[1:5000,]
pheno = pData(dat)
edata = exprs(dat)
batch = pheno$batch
mod = model.matrix(~as.factor(cancer),data=pheno)
那么如果我们只是单纯的引入所用的变量，而不是针对某个变量的虚拟变量，那么我们可以直接：
代码语言：javascript复制 mod0 = model.matrix(~1,data=pheno)
接下来我们看下其内置的函数：
1. num. sv 评估批次数量。其主要的参数评估的方法设置：be：置换检验，基本思想是：在H0假设成立的前提下，根据研究目的构造一个检验统计量，并利用样本数据，按排列组合的原理，导出检验统计量的理论分布，在实际中往往因为排列组合数太多，而模拟其近似分布，然后求出在该分布中出现观察样本及更极端样本的概率p，通过和0.05比较，做出统计推断。Leek：逼近思想，采用某种函数（线性函数或者非线性函数等等），基于小部分训练样本，来拟合价值函数的真值。实例：
代码语言：javascript复制n.sv =num.sv(edata,mod,method="leek")2. sva 消除除了批次之外的环境或者生物学实验差异。具体参数如下：
其中主要的参数：
Dat 表达数据矩阵
Mod 调整变量的虚拟变量模型；无调整变量时可以用mod0.
Control 设置基因的属性，1-控制变量，0-非控制变量。一般我们得到的矩阵都已经做了这一部分的处理，所以一般都是0
Method 迭代计算的方法，irw 可以自行评估control probes；supervised 对control probes已知矩阵选此项。
实例：
代码语言：javascript复制svobj = sva(edata,mod,mod0,n.sv=n.sv)
其中：
Sv: 每一批对应的样本的替换变量。
Pprob.gam 指的每个基因受到一个或者多个因素影响的概率。
Pprob.b 指的是每个基因受到调整因素影响的概率。
n.sv 指的评估出来的批次数。
3. f.pvalue 主要用来计算批次之间的差异性。同时还自带矫正的函数p.adjust。其使用很简单，我们直接看下实例：
modSv = cbind(mod,svobj$sv)
mod0Sv = cbind(mod0,svobj$sv)
pValuesSv = f.pvalue(edata,modSv,mod0Sv)
qValuesSv =p.adjust(pValuesSv,method="BH")
前面是对无法确定的批次的影响因素的批次处理，那么接下来我们看下对已知的批次数据如何进行合并，需要用到函数ComBat。
其中主要的参数：
mean.omly 一个逻辑值，可以判断是否进行标准化处理，还是只针对批次的均值影响。
Prior.plots 绘制通过核密度估计（黑）——是在概率论中用来估计未知的密度函数，属于非参数检验方法和参数估计（红）的结果比较。
代码语言：javascript复制ref.batch 可以参考的前面num.sv计算的批次数量。par.prior 判断用哪种估计方法，TRUE 参数估计，FALSE 非参数估计实例：# parametric adjustment
combat_edata1 = ComBat(dat=edata,batch=batch, mod=NULL, par.prior=TRUE, prior.plots=FALSE)
# non-parametric adjustment, mean-onlyversion
combat_edata2 = ComBat(dat=edata,batch=batch, mod=NULL, par.prior=FALSE, mean.only=TRUE)
# reference-batch version, with covariates
combat_edata3 = ComBat(dat=edata,batch=batch, mod=mod, par.prior=TRUE, ref.batch=3)
当然如果设置prior.plots=T，那么就会得到第一批的评估的样本分布图。除了核算法的选择外，带宽（bandwidth）也会影响密度估计，过大或过小的带宽值都会影响估计结果。实例
combat_edata1 = ComBat(dat=edata,batch=batch, mod=NULL, par.prior=TRUE,prior.plots=T)
以上就是整个过程里，通过上面的处理就可以将微阵列的测序数据进行多批次的合并，当然，个人觉得数据类型类似的其他实验数据也是可以用此包进行合并的。
