高通量数据中批次效应的鉴定和处理(六)- 直接校正表达矩阵
直接校正表达矩阵
处理批次因素最好的方式还是如前面所述将其整合到差异基因鉴定模型中,降低批次因素带来的模型残差的自由度。但一些下游分析,比如数据可视化,也需要直接移除效应影响的数据来展示,这时可以使用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_N061011
ENSG00000115414 18.02 18.62 17.83 18.45 17.95 18.54 18.15 18.50
ENSG00000011465 17.79 18.36 17.96 18.52 17.79 18.23 17.84 18.47
ENSG00000091986 17.15 17.74 16.41 17.59 17.27 17.79 16.88 17.76
ENSG00000103888 15.56 16.90 15.88 16.42 15.94 17.43 17.38 17.05
包含已知批次信息和预测的批次信息的样本属性文件 sampleFile2
Samp conditions individual sizeFactor SV1 SV2 SV3
untrt_N61311 untrt N61311 1.0211325 -0.10060313 -0.4943517 -0.31643389
untrt_N052611 untrt N052611 1.1803986 0.01827734 -0.1701068 0.58841464
untrt_N080611 untrt N080611 1.1796083 -0.42949195 0.3756338 -0.08929556
untrt_N061011 untrt N061011 0.9232642 0.53452392 0.2413738 -0.17649091
trt_N61311 trt N61311 0.8939275 -0.12535603 -0.4956603 -0.36550102
trt_N052611 trt N052611 0.6709229 0.03588273 -0.151201 0.5914179
trt_N080611 trt N080611 1.3967665 -0.46668403 0.4413431 -0.07016903
trt_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需要的是matrix
expr_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_N061011
ENSG00000109906 5.107739 5.201200 5.293888 5.092441
ENSG00000152583 7.002840 7.199577 7.281516 7.103559
ENSG00000224114 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需要的是matrix
expr_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)
封面来源于:https://www.pexels.com/zh-cn/photo/264537/
- Go语言的小细节--map
- RBD至FileStore之所见(原理篇)
- FileStore压缩存储(优化篇)
- 深入学习Golang—channel
- 当Python字符串遇上MySQL
- Spring Boot 2.0正式发布,升还是不升呢?
- Spring Cloud构建微服务架构:分布式服务跟踪(入门)
- Spring Cloud构建微服务架构:分布式服务跟踪(跟踪原理)
- Spring Cloud Gateway真的有那么差吗?
- Spring Cloud构建微服务架构:消息驱动的微服务(核心概念)【Dalston版】
- Golang语言社区--【基础知识】循环
- Spring Cloud构建微服务架构:消息驱动的微服务(消费组)【Dalston版】
- Spring Cloud构建微服务架构:消息驱动的微服务(消费分区)【Dalston版】
- Spring Boot中使用LDAP来统一管理用户信息
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- php基础教程 第九步 自定义函数
- php基础教程 第十步 阶段性知识补充
- php基础教程 第十一步 面向对象
- php基础教程 第十一步 面向对象补充
- Serverless|Framework——图文玩转 AWS Lambda
- C++入门指南及实战 第一步 概述及经典HelloWorld
- C++入门指南及实战 第二步 HelloWorld及扩展详解
- 依托于GitLab持续集成基础配置和使用
- C++入门指南及实战 第三步 基本变量
- 最全总结 | 聊聊 Python 数据处理全家桶(配置篇)
- FlexSDK工具包的介绍与编译使用
- 《零基础看得懂的C语言入门教程 》——(二)C语言没那么难简单开发带你了解流程
- Ansible自动化运维学习笔记5
- Ansible自动化运维学习笔记4
- 《零基础看得懂的C语言入门教程 》——(三)轻轻松松理解第一个C语言程序