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/

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@hpcf.upr.edu> >Subject: [BioC] decideTests and topTable >To: Bioconductor@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/