[BioC] Fold change GEO data

[Mark Robinson]

Hemant,

The reason you are getting that error is because you need to match
topTable's 'coef' argument with the column names of the your
design/contrast matrix.  Guessing from your output, probably you need to
say coef="fA...fB" instead of coef="fA-fB", and so on for the others. 
Alternatively, you can specify coef=1 for the first column, coef=2, etc.

Cheers,
Mark


> 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 at mail.nih.gov> wrote:
>
>> On Tue, Sep 9, 2008 at 1:31 PM, hemant ritturaj
>> <ritturajhemant at 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]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>