- CarieVinne
-
一、什么是批次效应
批次效应(batch effect),表示样品在不同批次中处理和测量产生的与试验期间记录的任何生物变异无关的技术差异。批次效应是高通量试验中常见的变异来源,受日期、环境、处理组、实验人员、试剂、平台等一些非生物因素的影响。
合并分析不同批次的数据时,平常的标准化方法不足以调整批次之间的差异。如果批次效应比较严重,这些差异就会干扰实验结果,我们就不能够判断得到的差异表达的基因是来源于想要研究的因素,还是和批次相关。
批次效应不能被消除,只有尽可能的降低。校正批次效应的目的是,减少批次之间的差异,尽量让多个批次的数据重新组合在一起,这样下游分析就可以只考虑生物学差异因素。
二、处理方法
目前已经有多种处理批次差异的方法。[2]
批次效应处理方法
三、哪种方法比较好
一项研究比较了6种去除批次效应的方法,其中包括ComBat方法(parametric prior method,ComBat_p和non-parametric method,ComBat_n)、代理变量法(Surrogate variable analysis,SVA)、基于比值的方法(Geometric ratio-based method,Ratio_G)、平均中心方法(Mean-centering,PAMR)和距离加权判别(Distance-weighted discrimination,DWD)方法,综合多种指标认为ComBat在精确性、准确性和整体性能方面(precision, accuracy and overall performance)总体优于其他方法。
四、ComBat方法的算法
模型的假设是基于位置和尺度(Location and scale,L/S)的调整。L/S调整可以定义为一系列广泛的调整,其中为数据在批次内的位置(均值)和/或规模(方差)假设了一个模型,然后调整批次以满足假设模型的规范。因此,L/S批次调整假设批次效应可以通过标准化批次之间的均值和方差来建模。这些调整可以从简单的基因范围的均值和方差标准化,到复杂的基因间线性或非线性调整。
模型:Yijg = αg + Xβg + γig + δigεijg
Yijg表示来自批次i的样品j的基因g的表达值。其中αg是基因g的平均表达值,X是样本条件的设计矩阵,βg是对应于X的回归系数向量。误差项εijg服从期望值为0和方差为σg的正态分布N(0,σg ),γig和δig表示批次i中基因g加法和乘法的批次效应。
算法分为3步: [4]
算法
五、推荐的解决方法
(1)根据分析目的确定批次效应处理方法:差异表达分析,在模型中添加批次因素;可视化,先对原数据进行校正,再使用校正后的数据进行分析。
(2)已知的批次,removeBatchEffect或ComBat;未知的批次,sva。
(3)removeBatchEffect和ComBat、sva输入的数据需要进行转换,例如取对数(rlog或logCPM)。
(4)read counts数据可使用ComBat-Seq或svaseq。
六、差异表达分析的批次校正
很多人以为去除批次效应是要改变你的表达矩阵,新的表达矩阵然后去走差异分析流程,其实大部分的差异分析流程包里面,人家内置好了考虑你的批次效应这样的混杂因素的函数用法设计。例如:
构建DESeq2对象时的设计公式: design = ~ batch + conditions
如果要合并不同批次的数据进行差异表达分析,建议直接把批次信息加入到构建模型里面。但是这种方法并没有改变原数据。如果你确实一定要亲眼看看批次效应到底是如何影响这个表达矩阵的,就需要对表达矩阵做另外的处理,例如removeBatchEffect或ComBat。但是处理之后会改变counts矩阵,之后就没办法走DESeq2差异分析流程啦,仅仅是为了拿到去除批次效应前后对比的表达矩阵而已。
PCA
七、使用limma的removeBatchEffect处理批次效应
removeBatchEffect这个函数用于进行聚类或无监督分析之前,移除与杂交时间或其他技术变异相关的批次效应。它是针对芯片设计的,因此不要直接使用read counts,数据需要经过一定的标准化操作,如log转化。
removeBatchEffect只用于衔接聚类、PCA等可视化展示,不要在线性建模之前使用。因为用矫正后的数据进行差异表达分析,有两个缺陷:一是批次因素和分组因素可能重叠,所以直接对原数据矫正批次可能会抵消一部分真实生物学因素;二是低估了误差。所以,如果想做差异表达分析,但数据中又有已知的批次问题,则最好把批次效应纳入线性模型中。
removeBatchEffect用法
removeBatchEffect
八、使用SVA的ComBat处理批次效应
sva这个R包可以处理已知的和未知的批次效应,sva函数可以通过构建高维数据集的代理变量,移除批次效应和其它所有不需要的变异。如果是芯片数据用sva,高通量测序数据使用svaseq。ComBat函数可以处理已知的批次效应。
ComBat用法
ComBat
九、使用sva处理未知的批次
SVA具有移除批次效应和高通量测序中其它不需要的变异的功能。使用sva识别和构建高维数据集代理变量(surrogate variables),代理变量是由高维数据直接构建的协变量(covariates),可用于后续分析,以调整未知、未建模或潜在的噪声源。
sva函数的输出本身就是代理变量。它们可以包含在全模型矩阵和空模型矩阵中,然后与数据矩阵一起传递给SVA包中的f.pvalue函数,以计算参数F检验p值,从而调整代理变量。
根据经验,当存在大量已知或未知的潜在混杂因素时,代理变量调整(sva)可能更合适。当一个或多个生物学分组已知是异质的,并且有已知的批次变量时,直接调整(ComBat)可能更合适。
sva考虑了两种类型的变量:调整变量和感兴趣的变量。例如,感兴趣变量(variables of interest)可能是癌症组与对照组;调整变量(adjustment variables)可能是病人的年龄、性别、测序时间。
建立两个模型矩阵:全模型(full model)和空模型(null model)。空模型包括所有调整变量,而不包含感兴趣的变量;全模型包括所有调整变量和感兴趣的变量。我们将试图分析感兴趣的变量与基因表达之间的关联,调整调整变量。模型矩阵可以使用model.matrix创建。sva的目标是消除所有不需要的变异来源,同时通过mod中包含的主要变量来检测对比。
注意:在我们最初的工作中,使用识别函数来测量在近似对称和连续尺度上的数据。对于通常表示为counts的测序数据,更合适的模型可能涉及使用适度的对数函数。例如,我们首先用log(counts+1)转换基因表达数据。
用法:
(1)使用sva得到代理变量
(2)使用f.pvalue函数调整代理变量(calculate parametric F-test P-values adjusted for surrogate variables)
(3)sva可以与差异表达分析程序一起使用,如limma、DESeq2。
References
[1] Chen C, Grennan K, Badner J, Zhang D, Gershon E, Jin L, et al. Removing batch effects in analysis of expression microarray data: an evaluation of six batch adjustment methods[J]. PloS One. 2011;6(2):e17238.
[2] 李飒,赵毅强.基因表达数据批次效应去除方法的研究进展[J].南京农业大学学报,2019,42(03):389-397.
[3] Chen C, Grennan K, Badner J, Zhang D, Gershon E, Jin L, et al. Removing batch effects in analysis of expression microarray data: an evaluation of six batch adjustment methods[J]. PloS One. 2011;6(2):e17238.
[4]陈天成,侯艳,李康.基因组学数据整合中的批次效应移除算法[J].中国卫生统计,2016,33(03):527-529+533.
[5] Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods[J]. Biostatistics. 2007;8(1):118-27.
[6] Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments[J]. Bioinformatics. 2012;28(6):882-3.