[BioC] How to perform differential expression analysis on a pre-normalized dataset?

[James W. MacDonald]

Hi Karthi,

On 9/29/2010 3:19 AM, Karthi Sivaraman wrote:
>    Dear List,
>
> I am sure this has been discussed before in the list, so I would be
> happy even if someone points me to the right previous discussion.
>
> I have a Nimblegen gene expression dataset as .pairs file. Using dummy
> 2nd channel reading, I have managed to normalize this dataset (vsn
> normalization). lmFit also works with the normalized dataset. However,
> when I try to perform eBayes fit, I am hitting a snag. eBayes fit
> returns a lot of error at this last stage.
>
> I have given the steps / results from each step in this workflow. If
> someone has already faced this issue and resolved it - i would be happy
> to try it. Also, if there is a basic mistake I am making, can someone
> please let me know. Thanks. -- k.
> ===========================================================*
>   >  library(limma)
>   >  target=readTargets('list')
>   >
> x=read.maimages(target$FileName,columns=list(R="PM",G="PM"),annotation=c("PROBE_ID","X","Y"))
> Read GSM519643.pair.val
> Read GSM519644.pair.val
> Read GSM519645.pair.val
> Read GSM519646.pair.val
> Read GSM519647.pair.val
> Read GSM519648.pair.val
> Read GSM519649.pair.val
> Read GSM519650.pair.val
> Read GSM519651.pair.val
>   >  x$R=NULL
>   >  xNorm=backgroundCorrect(x,method="normexp")
>   >  library(arrayQualityMetrics)
> Loading required package: affyPLM
> Loading required package: affy
> Loading required package: Biobase
>
> Welcome to Bioconductor
>
>     Vignettes contain introductory material. To view, type
>     'openVignette()'. To cite Bioconductor, see
>     'citation("Biobase")' and for packages 'citation(pkgname)'.
>
> Loading required package: gcrma
> Loading required package: preprocessCore
>
> Attaching package: 'affyPLM'
>
>
>       The following object(s) are masked from package:stats :
>
>        resid,
>        residuals,
>        weights
>
>   >  xEsetNorm=new("ExpressionSet",exprs=xNorm$G)
>   >  library(vsn)
>   >  xEsetVsn=vsn(xEsetNorm)
> Note:
> The function 'vsn' has been superseded by 'vsn2'.
> The function 'vsn' remains in the package for backward compatibility,
> but for new projects, please use 'vsn2'.

You might read (and follow) the note...

>
> vsn: 389307 x 9 matrix (1 stratum). 100% done.
>   >  design=model.matrix(~ -1+factor(c(1,1,1,2,2,2,3,3,3,))
> + )
> Error in factor(c(1, 1, 1, 2, 2, 2, 3, 3, 3, )) :
>     element 9 is empty;
>      the part of the args list of 'c' being evaluated was:
>      (1, 1, 2, 2, 2, 3, 3, 3, )
>   >  design=model.matrix(~ -1+factor(c(1,1,1,2,2,2,3,3,3)))
>   >  colnames(design)=c('Ctrl','Chlr1','Chlr2')
>   >  design
>     Ctrl Chlr1 Chlr2
> 1    1     0     0
> 2    1     0     0
> 3    1     0     0
> 4    0     1     0
> 5    0     1     0
> 6    0     1     0
> 7    0     0     1
> 8    0     0     1
> 9    0     0     1
> attr(,"assign")
> [1] 1 1 1
> attr(,"contrasts")
> attr(,"contrasts")$`factor(c(1, 1, 1, 2, 2, 2, 3, 3, 3))`
> [1] "contr.treatment"
>
>   >  xFit=lmFit(xEsetVsn,design)
>   >  contmat=makeContrasts(Chlr1-Ctrl,Chlr2-Ctrl,levels=design)
>   >  contmat
>          Contrasts
> Levels  Chlr1 - Ctrl Chlr2 - Ctrl
>     Ctrl            -1           -1
>     Chlr1            1            0
>     Chlr2            0            1
>   >  xFit2=eBayes(xFit
> xFit
>   >  xFit2=eBayes(xFit,contmat)

There are no arguments for eBayes() that take a contrasts matrix, so I 
am not sure what you are trying to do here. You want to do the usual:

xFit2 <- contrasts.fit(xFit, contmat)
xFit2 <- eBayes(xFit2)

Best,

Jim



> Error in log(proportion/(1 - proportion)) - log(r)/2 :
>     non-conformable arrays
> In addition: Warning messages:
> 1: In if (ntarget<  1) return(NA) :
>     the condition has length>  1 and only the first element will be used
> 2: In if (ntarget<  1) return(NA) :
>     the condition has length>  1 and only the first element will be used
> 3: In if (ntarget<  1) return(NA) :
>     the condition has length>  1 and only the first element will be used
> 4: In ebayes(fit = fit, proportion = proportion, stdev.coef.lim =
> stdev.coef.lim) :
>     Estimation of var.prior failed - set to default value
> 5: In log(proportion/(1 - proportion)) : NaNs produced
> *=========================================================================
>

-- 
James W. MacDonald, M.S.
Biostatistician
Douglas Lab
University of Michigan
Department of Human Genetics
5912 Buhl
1241 E. Catherine St.
Ann Arbor MI 48109-5618
734-615-7826
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues