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

生物信息学习的正确姿势

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)

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

首先加载所有的包

# 若缺少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：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    1

输入文件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    N061011

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

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

# Prefix for all output file
output_prefix = "ehbio.simplier"

# pipelineStar.sh或其它方式生成的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=1

数据读入和标准化

dds <- readscount2deseq(file, sampleFile, design=design, covariate = covariate)

normexpr <- deseq2normalizedExpr(dds, output_prefix=output_prefix)

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

样本聚类查看样品相似性，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="."))

主成分分析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)

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

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

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

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

# Prefix for all output file
output_prefix = "ehbio.simpler.batch"

# pipelineStar.sh生成的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=1

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

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

样本聚类查看样品相似性，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="."))

主成分分析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)

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

鉴定出差异基因，获得差异基因文件ehbio.simpler.batch.DESeq2.all.DE和其它可视化图表（暂时忽略）。

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

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

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

de_before_batch = sp_readTable("ehbio.simpler.DESeq2.all.DE", header=F)
de_before_batch$V2 = paste("Before_batch",de_before_batch$V2,sep="_")
table(de_before_batch$V2)

Beforebatch_untrt._higherThan.trt  Beforebatch_untrt._lowerThan.trt

                           398                                 466

de_after_batch = sp_readTable("ehbio.simpler.batch.DESeq2.all.DE", header=F)
de_after_batch$V2 = paste("After_batch",de_after_batch$V2,sep="_")
table(de_after_batch$V2)

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)

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_.trt

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

从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')

我们再选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')

DESeq2，edgeR和limma在考虑批次因素鉴定差异基因时基本操作是一致的，上面我们也完成和比较了已知批次的数据的差异基因鉴定。

后续还有2个问题：

1. DESeq2不校正表达矩阵自身的值，如果需要用到批次校正后的表达矩阵怎么做？

2. 如果不知道数据是否来源于同一个个体或是否有其他批次因素的影响，怎么处理？

测试数据和代码在 https://github.com/Tong-Chen/Bioinfo_course_R

送书

高级教程和生信基础知识在生信宝典往期推文中都有，赠送的书籍可以作为一类延伸阅读扩展知识面。

本期的留言主题是：您在生物信息学习过程中最喜欢哪种语言处理数据，为什么？欢迎转发朋友圈并留言评论。

留言得赞最高者将获得下面由北京大学出版社赞助的书籍（联系小编时请附上分享截图），欢迎分享，结果在下一期送书活动中公布：

Python数据分析圣经：通过核心概念剖析、45个“新手问答”、17个章节的“实训”、3个项目综合实战、50道Python面试题精选，教你轻松玩转数据分析与大数据处理。

往期精品(点击图片直达文字对应教程)

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集