生物信息学习的正确姿势

NGS系列文章包括NGS基础、转录组分析 （Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析 （ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程 （原理、代码和评述）)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集）等内容。

高通量数据中批次效应的鉴定和处理（一）

高通量数据中批次效应的鉴定和处理（二）

高通量数据中批次效应的鉴定和处理（三）- 如何设计尽量避免批次影响

高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应

高通量数据中批次效应的鉴定和处理（五）- 预测并校正可能存在的混杂因素

直接校正表达矩阵

处理批次因素最好的方式还是如前面所述将其整合到差异基因鉴定模型中，降低批次因素带来的模型残差的自由度。但一些下游分析，比如数据可视化，也需要直接移除效应影响的数据来展示，这时可以使用ComBat或removeBatchEffect函数来处理。

输入数据，标准化且log转换后的数据 ehbio.simpler.DESeq2.normalized.rlog.xls

id    untrt_N61311    untrt_N052611    untrt_N080611    untrt_N061011    trt_N61311    trt_N052611    trt_N080611    trt_N061011ENSG00000115414    18.02     18.62     17.83     18.45     17.95     18.54     18.15     18.50ENSG00000011465    17.79     18.36     17.96     18.52     17.79     18.23     17.84     18.47ENSG00000091986    17.15     17.74     16.41     17.59     17.27     17.79     16.88     17.76ENSG00000103888    15.56     16.90     15.88     16.42     15.94     17.43     17.38     17.05

包含已知批次信息和预测的批次信息的样本属性文件 sampleFile2

Samp    conditions    individual    sizeFactor    SV1    SV2    SV3untrt_N61311    untrt    N61311    1.0211325    -0.10060313    -0.4943517    -0.31643389untrt_N052611    untrt    N052611    1.1803986    0.01827734    -0.1701068    0.58841464untrt_N080611    untrt    N080611    1.1796083    -0.42949195    0.3756338    -0.08929556untrt_N061011    untrt    N061011    0.9232642    0.53452392    0.2413738    -0.17649091trt_N61311    trt    N61311    0.8939275    -0.12535603    -0.4956603    -0.36550102trt_N052611    trt    N052611    0.6709229    0.03588273    -0.151201    0.5914179trt_N080611    trt    N080611    1.3967665    -0.46668403    0.4413431    -0.07016903trt_N061011    trt    N061011    0.9462307    0.53345114    0.2529692    -0.16194213

加载需要的包

# 若缺少YSX包，则安装# BiocManager::install("Tong-Chen/YSX", update=F)suppressMessages(library(DESeq2))suppressMessages(library("RColorBrewer"))suppressMessages(library("gplots"))suppressMessages(library("amap"))suppressMessages(library("ggplot2"))suppressMessages(library("BiocParallel"))suppressMessages(library("YSX"))suppressMessages(library(sva))suppressMessages(library(ggfortify))suppressMessages(library(patchwork))suppressMessages(library(ggbeeswarm))suppressMessages(library(limma))

读入标准化后的表达矩阵和样品信息表

expr_File <- "ehbio.simpler.DESeq2.normalized.rlog.xls"expr_mat <- sp_readTable(expr_File, row.names=1)head(expr_mat)sampleFile <- "sampleFile2.txt"metadata <- sp_readTable(sampleFile, row.names=1)head(metadata)

使用ComBat校正

ComBat校正时考虑生物分组信息

biological_group = "conditions"batch = "individual"metadata[[biological_group]] <- factor(metadata[[biological_group]])metadata[[batch]] <- factor(metadata[[batch]])# 模型中引入关注的生物变量和其它非批次变量，保留生物差异和非批次差异modcombat = model.matrix(as.formula(paste('~', biological_group, sep=" ")), data=metadata)# ComBat需要的是matrixexpr_mat_batch_correct <- ComBat(dat=as.matrix(expr_mat), batch=metadata[[batch]], mod=modcombat)expr_mat_batch_correct <- as.data.frame(expr_mat_batch_correct)expr_mat_batch_correct[1:3,1:4]

untrt_N61311 untrt_N052611 untrt_N080611 untrt_N061011ENSG00000109906     5.107739      5.201200      5.293888      5.092441ENSG00000152583     7.002840      7.199577      7.281516      7.103559ENSG00000224114     6.428427      6.433242      6.236464      6.235308

校正后的PCA结果显示在PC1轴代表的差异变大了，PC2轴代表的差异变小了，不同来源的样本在PC2轴的分布没有规律了 (或者说成镜像分布了)。

sp_pca(expr_mat_batch_correct[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

ComBat校正时不考虑分组信息，也可以获得一个合理的结果，但是一部分组间差异被抹去了。

# ComBat需要的是matrixexpr_mat_batch_correct2 <- ComBat(dat=as.matrix(expr_mat), batch=metadata[[batch]], mod=NULL)expr_mat_batch_correct2 <- as.data.frame(expr_mat_batch_correct2)sp_pca(expr_mat_batch_correct2[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

关于运行ComBat时是否应该添加关注的生物分组信息，即mod变量，存在一些争议。反对添加mod的人的担心是这么处理后，是否会强化生物分组之间的差异。支持添加mod的人是担心如果不添加mod那么去除批次时可能也会去除样本组之间的差异，尤其是实验设计不合理时。

如果是非平衡实验，类似我们在实验设计部分时提到的方案2，则没有办法添加mod变量，程序会报出Design matrix is not full rank类似的错误，这时是不能区分差异是来源于批次还是来源于样本，强行移除批次时，也会移除一部分或者全部样本分组带来的差异。这个在第一篇帖子处有两位朋友的留言讨论可以参考。

ComBat只能处理批次信息为l离散型分组变量的数据，不能处理sva预测出的连续性混杂因素。

使用limma校正

如果批次信息有多个或者不是分组变量而是类似SVA预测出的数值混杂因素，则需使用limma的removeBatchEffect (这里使用的是SVA预测出的全部3个混杂因素进行的校正。)。

样品在PC1和PC2组成的空间的分布与ComBat结果类似，只是PC1能解释的差异略小一些。

SV = metadata[,c("SV1","SV2","SV3")]expr_mat_batch_correct_limma1 <- removeBatchEffect(expr_mat, covariates = SV, design=modcombat)sp_pca(expr_mat_batch_correct_limma1[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

removeBatchEffect运行时也可以不提供生物分组信息。

expr_mat_batch_correct_limma1 <- removeBatchEffect(expr_mat, covariates = SV)sp_pca(expr_mat_batch_correct_limma1[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

removeBatchEffect也可以跟ComBat一样，对给定的已知一个或多个定性批次信息进行校正。

expr_mat_batch_correct_limma2 <- removeBatchEffect(expr_mat, batch=metadata[[batch]], design=modcombat)sp_pca(expr_mat_batch_correct_limma2[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

不指定目标分组变量，结果也不受影响。

expr_mat_batch_correct_limma2 <- removeBatchEffect(expr_mat, batch=metadata[[batch]])sp_pca(expr_mat_batch_correct_limma2[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

同时考虑批次、混杂因素和生物分组信息进行校正，校正后差异就全部集中在生物分组信息水平 (PC1)上了 (PC1 variance=100)，应该是过拟合了，每组样本的基因表达都一致了。

expr_mat_batch_correct_limma3 <- removeBatchEffect(expr_mat, batch=metadata[[batch]], covariates = SV, design=modcombat)sp_pca(expr_mat_batch_correct_limma3[1:5000,], metadata,color_variable="conditions", shape_variable = "individual", coord_fixed_ratio=0) +aes(size=1) + guides(size = F)