[BioC] decideTests and topTable

[Gordon Smyth]

Are you clear on what an adjusted p-value is? When you used topTable() you 
used the default adjustment method which is "holm", a very conservative 
method. When you used decideTests() you used adjustment method "fdr", a 
much less conservative method. When you looked directly at the fitted model 
object, the p-value is not adjusted at all. It is not therefore not 
surprising that the p-values are different in the three cases.

Gordon

>Date: Tue, 19 Apr 2005 10:53:31 -0400
>From: Humberto Ortiz Zuazaga <humberto at hpcf.upr.edu>
>Subject: [BioC] decideTests and topTable
>To: Bioconductor at stat.math.ethz.ch
>
>I've got some Agilent microarray data for 3 different groups, and am 
>trying to
>use limma to select differentially expressed genes.
>
>Here's the steps I've taken so far:
>
> > allone <-function (qta)
>+ {
>+      1
>+ }
> > targets <- readTargets("new-cort-targets.txt")
> > RG <- read.maimages(targets$FileName,source="agilent",wt.fun=allone)
>Read 1/Cort(N-C)2.txt
>Read 2/cort(N-E)2.txt
>Read 2/cort(N-E)3.txt
>Read 5/dye swap/cort(C-N).txt
>Read 5/dye swap/Cort(E-N).txt
> > types <- readSpotTypes("spottype.txt")
> > status <- controlStatus(types, RG)
>Matching patterns for: ControlType
>Found 22393 other
>Found 913 Positive
>Found 162 Negative
>Found 21318 gene
>Setting attributes: values Color ID Name
> > RG$genes$Status <- status
> > RG <- backgroundCorrect(RG, method="none")
> > weights <- modifyWeights(RG$weights, status,
>+                          values=c("Positive","Negative"),multipliers=0)
> > MA <- normalizeWithinArrays(RG,method="loess",weights=weights)
> > MA.b <- normalizeBetweenArrays(MA,method="scale")
> > design <- modelMatrix(targets, ref="Naive")
>Found unique target names:
>  Control Enriched Naive
> > fit <- lmFit(MA.b,design,weights=weights)
> > fit.b <- eBayes(fit)
> > table <- topTable(fit.b,coef=1,number=100)
> > write.table(table,file="table.txt", sep="\t", col.names = NA)
>
> > calls.strict <- decideTests(fit.b,adjust.method="fdr")
> > write.fit(fit.b,results=calls.strict,file="fit.txt",digits=3)
>
>My question is, when I look at the top table, my best candidate is
>
>""      "Row"   "Col"   "ProbeUID"      "ControlType"   "ProbeName" 
>"GeneName"      "Description"
>"Status"        "M"     "A"     "t"     "P.Value"       "B"
>"10984" 
>52      106     10181   0       "A_51_P443387"  "AJ276707"      "Mus 
>musculus partial mRNA
>for WTAP 
>protein"       "gene"  -0.9176259      9.403058        -11.262342 
>0.2490771       -2.218647
>
>Which has an adjusted p value of 0.2490771
>
>The fit object also has a p-value column, and it is adjusted in write.fit, 
>but
>the corresponding line from the fit is:
>
>A       Control Enriched        Control Enriched        Control 
>Enriched        Control Enriched        Row     Col
>ProbeUID        ControlType     ProbeName       GeneName 
>Description     Status
>  9.40   -0.918  -0.267  -11.26  -4.02   0.00001 0.00533 
> -1      0       52      106     10181   0
>A_51_P443387    AJ276707        Mus musculus partial mRNA for WTAP 
>protein      gene
>
>The p value for contrast 1 is 0.00001.
>
>Why are the p values so different?
>
>Can I say this gene is or is not differentially expressed? Note that the
>decideTests result for contrast 1 is -1, so I understand that decideTests
>thinks it is differentially expressed. Looking at the topTable output,
>however,  makes it unlikely to be differentially expressed.
>
>--
>Humberto Ortiz Zuazaga
>Programmer-Archaeologist
>High Performance Computing facility
>University of Puerto Rico
>http://www.hpcf.upr.edu/~humberto/