[BioC] Limma linear models

Wu, Di dwu at fas.harvard.edu
Wed Apr 3 21:09:57 CEST 2013


Hi Joseph,

It's correct that "the fold changes would be the "Coef" column, not the A column, which is the average across all of the arrays". These are in log2 scale.

Regarding other annotation of probe IDs, I will firstly check "names(fit)".
If there is  fit$genes or fit$genes$IDs, and the anno file has the columns of IDs, GeneNames, and ReSeq, code as the following.

m<-match(fit$genes$IDs, anno$IDs)

fit$genes$symb<- GeneNames[m]
fit$genes$refseq<- ReSeq[m]
....

data.frame(fit) [1:2,]

# this dataFrame has all anno information you have added to fit$genes.

Hope this help.
Cheers,
Di

----
Di Wu
Postdoctoral fellow
Harvard University, Statistics Department
Harvard Medical School
Science Center, 1 Oxford Street, Cambridge, MA 02138-2901 USA

________________________________________
From: bioconductor-bounces at r-project.org [bioconductor-bounces at r-project.org] on behalf of Cornish, Joseph (NIH/NIAID) [F] [joseph.cornish at nih.gov]
Sent: Wednesday, April 03, 2013 1:35 PM
To: bioconductor at r-project.org
Subject: [BioC] Limma linear models

I had a few quick questions about the limma models. The first is when writing these to a file, the fold changes would be the "Coef" column, not the A column, which is the average across all of the arrays? The second is when I write my bayes model to a file, it only attaches the probe ID, I wanted to see if there was a way to attach the remaining gene annotation but the limma documentation isn't clear on how to accomplish that.

Here is how I am currently performing the fitting:
    contmat   <- makeContrasts(d-s,
                               levels = design)
    fit1      <- lmFit(x[[1]]$avg$E, design)
    fit2      <- eBayes(contrasts.fit(fit1, contmat))
    diff      <- topTable(fit2,
                          coef     = 1,
                          adjust   = "BH",
                          genelist = x[[1]]$avg$genes )
    res       <- decideTests(fit2)
    comp      <- vennDiagram(res)
    write.fit(fit1,
              results  = res,
              sep      = sep,
              file     = paste(c(x[[1]]$bg, x[[1]]$norm, "lmfit.csv"),
              collapse = ""))

    write.fit(fit2,
              results  = res,
              sep      = sep,
              file     = paste(c(x[[1]]$bg, x[[1]]$norm, "be.csv"),
              collapse = ""))

where x[[1]]$avg is the corrected, normalized, replicate averaged array.

Here is some example output from write.fit(fit2...):

A

Coef

t

p.value

F

F.p.value

d - s

10.99

5.352

6.82

6.00E-04

46.55

6.00E-04

0

12.48

5.347

14.13

1.00E-05

199.78

1.00E-05

1

12.84

4.775

11.99

3.00E-05

143.65

3.00E-05

1

13

4.646

10.98

5.00E-05

120.49

5.00E-05

1



Joseph Cornish
Post-Bac IRTA
NIAID EVPS


        [[alternative HTML version deleted]]

_______________________________________________
Bioconductor mailing list
Bioconductor at r-project.org
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor



More information about the Bioconductor mailing list