【原】使用sva包处理批次效应

生信修炼手册 2022-04-22

展开全文

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./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 cancerGSM71019.CEL      1  Normal     3 NormalGSM71020.CEL      2  Normal     2 NormalGSM71021.CEL      3  Normal     2 NormalGSM71022.CEL      4  Normal     3 NormalGSM71023.CEL      5  Normal     3 NormalGSM71024.CEL      6  Normal     3 Normal# 样本表达量> edata[1:5, 1:5]          GSM71019.CEL GSM71020.CEL GSM71021.CEL GSM71022.CEL GSM71023.CEL1007_s_at    10.115170     8.628044     8.779235     9.248569    10.2568411053_at       5.345168     5.063598     5.113116     5.179410     5.181383117_at        6.348024     6.663625     6.465892     6.116422     5.980457121_at        8.901739     9.439977     9.540738     9.254368     8.7980861255_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)Found5batchesAdjusting for0covariate(s) or covariate level(s)Standardizing Data across genesFitting L/S model and finding priorsFinding parametric adjustmentsAdjusting the Data# 校正后的表达量> combat_edata[1:5, 1:5]          GSM71019.CEL GSM71020.CEL GSM71021.CEL GSM71022.CEL GSM71023.CEL1007_s_at    10.086492     8.593022     8.735540     8.904952    10.2796501053_at       5.519207     5.113332     5.158352     5.262212     5.265272117_at        6.561742     6.432693     6.270198     6.220636     6.020385121_at        8.789976     9.223917     9.316884     9.194757     8.6709931255_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_at0.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_at0.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:  2Iteration (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_at2.003121e-03 5.984151e-05 5.127255e-01 9.139991e-05 1.112579e-08     1294_at1.218335e-02>> head(qValuesSv)   1007_s_at      1053_at       117_at       121_at    1255_g_at4.015794e-03 1.678983e-04 5.652893e-01 2.466591e-04 8.168562e-08     1294_at2.035400e-02