使用sva包处理批次效应

生信修炼手册

发布于 2022-06-09 17:50:26

2.5K0

发布于 2022-06-09 17:50:26

文章被收录于专栏：生信修炼手册

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·

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-04-22，如有侵权请联系 cloudcommunity@tencent.com 删除

批量计算