Dear Sean,

Thanks for replying.

It really worked the way you suggested.

I have another question related to logFC

I am getting the table
               fA...fB   fA...fC      fA...fD  AveExpr        F      P.Value
35463 -7.708127 -7.938576  0.599991797 4.975526 543.0372 6.193852e-13
18194 -9.328894 -8.614855  0.071930586 7.696517 482.1491 1.244436e-12
23024 -8.227234 -6.507593  0.049250068 7.259353 425.1816 2.600162e-12
16492 -8.456246 -7.544640 -0.990340311 6.837342 404.4312 3.485316e-12

when I try to take any specific coefficient  using

topTable(fit2,number=10,coef="fA-fB")
        fA...fB   fA...fC      fA...fD  AveExpr        F      P.Value
35463 -7.708127 -7.938576  0.599991797 4.975526 543.0372 6.193852e-13
18194 -9.328894 -8.614855  0.071930586 7.696517 482.1491 1.244436e-12
23024 -8.227234 -6.507593  0.049250068 7.259353 425.1816 2.600162e-12
16492 -8.456246 -7.544640 -0.990340311 6.837342 404.4312 3.485316e-12

I get this error

Error in toptable(fit = fit[c("coefficients", "stdev.unscaled")], coef =
coef,  :
  subscript out of bounds

Does that mean fA...fB is the logFC that i can consider directly ??

Is there any way that i can sort out the topTable error to figure out logFC
for specific coeff ??

Regards

On Wed, Sep 10, 2008 at 12:21 AM, Sean Davis <sdavis2@mail.nih.gov> wrote:

> On Tue, Sep 9, 2008 at 1:31 PM, hemant ritturaj
> <ritturajhemant@gmail.com> wrote:
> > Dear Sean,
> >
> > My apologies for the errror.
> >
> > Here I am giving the complete error which i am getting
> >
> >> gse <- getGEO('gse6901')[[1]]
> >> exprs(gse)
> >> log2(exprs(gse))
> >> set <- log2(exprs(gse))
> >> targets <- c('A','A','A','B','B','B','C','C','C','D','D','D')
> >> f <- factor(targets,levels=c('A','B','C','D'))
> >> design <- model.matrix(~0+f)
> >> contrast.matrix <- makeContrasts(fA-fB,fA-fC,fA-fD,levels=design)
> >> fit <- lmFit(set,design,method="robust")
> >> fit2 <- contrasts.fit(fit,contrast.matrix)
> >> fit2 <- ebayes(fit2)
> >
> >> names(fit2)
> > [1] "df.prior"  "s2.prior"  "s2.post"   "t"         "p.value"
> "var.prior"
> > [7] "lods"
> >
> > topTable(fit2,n=10)
> > Error in if (ncol(fit) > 1) return(topTableF(fit, number = number,
> genelist
> > = genelist,  :
> >  argument is of length zero
>
> This works for me:
>
> gse <- getGEO('gse6901')[[1]]
> exprs(gse) <- log2(exprs(gse)+1) # avoid log2(0) since we are not
> normalizing
> targets <- c('A','A','A','B','B','B','C','C','C','D','D','D')
> f <- factor(targets,levels=c('A','B','C','D'))
> design <- model.matrix(~0+f)
> contrast.matrix <- makeContrasts(fA-fB,fA-fC,fA-fD,levels=design)
> fit <- lmFit(gse,design) # you can use an ExpressionSet just fine in limma
> fit2 <- contrasts.fit(fit,contrast.matrix)
> fit2 <- eBayes(fit2)
> names(fit2)
>
>  [1] "coefficients"     "rank"             "assign"           "qr"
>  [5] "df.residual"      "sigma"            "cov.coefficients"
> "stdev.unscaled"
>  [9] "genes"            "Amean"            "method"           "design"
> [13] "contrasts"        "df.prior"         "s2.prior"         "var.prior"
> [17] "proportion"       "s2.post"          "t"                "p.value"
> [21] "lods"             "F"                "F.p.value"
>
> Again, start with a clean workspace, load GEOquery and limma, then
> paste the above commands into the window.  Let us know if you have
> further questions.
>



-- 
Hemant Ritturaj Kushwaha
PhD Student
Center for Computational Biology and Bioinformatics
Jawaharlal Nehru University
New Delhi-110067
Phone: 09868801604

	[[alternative HTML version deleted]]

