Question: about SVA: matrix of corrected expression values
0
gravatar for Bogdan
2.3 years ago by
Bogdan530
Palo Alto, CA, USA
Bogdan530 wrote:

Dear all,

I would need a bit of help please on SVA package for batch correction :

https://www.bioconductor.org/packages/release/bioc/html/sva.html

Assuming I have a matrix (let's call it DATA) of gene expression having 3 treatments (lets's say, CTRL, TREAT1, TREAT2), after correcting for batch effects in SVA, how could I obtain the matrix of CORRECTED VALUES ?

A piece of R code is :

DATA= as.matrix(gene.expression) ### here is the MATRIX DATA ;)

group = as.factor(rep(c("CTRL", "TREAT1", "TREAT2"), each=3))

mod1 = model.matrix(~group)
mod0 = cbind(mod1[,1])

svseq = svaseq(DATA, mod1, mod0, n.sv=1) 

and I would need to obtain the matrix of CORRECTED_DATA ? thanks !

combat sva • 1.1k views
ADD COMMENTlink modified 2.3 years ago • written 2.3 years ago by Bogdan530
Answer: about SVA: matrix of corrected expression values
2
gravatar for Jeff Leek
2.3 years ago by
Jeff Leek530
United States
Jeff Leek530 wrote:

Hi 

You would need to do this by performing a regression analysis on the svs estimated from svaseq and getting the residuals. If your downstream goal involves any type of significance test (differential analysis between groups for example) then we _do not_ recommend cleaning your data and then performing analysis on that cleaned data. This will cause problems with the calculation of the degrees of freeom for the statistical test. 

Instead we suggest incorporating the estimated surrogate variables as covariates in your analyses with limma, egdeR or DESeq after performing your svaseq analysis. 

 

Best

 

Jeff

ADD COMMENTlink written 2.3 years ago by Jeff Leek530
Answer: about SVA: matrix of corrected expression values
2
gravatar for bioinfo20014
2.3 years ago by
bioinfo2001420
bioinfo2001420 wrote:

I once found the function below to do what you want:

  # function from Back-estimating batch variables from SVA for ComBat?
  cleaningY = function(y, mod, svaobj) {
         X = cbind(mod, svaobj$sv)
           Hat = solve(t(X)%*%X)%*%t(X)
           beta = (Hat%*%t(y))
           P = ncol(mod)
           cleany = y-t(as.matrix(X[,-c(1:P)])%*%beta[-c(1:P),])
           return(cleany)
  }

You can also use the removeBatchEffect() function from the limma package and give the SVs as covariates:

svseq <- svaseq(counts_input, mod, mod0, n.sv=nsv) 

cov = cbind(svseq$sv[,1:nsv])

counts_deseq_sva <- removeBatchEffect(counts_input, covariates = cov)

I notice, however, that the results do differ between the two methods (the function and removeBatchEffect()). The function seems to work almost miraculously, but I think it makes sense, since you are directly regressing the SVs. removeBatchEffect() doesn't always give the best results in terms of "fixing" clustering issues.

 

See also this post DESeq2 - Acquiring batch-corrected values for PCA and hierarchical clustering about using rlog or vst transformed data as input to your batch effect removal.

ADD COMMENTlink written 2.3 years ago by bioinfo2001420
Answer: about SVA: matrix of corrected expression values
0
gravatar for Bogdan
2.3 years ago by
Bogdan530
Palo Alto, CA, USA
Bogdan530 wrote:

Thank you gentlemen, very helpful !

ADD COMMENTlink written 2.3 years ago by Bogdan530
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 152 users visited in the last hour