送书

发布时间：2026-01-11 06:37

生物信息学习的正确姿势

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

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

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

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

如何在差异基因鉴定过程中移除批次效应

在我们之前的文章 DESeq2差异基因分析和批次效应移除中也提到了用如下方式构建设计矩阵，以便在差异基因分析过程中移除批次效应的影响。

ddsFullCountTable <- DESeqDataSetFromMatrix(countData = data, colData = sample, design= ~ batch + conditions) dds <- DESeq(ddsFullCountTable)1.2.3.4.

下面我们以一个具体例子实战（配对样品处理前后基因表达的变化）和检验下效果。为了演示批次效应的影响，大部分代码做了封装，我们只关心核心的地方。如果自己对封装的代码赶兴趣，可以自行查看函数源码。

首先加载所有的包

# 若缺少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))1.2.3.4.5.6.7.8.9.10.11.12.13.

输入文件1：reads count矩阵 (ehbio_trans.Count_matrix.xls)，格式如下：

ENSG untrt_N61311 untrt_N052611 untrt_N080611 untrt_N061011 trt_N61311 trt_N052611 trt_N080611 trt_N061011 ENSG00000223972 1 0 0 0 0 1 0 0 ENSG00000227232 13 25 23 24 12 12 22 22 ENSG00000278267 0 5 3 4 2 4 3 11.2.3.4.

输入文件2：实验设计信息表 (sampleFile): conditions为处理条件（untrt是对照, trt是加药处理），individual标记样品的个体来源 (4个个体：N61311、N052611、N080611、N061011)。

Samp conditions individual untrt_N61311 untrt N61311 untrt_N052611 untrt N052611 untrt_N080611 untrt N080611 untrt_N061011 untrt N061011 trt_N61311 trt N61311 trt_N052611 trt N052611 trt_N080611 trt N080611 trt_N061011 trt N0610111.2.3.4.5.6.7.8.9.

不考虑批次因素直接进行差异基因分析

初始化，定义输入、输出和参数

# Prefix for all output file output_prefix = "ehbio.simplier" # 或其它方式生成的reads count 文件，行为基因，列为样品 file = "ehbio_trans.Count_matrix.xls" # 分组信息表 sampleFile = "sampleFile" # 分组信息所在列名字 covariate = NULL # covariate = "batch" design="conditions" # 输入数据类型，salmon结果或reads count 矩阵 type="readscount" # 差异基因参数 padj=0.05 log2FC=11.2.3.4.5.6.7.8.9.10.11.12.13.14.15.16.

数据读入和标准化

dds <- readscount2deseq(file, sampleFile, design=design, covariate = covariate) normexpr <- deseq2normalizedExpr(dds, output_prefix=output_prefix)1.2.3.

[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 factors

estimating dispersions

gene-wise dispersion estimates

mean-dispersion relationship

final dispersion estimates

fitting model and testing

[1] “Output normalized counts”

[1] “Output rlog transformed normalized counts”

检查数据标准化效果: 标准化后基因在不同样品的表达分布越均一越好。从下图看不出存在批次效应的影响。

# normalizedExpr2DistribBoxplot(normexpr, # saveplot=paste(output_prefix, "DESeq2.normalizedExprDistrib.pdf", sep=".")) normalizedExpr2DistribBoxplot(normexpr)1.2.3.

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._人工智能

样本聚类查看样品相似性，trt组和untrt组区分明显 (聚类采用的不同基因数目、聚类参数都可能影响聚类结果)

# clusterSampleHeatmap2(normexpr$rlog, # cor_file=paste(output_prefix, "DESeq2.sampleCorrelation.txt", sep="."), # saveplot=paste(output_prefix, "DESeq2.sampleCorrelation.pdf", sep=".")) # 根据前5000个表达变化幅度最大的基因进行聚类分析 clusterSampleHeatmap2(normexpr$rlog[1:5000,], cor_file=paste(output_prefix, "DESeq2.sampleCorrelation.txt", sep=".")) clusterSampleUpperTriPlot(normexpr$rlog[1:5000,], cor_file=paste(output_prefix, "DESeq2.sampleCorrelation.txt", sep="."))1.2.3.4.5.6.

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._数据可视化_02

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._人工智能_03

主成分分析PCA查看样品相似性，发现在PC1轴上，样品按处理条件区分开；在PC2轴上，样品按个体区分开，不同的个体是影响样品基因表达差异的一个重要因素。

metadata = as.data.frame(colData(dds)) sp_pca(normexpr$rlog[1:5000,], metadata, color_variable="conditions", shape_variable = "individual") + aes(size=1) + guides(size = F)1.2.

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._数据可视化_04

先鉴定出差异基因，获得差异基因文件和其它可视化图表（暂时忽略）。

multipleGroupDEgenes(dds, design=design, output_prefix=output_prefix, padj=padj, log2FC=log2FC)1.

考虑已知的批次因素进行差异基因分析

初始化，定义输入、输出和参数（注意covariate变量使用individual列作为了批次信息）

# Prefix for all output file output_prefix = "ehbio.simpler.batch" # 生成的reads count 文件，行为基因，列为样品 file = "ehbio_trans.Count_matrix.xls" # 分组信息表 sampleFile = "sampleFile" # 分组信息所在列名字 # covariate = NULL # ********* covariate = "individual" design="conditions" # 输入数据类型，salmon结果或reads count 矩阵 type="readscount" # 差异基因参数 padj=0.05 log2FC=11.2.3.4.5.6.7.8.9.10.11.12.13.14.15.16.17.

数据读入和标准化，并检查数据标准化效果: 标准化后基因在不同样品的表达分布越均一越好（此图略过，与上面的表达分布图一致）。

dds <- readscount2deseq(file, sampleFile, design=design, covariate = covariate) normexpr <- deseq2normalizedExpr(dds, output_prefix=output_prefix) normalizedExpr2DistribBoxplot(normexpr)1.2.3.

样本聚类查看样品相似性，trt组和untrt组区分明显 (此部分结果略过，与上面的聚类结果一致)

主成分分析PCA查看样品相似性，发现在PC1轴上，样品按处理条件区分开；在PC2轴上，样品按个体区分开，表明不同的个体可能会对后续的差异基因分析造成影响。这个结果也与我们前面不考虑批次因素的结果是一样的。

metadata = as.data.frame(colData(dds)) sp_pca(normexpr$rlog[1:5000,], metadata, color_variable="conditions", shape_variable = "individual") + aes(size=1) + guides(size = F)1.2.

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._人工智能_05

是不是批次变量加错了呢，还是添加的批次变量未生效？可以说都不是，操作没问题，只是DESeq2处理时只在差异分析模型中考虑批次效应信息，而不会直接校正表达矩阵。那我们先看下加了批次后差异分析的结果怎样，后续我们再讲如何校正表达矩阵。

鉴定出差异基因，获得差异基因文件和其它可视化图表（暂时忽略）。

multipleGroupDEgenes(dds, design=design, output_prefix=output_prefix, padj=padj, log2FC=log2FC)1.

比较批次校正前后差异基因变化

校正后，差异基因数目变多了，上调多了99个，下调多了61个。不过数目变化，也说明不了太多问题。

de_before_batch = sp_readTable("", header=F) de_before_batch$V2 = paste("Before_batch",de_before_batch$V2,sep="_") table(de_before_batch$V2)1.2.3.

Beforebatch_untrt._higherThan.trt Beforebatch_untrt._lowerThan.trt

398 466

de_after_batch = sp_readTable("", header=F) de_after_batch$V2 = paste("After_batch",de_after_batch$V2,sep="_") table(de_after_batch$V2)1.2.3.

Afterbatch_untrt._higherThan.trt Afterbatch_untrt._lowerThan.trt

497 527

画个Venn图，看下哪些基因是新增的差异基因，哪些基因批次校正后没差异了。这里就不写代码了，采用在线工具http://www.ehbio.com/test/venn/#/ 来做，准备在线工具所需的文件，一个两列格式的文件：第一列为基因名，第二列为基因的上下调状态。

all_de = rbind(de_before_batch, de_after_batch) # 随机查看6行，信息代表更全面 all_de[sample(1:nrow(all_de),6),] # 结果存储到文件中 sp_writeTable(all_de, file="Compare_de_gene_beofore_and_after_batch.txt", keep_rownames = F, col.names = F)1.2.3.4.5.

ENSG00000114270 After_batch_untrt._lowerThan_.trt ENSG00000102935 After_batch_untrt._higherThan_.trt ENSG00000131370 Before_batch_untrt._higherThan_.trt ENSG00000174944 After_batch_untrt._lowerThan_.trt ENSG00000157510 After_batch_untrt._lowerThan_.trt ENSG00000116711 After_batch_untrt._higherThan_.trt1.2.3.4.5.6.

拷贝文件数据到网站数据输入处：

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._数据可视化_06

从Venn图可以看出，批次校正后既有新增的差异基因，又丢失了之前的一部分差异基因，那么哪个方式更合理呢？

选择1个批次校正后检测为上调的基因和1个批次校正后检测为下调的基因，观察下其表达模式。从下图可以看出，这些基因具有明显的个体表达一致性。ENSG00000163394基因在每个个体来源的样本中处理后表达都上调了近4倍，但是其本底表达在不同个体中却差异较大。如其在N080611个体（蓝色线）中表达整体偏低，药物处理后表达虽然有上调但表达值却低于其在N061011个体（绿色线）处理前的表达。从这两个例子可以看出，考虑到每个个体的基准表达水平不同，最终获得的差异倍数会有较高的方差。批次校正后解决了样品个体来源基因本底表达差异的影响，获得的差异基因倍数方差会变小，所以检测出更多差异基因，理论上也是更可靠的方式。（这个在之前文章典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集也有阐述。）

ENSG00000163394 = data.frame(Expr=normexpr$rlog["ENSG00000163394",]) p1 <- sp_boxplot(ENSG00000163394, melted=T, metadata=metadata, xvariable = "conditions", yvariable = "Expr", jitter_bp = T, group_variable_for_line = "individual") ENSG00000221866 = data.frame(Expr=normexpr$rlog["ENSG00000221866",]) p2 <- sp_boxplot(ENSG00000221866, melted=T, metadata=metadata, xvariable = "conditions", yvariable = "Expr", jitter_bp = T, group_variable_for_line = "individual") p1 + p2 + plot_layout(guide = 'collect')1.2.3.4.5.6.7.

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._编程语言_07

我们再选2个批次校正前鉴定为有差异、批次校正后鉴定为无差异的基因观察下其表达模式。这两个基因的表达模式没看出存在个体本底的一致变化差异。处理前后在不同个体中变化幅度不一，可能是被动变化。但这些基因一定是没有差异吗？我个人也下不出结论，后续得结合其功能再做判断了。

ENSG00000109689 = data.frame(Expr=normexpr$rlog["ENSG00000109689",]) p1 <- sp_boxplot(ENSG00000109689, melted=T, metadata=metadata, xvariable = "conditions", yvariable = "Expr", jitter_bp = T, group_variable_for_line = "individual", title="ENSG00000109689") ENSG00000137124 = data.frame(Expr=normexpr$rlog["ENSG00000137124",]) p2 <- sp_boxplot(ENSG00000137124, melted=T, metadata=metadata, xvariable = "conditions", yvariable = "Expr", jitter_bp = T, group_variable_for_line = "individual", title="ENSG00000137124") p1 + p2 + plot_layout(guide = 'collect')1.2.3.4.5.6.7.

送书|高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应..._数据可视化_08