[BioC] Understanding limma, fdr and topTable

Louisa A Rispoli/AS/EXP/UTIA larispoli at mail.ag.utk.edu
Mon Jul 7 22:35:38 CEST 2008


To all-

I am a newbie trying to analyze microarray with minimal help and thought
that I had figured out all. We have a simple task of comparing two groups
with 8 replicates on the bovine genechip. I am attempting to understand the
results that I obtain using limma and adjusting for fdr. I have tried
reading the help vignettes on p.adjust and topTablle but no closer to
understanding if the adjusted p-value represents the fdr or the q-value or
something altogether different. Based on some recent work in another lab
(by abstract) I know that I may need to use a less stringent fdr then 5%
but I am unsure in limma how to change that value (or if it is feasible at
all). Looking at the results that I obtained so far, I do not have any
differentially expressed genes with the fdr adjustment. But I could be
interpreting that wrong, since I was looking for values of the adjusted p
to be lower then 0.05 and the smallest that I see is 0.7816, If someone
could look at this and assist, any help, advice or reprimmand would be
appreciated.


Thanks

Louisa

"If we knew what we were doing, it wouldn't be called Research." - Albert
Einstein

Louisa Rispoli, Ph.D. Reproductive Physiology
Department of Animal Science
University of Tennessee, Knoxville
A105 Johnson Animal Research and Teaching Unit
1750 Alcoa Highway
Knoxville, TN 37920
phone:(865) 946-1874
fax:(865) 946-1010
email: lrispoli at utk.edu

R version 2.7.1 (2008-06-23)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
> library(affycoretools)
Loading required package: affy
Loading required package: Biobase
Loading required package: tools

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: affyio
Loading required package: preprocessCore
Loading required package: limma
Loading required package: GOstats
Loading required package: graph
Loading required package: GO.db
Loading required package: AnnotationDbi
Loading required package: DBI
Loading required package: RSQLite
Loading required package: annotate
Loading required package: xtable
Loading required package: RBGL
Loading required package: Category
Loading required package: genefilter
Loading required package: survival
Loading required package: splines
Loading required package: biomaRt
Loading required package: RCurl

Attaching package: 'biomaRt'


        The following object(s) are masked from package:annotate :

         getGO

Loading required package: gcrma
Loading required package: matchprobes
Loading required package: annaffy
Loading required package: KEGG.db

Attaching package: 'annaffy'


        The following object(s) are masked from package:RCurl :

         getURL

> pd <- read.AnnotatedDataFrame("pData.txt", sep="\t", header=TRUE,
row.names=1)
> data <- ReadAffy(phenoData=pd)
> pData(data)
               amp  trt
PolyC-1.CEL  PolyA Ctrl
PolyC-2.CEL  PolyA Ctrl
PolyC-3.CEL  PolyA Ctrl
PolyC-4.CEL  PolyA Ctrl
PolyC-5.CEL  PolyA Ctrl
PolyC-6.CEL  PolyA Ctrl
PolyC-7.CEL  PolyA Ctrl
PolyC-8.CEL  PolyA Ctrl
PolyHS-1.CEL PolyA   HS
PolyHS-2.CEL PolyA   HS
PolyHS-3.CEL PolyA   HS
PolyHS-4.CEL PolyA   HS
PolyHS-5.CEL PolyA   HS
PolyHS-6.CEL PolyA   HS
PolyHS-7.CEL PolyA   HS
PolyHS-8.CEL PolyA   HS
WTC-1.CEL       WT Ctrl
WTC-2.CEL       WT Ctrl
WTC-3.CEL       WT Ctrl
WTC-4.CEL       WT Ctrl
WTC-5.CEL       WT Ctrl
WTC-6.CEL       WT Ctrl
WTC-7.CEL       WT Ctrl
WTC-8.CEL       WT Ctrl
WTHS-1.CEL      WT   HS
WTHS-2.CEL      WT   HS
WTHS-3.CEL      WT   HS
WTHS-4.CEL      WT   HS
WTHS-5.CEL      WT   HS
WTHS-6.CEL      WT   HS
WTHS-7.CEL      WT   HS
WTHS-8.CEL      WT   HS
> Poly.rma <- rma(data[,1:16])
Background correcting
Normalizing
Calculating Expression
> pData(Poly.rma)
               amp  trt
PolyC-1.CEL  PolyA Ctrl
PolyC-2.CEL  PolyA Ctrl
PolyC-3.CEL  PolyA Ctrl
PolyC-4.CEL  PolyA Ctrl
PolyC-5.CEL  PolyA Ctrl
PolyC-6.CEL  PolyA Ctrl
PolyC-7.CEL  PolyA Ctrl
PolyC-8.CEL  PolyA Ctrl
PolyHS-1.CEL PolyA   HS
PolyHS-2.CEL PolyA   HS
PolyHS-3.CEL PolyA   HS
PolyHS-4.CEL PolyA   HS
PolyHS-5.CEL PolyA   HS
PolyHS-6.CEL PolyA   HS
PolyHS-7.CEL PolyA   HS
PolyHS-8.CEL PolyA   HS
> treatment <-c("C",
"C","C","C","C","C","C","C","HS","HS","HS","HS","HS","HS","HS","HS")
> design <- model.matrix(~factor(treatment))
> colnames(design) <- c("Ctrl", "CvsHS")
> design
   Ctrl CvsHS
