SVA package - for removing batch effect after DESeq - error in ComBat
3
0
Entering edit mode
GFM ▴ 20
@gfm-8326
Last seen 10 months ago
European Union

Hello,
I am trying to use the SVA package to obtain values after batch effect removal in RNA seq data (using DESeq rld counts),

My design:

design
                   strain batch
Late_mut_count   Late_mut     a
mid_mut_count     mid_mut     a
mid_WT_count       mid_WT     a
Young_mut_count Young_mut     a
Young_WT_count   Young_WT     a
mid2_mut_count    mid_mut     b
mid2_WT_count      mid_WT     b
mud_count             mud     b
y2_mut_count    Young_mut     b

y2_WT_count      Young_WT     b

 

mod = model.matrix(~strain, data=design)

> mod
                (Intercept) strainWT
Late_mut_count            1        0
mid_mut_count             1        0
mid_WT_count              1        1
Young_mut_count           1        0
Young_WT_count            1        1
mid2_mut_count            1        0
mid2_WT_count             1        1
y2_mut_count              1        0
y2_WT_count               1        1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$strain
[1] "contr.treatment"

 

I am using the DESeq's rld transformation: For the first columns, it looks like this:

Late_mut_count mid_mut_count mid_WT_count Young_mut_count Young_WT_count

gene1       2.719515      2.722055     2.742202        2.673428       2.678986
batch = design$batch



Now, I am trying to use ComBat:

combat_edata = ComBat(dat=rldCounts, batch=batch, mod=mod, numCovs=NULL, par.prior=TRUE)

But I get the following error:

Error in ComBat(dat = rldCounts, batch = batch, mod = mod, numCovs = NULL,  : 
  unused argument (numCovs = NULL)

 

Can you possibly help with this error please?

Thanks

sva combat • 2.6k views
ADD COMMENT
1
Entering edit mode

Also: how did you end up filling in the code in this question? I mean, where did you copy/paste from?

I went back in to try and clean up its formatting a bit, and it's all sorts of crazy in there ...

ADD REPLY
0
Entering edit mode
@steve-lianoglou-2771
Last seen 3 days ago
Denali

Don't forget to include the output of sessionInfo() so we know what version you're working with, but anyway: the error is telling you exactly what the problem is.

When I do ?ComBat using the latest version of sva, the help page will show you that the ComBat function only takes five arguments:

  • dat
  • batch
  • mod
  • par.prior
  • prior.plots

None of which is numCovs.

ADD COMMENT
0
Entering edit mode
GFM ▴ 20
@gfm-8326
Last seen 10 months ago
European Union

Thanks a lot!!!

Now I get a different error:

combat_edata = ComBat(dat=rldCounts, batch=batch, mod=mod, par.prior=TRUE)
Found 2 batches
Adjusting for 2 covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Error in while (change > conv) { : missing value where TRUE/FALSE needed

Can you help with this please?

===========================================

Here is the sessionInfo:

sessionInfo()
R version 3.2.1 (2015-06-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8 x64 (build 9200)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DESeq2_1.8.1              RcppArmadillo_0.5.200.1.0 Rcpp_0.11.6               GenomicRanges_1.20.5     
 [5] GenomeInfoDb_1.4.1        IRanges_2.2.5             S4Vectors_0.6.1           BiocGenerics_0.14.0      
 [9] BiocInstaller_1.18.3      sva_3.14.0                genefilter_1.50.0         mgcv_1.8-6               
[13] nlme_3.1-121             
ADD COMMENT
0
Entering edit mode
GFM ▴ 20
@gfm-8326
Last seen 10 months ago
European Union

OK solved this by filtering out genes with zero counts.
Thanks a lot for the help.

ADD COMMENT

Login before adding your answer.

Traffic: 378 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6