专栏首页生信修炼手册使用sva包处理批次效应

使用sva包处理批次效应

SVA适用于高维数据的批次效应校正,支持以下数据

1. 基因芯片

2. RNA-seq

3. 甲基化表达谱

4. 其他表达量数据

提供了两种方法来处理不同的批次效应

1. 直接校正已知的batch effect, 使用ComBat 函数

2. 识别未知的batch effect,并校正,使用sva函数

需要注意的是,在校正批次效应之前,表达量数据必须经过归一化操作,而且去除了缺失的基因,比如在80%的样本中没有表达量的基因。

除了表达量矩阵外,还要求提供以下两种设计矩阵model matrix

1. null model,只包含adjustment variable, 即需要校正批次效应的变量

2. full model, 包含了adjustment variable + interest variable,包含所有的变量

两个函数的具体用法如下

1. ComBat

ComBat函数提供了基于经验贝叶斯框架的校正模型,具体的原理可以查看对应的文章

https://academic.oup.com/biostatistics/article/8/1/118/252073?login=false

对于该函数的使用,完整代码如下

> library(sva)
> library(bladderbatch)
> data(bladderdata)
> pheno = pData(bladderEset)
> edata = exprs(bladderEset)
> batch = pheno$batch
# 样本的metadata
> head(pheno)
             sample outcome batch cancer
GSM71019.CEL      1  Normal     3 Normal
GSM71020.CEL      2  Normal     2 Normal
GSM71021.CEL      3  Normal     2 Normal
GSM71022.CEL      4  Normal     3 Normal
GSM71023.CEL      5  Normal     3 Normal
GSM71024.CEL      6  Normal     3 Normal
# 样本表达量
> edata[1:5, 1:5]
          GSM71019.CEL GSM71020.CEL GSM71021.CEL GSM71022.CEL GSM71023.CEL
1007_s_at    10.115170     8.628044     8.779235     9.248569    10.256841
1053_at       5.345168     5.063598     5.113116     5.179410     5.181383
117_at        6.348024     6.663625     6.465892     6.116422     5.980457
121_at        8.901739     9.439977     9.540738     9.254368     8.798086
1255_g_at     3.967672     4.466027     4.144885     4.189338     4.078509
# combat校正
> combat_edata = ComBat(dat=edata, batch=batch, mod=NULL, par.prior=TRUE)
Found5batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Adjusting the Data
# 校正后的表达量
> combat_edata[1:5, 1:5]
          GSM71019.CEL GSM71020.CEL GSM71021.CEL GSM71022.CEL GSM71023.CEL
1007_s_at    10.086492     8.593022     8.735540     8.904952    10.279650
1053_at       5.519207     5.113332     5.158352     5.262212     5.265272
117_at        6.561742     6.432693     6.270198     6.220636     6.020385
121_at        8.789976     9.223917     9.316884     9.194757     8.670993
1255_g_at     3.899327     4.403982     4.092681     4.221116     4.060228
# 差异分析
# full model
> mod = model.matrix(~as.factor(cancer), data=pheno)
# null model
> mod0 = model.matrix(~1,data=pheno)
> pValuesComBat = f.pvalue(combat_edata,mod,mod0)
> qValuesComBat = p.adjust(pValuesComBat,method="BH")
> head(pValuesComBat)
  1007_s_at     1053_at      117_at      121_at   1255_g_at     1294_at
0.003051637 0.002971959 0.683888502 0.002948907 0.000299021 0.063957980
> head(qValuesComBat)
  1007_s_at     1053_at      117_at      121_at   1255_g_at     1294_at
0.009236571 0.009053201 0.755944614 0.009001439 0.001684293 0.109317762

2. SVA

SVA用于处理未知的batch effect, 完整的代码如下

> library(sva)
> library(bladderbatch)
> data(bladderdata)
> pheno = pData(bladderEset)
> edata = exprs(bladderEset)
> mod = model.matrix(~as.factor(cancer), data=pheno)
> mod0 = model.matrix(~1,data=pheno)
> n.sv = num.sv(edata,mod,method="leek")
> n.sv
[1] 2
> svobj = sva(edata,mod,mod0,n.sv=n.sv)
Number of significant surrogate variables is:  2
Iteration (out of 5 ):1  2  3  4  5 
> modSv = cbind(mod,svobj$sv)
> mod0Sv = cbind(mod0,svobj$sv)
> pValuesSv = f.pvalue(edata,modSv,mod0Sv)
> qValuesSv = p.adjust(pValuesSv,method="BH")
> head(pValuesSv)
   1007_s_at      1053_at       117_at       121_at    1255_g_at
2.003121e-03 5.984151e-05 5.127255e-01 9.139991e-05 1.112579e-08
     1294_at
1.218335e-02
>
> head(qValuesSv)
   1007_s_at      1053_at       117_at       121_at    1255_g_at
4.015794e-03 1.678983e-04 5.652893e-01 2.466591e-04 8.168562e-08
     1294_at
2.035400e-02

