高通量数据中批次效应的鉴定和处理(六)- 直接校正表达矩阵
直接校正表达矩阵
处理批次因素最好的方式还是如前面所述将其整合到差异基因鉴定模型中,降低批次因素带来的模型残差的自由度。但一些下游分析,比如数据可视化,也需要直接移除效应影响的数据来展示,这时可以使用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/
- 洛谷P3391 【模板】文艺平衡树(Splay)(FHQ Treap)
- 12.python进程协程异步IO
- POJ3622 Gourmet Grazers(FHQ Treap)
- 洛谷P3201 [HNOI2009]梦幻布丁
- 洛谷P3374 【模板】树状数组 1(CDQ分治)
- 自然语言处理基础知识1. 分词(Word Cut)2. 词性标注(POS Tag)3.自动标注4.文本分类5.评估6.从文本提取信息7.分析句子结构《python自然语言处理》各章总结:
- 洛谷P3384 【模板】树链剖分
- 洛谷P2147 [SDOI2008]Cave 洞穴勘测
- linux基础
- 洛谷P3178 [HAOI2015]树上操作
- Numpy 修炼之道 (6)—— 复制和视图
- 事务日志已满,原因为“ACTIVE_TRANSACTION”
- 【 关关的刷题日记46】Leetcode 28. Implement strStr()
- Python的机器学习库之Sklearn快速入门1.基本概述2.入门实践3.部分结果
- 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 数组属性和方法
- Hive on spark下insert overwrite partition慢的优化
- 系统学习Lambda表达式
- 一文搞懂 Flink Kafka Consumer 类两阶段提交
- 在 Nowin 下运行 ASP.NET 5 Beta 2
- Bytom侧链Vapor源码浅析-节点出块过程
- Kubernetes Pod OOM 排查日记
- Netty之旅:你想要的NIO知识点,这里都有!
- (数据科学学习手札92)利用query()与eval()优化pandas代码
- Spring Boot 集成 Elasticsearch 实战
- Python之错误和异常、模块(基础系列第四篇)
- Spark存储Parquet数据到Hive,对map、array、struct字段类型的处理
- 为什么这条异常没有上报? HTTP 429
- 三问Spring事务:解决什么问题?如何解决?存在什么问题?
- 从 OAuth2 服务器获取授权授权
- NHibernate 代码映射实体类