批次效应表示样品在不同的批次处理和测量时引入的与生物状态不相关的系统性的技术偏差。很多因素都可能导致批次效应的产生，如不同实验条件、不同操作者、不同公司的试剂、不同批的试剂、实验开展的时间、检测设备、不同的测序批次等。
2014年生信领域的大牛Michael P Snyder在PNAS上发表了一篇文章Comparison of the transcriptional landscapes between human and mouse tissues，比较了人和小鼠不同组织和器官中表达谱的异同。研究发现不同物种之间组织特异表达的基因是一致的，但很多基因在同一物种不同组织的表达相似度大于它们在不同物种同一组织的表达相似度。“我”来引申下 （原文并没有这么直接说），大体可以理解为小鼠的脑与小鼠的肾脏的相似性大于小鼠的脑与人的脑的相似性。“我”得出的这个结论是有一些颠覆认知的，如果这样，用小鼠做为模式动物是否会对人的研究给出相似性的推导？
这篇PNAS文章发出后，芝加哥大学的Yoav Gilad在F1000上发表了一篇文章A reanalysis of mouse ENCODE comparative gene expression data来讨论这个不同于以往认知的研究项目的设计和分析的合理性。
首先作者从FASTQ数据的序列名字的ID中提取出对应测序数据来源的测序仪设备ID和测序通道信息，发现所有数据来源于5个批次，如下图所示，只有最后一个批次同时包含了人和小鼠的器官，其它批次都只包含了人的器官或小鼠的器官。
重现者Yoav Gilad等通过对数据进行重分析，重现了类似于原文中的结果。不论是PCA还是Heatmap的结果，都展示出来源于同一物种的组织或器官倾向于聚类到一起。
重现者Yoav Gilad等采用ComBat移除批次带来的影响，再次绘制PCA和Heatmap，结果显示表达谱按组织类型而非物种聚在了一起。
通常我们在整合多套数据集进行展示时也会加上数据来源信息以展示自己的分析结果未受批次等因素影响。如下图每一列是一个样品，每一行是一个菌群；列注释中有一行为Dataset指示样品来源于 2 个数据集，并且聚类结果没有明显受到数据集来源的影响（四个大的聚类分支中样品来源分布没有明显偏好性）；
通过样本整体表达分布查看有无批次影响。不同来源的样本一般是各自进行标准化（尤其是芯片数据），合并在一起后，可以简单的从整体表达分布来查看是否存在明显的偏移。如下左图存在明显的偏移，则提示有批次效应的存在。校正后，如右图，看上去样本的整体表达分布均一了。但是否批次影响就被移除了，却很难据此下结论。
通过部分基因集的表达变化查看有无批次效应影响。不同来源的数据一起标准化之后，如果标准化效果好的话，样品整体表达分布也会是均一的（如下面左数第二幅图）。但从中随机抽取数百基因却发现其表达收到了批次的影响（如下面左数第三幅图，只展示了数个基因），而且聚类结果也把两组正常样品分到了各自来源相对应的分支上。
合理的实验设计和一致的实验操作是避免批次效应的最好方式。如下图所示如何通过合理的设计来避免检测批次带来的影响。
上面也只是从检测方式的角度设计出一个合理的试验模式降低批次效应的影响，但除了检测方式，还会有很多不可控的因素也会影响到批次，如不同操作人、不同操作时间等客观因素，还有如配对设计实验中不同的个体自身也是批次的因素。
所以需要有个方式去检测和尽量降低批次效应带来的数据偏差的影响。
下面我们以一个具体例子实战（配对样品处理前后基因表达的变化）和检验下效果。为了演示批次效应的影响，大部分代码做了封装，我们只关心核心的地方。如果自己对封装的代码赶兴趣，可以自行查看函数源码。
首先加载所有的包
输入文件2：实验设计信息表 (sampleFile):conditions为处理条件（untrt是对照, trt是加药处理 ），individual标记样品的个体来源 (4个个体：N61311、N052611、N080611、N061011)。
初始化，定义输入、输出和参数
数据读入和标准化
[1] “Read in 32799 genes”[1] “23936 genes remained after filtering of genes with all counts less than 4 in all samples.”[1] “Perform DESeq on given datasets.”estimating size factorsestimating dispersionsgene-wise dispersion estimatesmean-dispersion relationshipfinal dispersion estimatesfitting model and testing[1] “Output normalized counts”[1] “Output rlog transformed normalized counts”
检查数据标准化效果: 标准化后基因在不同样品的表达分布越均一越好。从下图看不出存在批次效应的影响。
初始化，定义输入、输出和参数 （注意covariate变量使用individual列作为了批次信息）
数据读入和标准化，并检查数据标准化效果: 标准化后基因在不同样品的表达分布越均一越好 （此图略过，与上面的表达分布图一致）。
样本聚类查看样品相似性，trt组和untrt组区分明显 (此部分结果略过，与上面的聚类结果一致)
主成分分析PCA查看样品相似性，发现在PC1轴上，样品按处理条件区分开；在PC2轴上，样品按个体区分开，表明不同的个体可能会对后续的差异基因分析造成影响。这个结果也与我们前面不考虑批次因素的结果是一样的。
是不是批次变量加错了呢，还是添加的批次变量未生效？可以说都不是，操作没问题，只是DESeq2处理时只在差异分析模型中考虑批次效应信息，而不会直接校正表达矩阵。那我们先看下加了批次后差异分析的结果怎样，后续我们再讲如何校正表达矩阵。
校正后，差异基因数目变多了，上调多了99个，下调多了61个。不过数目变化，也说明不了太多问题。
拷贝文件数据到网站数据输入处：
从Venn图可以看出，批次校正后既有新增的差异基因，又丢失了之前的一部分差异基因，那么哪个方式更合理呢？
我们再选2个批次校正前鉴定为有差异、批次校正后鉴定为无差异的基因观察下其表达模式。这两个基因的表达模式没看出存在个体本底的一致变化差异。处理前后在不同个体中变化幅度不一，可能是被动变化。但这些基因一定是没有差异吗？我个人也下不出结论，后续得结合其功能再做判断了。
DESeq2，edgeR和limma在考虑批次因素鉴定差异基因时基本操作是一致的，上面我们也完成和比较了已知批次的数据的差异基因鉴定。
后续还有2个问题：
前面文章讲述了批次信息已知时，在差异基因分析中考虑批次效应的影响可以移除部分基因在个体中不同本底表达水平差异的影响，获得的差异基因倍数方差会变小，可以检测出更多差异基因，理论上也是更可靠的方式。那么如果批次信息未知或记录不完善时怎么处理呢？
这里我们就用到了另一个 R 包sva帮助从数据中预测可能存在的混杂因素包括但不限于批次效应的影响。下面我们实际看下这个包鉴定出的混杂因素与批次效应变量之间是否存在关联？利用预测出的混杂因素作为批次信息校正后结果会有什么变化？
首先加载所有的包
输入文件2：实验设计信息表 (sampleFile):conditions为处理条件（untrt是对照, trt是加药处理 ），individual标记样品的个体来源 (4个个体：N61311、N052611、N080611、N061011)。
初始化，定义输入、输出和参数
数据读入和标准化
[1] “Read in 32799 genes”[1] “23936 genes remained after filtering of genes with all counts less than 4 in all samples.”[1] “Perform DESeq on given datasets.”estimating size factorsestimating dispersionsgene-wise dispersion estimatesmean-dispersion relationshipfinal dispersion estimatesfitting model and testing[1] “Output normalized counts”[1] “Output rlog transformed normalized counts”
下面的方式也可以 (svaseq 是在 sva 的基础上对数据做了一个 log转换；如果处理的是芯片数据，通常已经做过 log 换，直接使用 sva 即可)。
添加预测出的Surrogate variable属性到dds对象
可视化展示预测出的Surrogate variable属性与已知的批次信息的关系
从下图可以看出，预测出的混杂因素SV1, SV2与样品来源的个体信息 (individual)还是比较一致的。
根据已知批次信息校正后差异基因数目变多了，上调了99个，下调多了61个。根据预测的混杂因素校正后，上调多了52个，下调少了1个。(这个是只保留2个混杂因素的结果)
根据已知批次信息校正后差异基因数目变多了，上调了99个，下调多了61个。根据预测的混杂因素校正后，上调多了96个，下调多了55个。从数目上看，根据已知批次和预测的混杂因素获得的差异基因是基本已知的。
画个Venn图，看下哪些基因是新增的差异基因，哪些基因批次校正后没差异了。
一个方式是采用代码，直接出图
从untrt上调基因Venn图可以看出，校正已知批次信息后既有新增的untrt上调差异基因，又丢失了之前的一部分untrt上调差异基因；校正预测的混杂因素后，大部分新增差异基因都与已知批次信息校正后的结果一致，但新增untrt上调差异基因少。
从untrt下调基因Venn图可以看出，校正预测的混杂因素后，新增39个差异基因；批次校正前鉴定为存在差异的40个基因在校正后被认为是非差异显著基因。
下面还是从这些基因的表达模式上看是否可以找到一些线索？
下图比对绘出了7种不同类型untrt上调的差异基因中随机选取1个绘制的表达模式比较图。
All_common代表了差异倍数特别大的基因，不论是否校正都可以检测出差异；不同类型批次信息校正后被检测视为差异的基因都有表达的本底差异；Uncorrect_specific的基因本底表达无固定模式。
下图比对绘出了7种不同类型untrt下调的差异基因表达分布，基本结论与上图类似。
另外一个导致SVA预测的批次与已知的批次效应校正后结果不同的原因也可能是我们只让SVA预测了2个混杂因素。留下2个去探索的问题，欢迎留言或投稿讨论：
处理批次因素最好的方式还是如前面所述将其整合到差异基因鉴定模型中，降低批次因素带来的模型残差的自由度。但一些下游分析，比如数据可视化，也需要直接移除效应影响的数据来展示，这时可以使用ComBat或removeBatchEffect函数来处理。
包含已知批次信息和预测的批次信息的样本属性文件 sampleFile2
加载需要的包
读入标准化后的表达矩阵和样品信息表
ComBat校正时考虑生物分组信息
校正后的PCA结果显示在PC1轴代表的差异变大了，PC2轴代表的差异变小了，不同来源的样本在PC2轴的分布没有规律了 (或者说成镜像分布了)。
ComBat校正时不考虑分组信息，也可以获得一个合理的结果，但是一部分组间差异被抹去了。
关于运行ComBat时是否应该添加关注的生物分组信息，即mod变量，存在一些争议。反对添加mod的人的担心是这么处理后，是否会强化生物分组之间的差异。支持添加mod的人是担心如果不添加mod那么去除批次时可能也会去除样本组之间的差异，尤其是实验设计不合理时。
如果是非平衡实验，类似我们在实验设计部分时提到的方案2，则没有办法添加mod变量，程序会报出Design matrix is not full rank类似的错误，这时是不能区分差异是来源于批次还是来源于样本，强行移除批次时，也会移除一部分或者全部样本分组带来的差异。这个在第一篇帖子处有两位朋友的留言讨论可以参考。
ComBat只能处理批次信息为l离散型分组变量的数据，不能处理sva预测出的连续性混杂因素。
如果批次信息有多个或者不是分组变量而是类似SVA预测出的数值混杂因素，则需使用limma的removeBatchEffect(这里使用的是SVA预测出的全部3个混杂因素进行的校正。)。
样品在PC1和PC2组成的空间的分布与ComBat结果类似，只是PC1能解释的差异略小一些。
removeBatchEffect运行时也可以不提供生物分组信息。
removeBatchEffect也可以跟ComBat一样，对给定的已知一个或多个定性批次信息进行校正。
不指定目标分组变量，结果也不受影响。
同时考虑批次、混杂因素和生物分组信息进行校正，校正后差异就全部集中在生物分组信息水平 (PC1)上了 (PC1 variance=100)。