通过这两种方法,可以快速的校正数据中存在的批次效应,便于下游的表达量分析和差异分析。

·end·

文章分享自微信公众号:
生信修炼手册

本文参与 腾讯云自媒体分享计划 ,欢迎热爱写作的你一起参与!

作者:lzyg
原始发表时间:2022-04-22
如有侵权,请联系 cloudcommunity@tencent.com 删除。
登录 后参与评论
0 条评论

相关文章

  • 高通量数据中批次效应的鉴定和处理(二)

    通过样品的层级聚类热图+样品属性信息的注释来展示样品聚类结果有无受批次效应的影响。如下面右图中可见WT_1样品在聚类分支上与其它样品处于不同的分支,而从列注释图...

    生信宝典
  • 使用scran包的MNN算法来去除多个单细胞转录组数据批次效应

    这里我们使用一下scran包的 mutual nearest neighbors (MNNs)方法吧,主要就是读文档而已:https://bioconducto...

    生信技能树jimmy
  • 一文搞定高通量数据整合分析中批次效应的鉴定和处理

    批次效应表示样品在不同的批次处理和测量时引入的与生物状态不相关的系统性的技术偏差。很多因素都可能导致批次效应的产生,如不同实验条件、不同操作者、不同公司的试剂、...

    生信宝典
  • 多种批次效应去除的方法比较

    前面我在生信技能树推文:你确定你的差异基因找对了吗? 提出了文章的转录组数据的60个样品并没有按照毒品上瘾与否这个表型来区分,而是不同人之间的异质性非常高,这个...

    生信技能树
  • GEO数据库挖掘之多个芯片数据集的合并

    生信技能树
  • 去除批次效应好,还是RobustRankAggreg优?

    小洁老师在去除批次效应的探索文件里给出了两种方法,一个是用R包limma中的函数removeBatchEffect(),另一个是R包SVA中的函数ComBat(...

    生信技能树
  • 手把手教你多套GEO数据集合并

    各位科研芝士的小伙伴,TCGA、GEO数据库的挖掘现如今已经十分火爆。不可否认,现如今各种培训层出不穷,几乎都是给你一个代码让你去跑,却并没有让你真正懂其精髓。...

    百味科研芝士
  • 高通量数据中批次效应的鉴定和处理(六)- 直接校正表达矩阵

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

    生信宝典
  • 高通量数据中批次效应的鉴定和处理(五)- 预测并校正可能存在的混杂因素

    下面的方式也可以 (svaseq 是在 sva 的基础上对数据做了一个 log 转换;如果处理的是芯片数据,通常已经做过 log 换,直接使用 sva 即可)。

    生信宝典
  • 批次效应去除工具

    我们在进行公共数据挖掘的时候,经常会碰到要对多个数据集联合分析的时候,如果想要把这些数据放到一起进行分析的话,那么首先还是需要先去除批次效应才能进行分析的。之前...

    医学数据库百科
  • 药物预测之相关性为什么不行

    现在可以尝试一下理解药物预测的原理啦,首先我们看看仅仅是相关性,能不能搞定药物预测!

    生信技能树
  • R语言批数处理

    在很多实验的时候都会遇到不同批次的数据整合的情况,那么今天就给大家介绍一个测序数据的批次数据分析的R包sva。首先我们看下包的安装,以及内置数据的安装:

    一粒沙
  • 3分+非肿瘤纯生信免疫浸润分析还可以这么做!

    大家好,今天小编分享的是今年3月份发表在Diagnostics (Basel)(IF:3.1)的一篇非肿瘤纯生信文章,本文作者通过挖掘GEO数据库中骨关节炎相关...

    百味科研芝士
  • 使用DEseq2做转录组测序差异分析的时候顺便去除批次效应

    昨天的讨论:TCGA等大样本量差异分析该使用DEseq2还是edgeR呢? 让大家印象深刻,也有不少留言问到如果转录组测序数据集有批次效应该怎么办。所以我打个补...

    生信技能树
  • OSCA单细胞数据分析笔记12—Intergrating Datasets

    无论是scRNA-seq,还是Bulk RNA-seq,批次效应都是一个很头疼的问题,如何有效地校正、并且正确地使用校正后的数据是很值得讨论的分析点。

    生信技能树jimmy
  • python 使用 asyncio 包处理并发

    适合 asyncio 的协程要由调用方驱动,并由调用方通过 yield from 调用(语法过时了,新版的用 async / await ) 或者把协程传给 ...

    Michael阿明
  • 跟着全网第一个单细胞视频课程和配套习题学是最佳策略

    第 11-20题是使用3大R包从细胞的降维和分群到每个群细胞的功能注释,最后到公共数据库的注释,生存分析看不同细胞亚群的临床意义。

    生信技能树
  • ChAMP 包分析450K甲基化芯片数据(一站式)

    就有非常棒的一站式教程投稿,也因此我结识了优秀的六六,以及其教程大力推荐的R包作者,见:

    生信技能树

扫码关注腾讯云开发者

领取腾讯云代金券