1     1     0
2     1     0
3     1     0
4     1     0
5     1     0
6     1     0
7     1     0
8     1     0
9     1     1
10    1     1
11    1     1
12    1     1
13    1     1
14    1     1
15    1     1
16    1     1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$`factor(treatment)`
[1] "contr.treatment"
> options(digits=4)
> topTable(fit2, coef=2, n=25, sort.by="B", adjust="fdr")
                      ID   logFC AveExpr      t   P.Value adj.P.Val      B
16088   Bt.27852.2.S1_at -0.4764   3.621 -5.274 5.709e-05    0.7816 -1.722
12937   Bt.24859.1.A1_at  0.6632   6.940  5.152 7.374e-05    0.7816 -1.790
2853    Bt.13563.2.S1_at -0.5288   6.747 -4.703 1.922e-04    0.7816 -2.061
2474    Bt.13162.1.S1_at -0.6424   5.964 -4.574 2.534e-04    0.7816 -2.142
10402   Bt.22107.1.S1_at -0.3261  10.358 -4.475 3.143e-04    0.7816 -2.207
10654   Bt.22355.1.S1_at -0.3533   8.650 -4.316 4.439e-04    0.7816 -2.312
8883    Bt.20584.1.S1_at -0.3671   8.106 -4.225 5.417e-04    0.7816 -2.374
2851    Bt.13562.1.S1_at  0.6901   8.865  4.190 5.850e-04    0.7816 -2.399
5772     Bt.1754.1.S1_at -0.5986   9.348 -4.187 5.885e-04    0.7816 -2.400
19715    Bt.4440.1.A1_at -0.6881   2.526 -4.160 6.242e-04    0.7816 -2.419
3163    Bt.13933.1.S1_at  0.4519   7.798  4.138 6.549e-04    0.7816 -2.434
14921   Bt.26765.1.S1_at -0.3673   7.846 -4.110 6.962e-04    0.7816 -2.454
10751   Bt.22445.1.S1_at -0.3222   8.622 -4.108 7.002e-04    0.7816 -2.456
2906    Bt.13629.1.A1_at  0.2830   2.294  4.095 7.196e-04    0.7816 -2.464
5721     Bt.1749.1.S1_at  0.6030   3.839  4.072 7.578e-04    0.7816 -2.481
21192    Bt.6078.1.S1_at -0.4385   8.564 -4.054 7.879e-04    0.7816 -2.493
800     Bt.11145.1.S1_at -0.5060   2.461 -4.030 8.312e-04    0.7816 -2.511
27    AFFX-Bt-ECOLOXL_at  0.3410   1.290  4.022 8.443e-04    0.7816 -2.516
22515    Bt.7980.1.S1_at -0.3734   7.622 -3.983 9.196e-04    0.7816 -2.543
4609    Bt.16378.1.A1_at -0.3775   5.071 -3.964 9.607e-04    0.7816 -2.557
13765   Bt.25669.1.S1_at -0.5523   8.197 -3.948 9.933e-04    0.7816 -2.568
12023   Bt.24033.1.A1_at -0.4962   6.256 -3.917 1.065e-03    0.7816 -2.591
20843    Bt.5578.1.S1_at -0.4133   7.904 -3.913 1.073e-03    0.7816 -2.594
17021   Bt.28620.1.S1_at  0.4720   5.181  3.899 1.106e-03    0.7816 -2.603
23922    Bt.9791.1.S1_at -0.3420   9.963 -3.893 1.122e-03    0.7816 -2.608

> sessionInfo()
R version 2.7.1 (2008-06-23)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United

States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] splines   tools     stats     graphics  grDevices utils     datasets
methods   base

other attached packages:
 [1] bovinecdf_2.2.0      affycoretools_1.12.0 annaffy_1.12.1
KEGG.db_2.2.0        gcrma_2.12.1         matchprobes_1.12.0
 [7] biomaRt_1.14.0       RCurl_0.9-3          GOstats_2.6.0
Category_2.6.0       genefilter_1.20.0    survival_2.34-1
[13] RBGL_1.16.0          annotate_1.18.0      xtable_1.5-2
GO.db_2.2.0          AnnotationDbi_1.2.2  RSQLite_0.6-9
[19] DBI_0.2-4            graph_1.18.1         limma_2.14.5
affy_1.18.2          preprocessCore_1.2.0 affyio_1.8.0
[25] Biobase_2.0.1

loaded via a namespace (and not attached):
[1] cluster_1.11.11 XML_1.95-3



More information about the Bioconductor mailing list