Hi Ying,

The problem is that model.matrix uses only the complete cases when building
the model matrix. So any of the samples that have an NA for even one
variable are removed from that model matrix. Since there are no NAs for the
batch variable, none of the variables are removed.  whereas the batch
variable has no NA values so none of the samples are removed.

One way to fix this problem is to choose a subset of variables for which
you have a set of samples with no NAs and run ComBat or SVA on only those
samples.

An alternative is to remove the batch variable using a linear model (see
the section in the SVA vignette for example on how to do this). Then
analyze the residuals after subtracting out the batch effect. You should
exercise caution when using this approach though, since you will be
analyzing a set of data that you have already done one regression on - this
reduces variability in the expression values and can lead to overestimated
significance.

Best,

Jeff



On Thu, Mar 15, 2012 at 3:24 AM, ying chen <ying_chen@live.com> wrote:

>
> Hi guys,
> I try to use ComBat function in SVA package to correct batch effect of my
> gene expression dataset, but I keep getting error messages. I think my
> problem may be that the model.matrix is not set correctly.
> I compiled several breast cancer datasets together and ended up with a big
> ExpressionSet.
> The pheno data is like this:
>                    Batch      Mol_type  ER     PR    Her2     TN
> Histologyarr01.cel    GSE01   Basal           ER-    PR-   Her2-   Yes
>  Breast_tumor_NOSarr02.cel    GSE02   LumA         ER+   NA   NA       No
>
> Breast_tumor_NOS.............................................................................................................................................................................................................................................
> There are many NAs. The only columns that do not have a single NA are
> Batch and Histology columns.
> I tried to set model.matrix in three ways and all got errors.
>
>
>
>
>
>
> > pheno <-
> pData(MyBreastExpressionSet)
>
> > edata <-
> exprs(MyBreastExpressionSet)
>
> # 1.  use two variables of interest: Mol_type and Histology
> > mod01 <-
> model.matrix(~as.factor(pheno$Mol_type)+as.factor(pheno$Histology))
> > batch <- pheno$Batch
>
> > combat_edata <-
> ComBat(dat=edata,batch=batch,mod=mod01,numCovs=NULL,par.prior=TRUE)
>
> Found 21 batches
>
> Found 48  categorical covariate(s)
>
> Error in tmp[, i] <- vec ==
> levels(vec)[i + start - 1] :
>
>  subscript out of bounds
>
> In addition: Warning message:
>
> In cbind(mod, batch) :
>
>  number of rows of result is not a
> multiple of vector length (arg 2)
>
> # 2.  use one variable of interest: Mol_type
>
>
>
>
>
> > mod02 <-
> model.matrix(~as.factor(pheno$Mol_type))
>
> > combat_edata <-
> ComBat(dat=edata,batch=batch,mod=mod02,numCovs=NULL,par.prior=TRUE)
>
> Found 21 batches
>
> Found 8  categorical covariate(s)
>
> Standardizing Data across genes
>
> Error in solve(t(design) %*% design)
> %*% t(design) %*% t(as.matrix(dat)) :
>
>  non-conformable arguments
>
> In addition: Warning message:
>
> In cbind(mod, batch) :
>
>  number of rows of result is not a
> multiple of vector length (arg 2)
>
> # 3. use one variable of interest: Histology
>
>
>
>
>
> > mod03 <-
> model.matrix(~as.factor(pheno$Histology))
> > combat_edata <-
> ComBat(dat=edata,batch=batch,mod=mod03,numCovs=NULL,par.prior=TRUE)
>
> Found 47 batches
>
> Found 40  categorical covariate(s)
>
> Standardizing Data across genes
>
> Error in solve.default(t(design) %*%
> design) :
>
>  Lapack routine dgesv: system is
> exactly singular
> Only during the last try the program found the correct number of batches
> (47). I think it may because of the NAs in Mol_type column. But I am not
> sure.
> What did I do wrong?
> Any suggestion is highly appreciated.
> Thanks,
>
> Ying
>
>
>
>
>
>
> >
> sessionInfo()
>
> R version 2.14.2 (2012-02-29)
>
> Platform: x86_64-pc-linux-gnu (64-bit)
>
>
>
> locale:
>
>  [1] LC_CTYPE=en_US.UTF-8
> LC_NUMERIC=C
>
>  [3] LC_TIME=en_US.UTF-8
> LC_COLLATE=en_US.UTF-8
>
>  [5] LC_MONETARY=en_US.UTF-8
> LC_MESSAGES=en_US.UTF-8
>
>  [7] LC_PAPER=C
> LC_NAME=C
>
>  [9] LC_ADDRESS=C
> LC_TELEPHONE=C
>
> [11] LC_MEASUREMENT=en_US.UTF-8
> LC_IDENTIFICATION=C
>
>
>
> attached base packages:
>
> [1] splines   stats     graphics
> grDevices utils     datasets  methods
>
> [8] base
>
>
>
> other attached packages:
>
> [1] pamr_1.54        survival_2.36-12
> cluster_1.14.2   sva_3.0.3
>
> [5] mgcv_1.7-13      corpcor_1.6.2
> Biobase_2.14.0
>
>
>
> loaded via a namespace (and not
> attached):
>
> [1] grid_2.14.2    lattice_0.20-6
> Matrix_1.0-4   nlme_3.1-103
>
> >
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

