Dear Bioconductor users, We are submitting our results to a journal and the reviewers ask us to reanalyze our data showing which p adjusted values have the genes that we selected. Our data have four different experimental groups (GEO number GDS1883), with three arrays for each group: -Healthy untreated animals: HU1, HU2, HU3 -Diabetic untreated animals: DU1, DU2, DU3 -Healthy treated animals: HT1, HT2, HT3 -Diabetic treated animals: DT1, DT2, DT3 In the past (2004) we selected those genes with a B value higher than 0. Now we are interested in those genes with a p value lower than 0.05. We have analyzed our data with the limma package, considering this model matrix: model.matrix( ~ DIABETES*TREATMENT) The detailed code is below. I have two concerns: 1.- I have analyzed the data, correcting for multiple testing with two methods, "separate" and "global". I have noticed that the "raw p value" is different in both approximations. Is it correct? I thought that only the adjusted p value would be modified. 2.- In the correction for multiple testing with method global, with an adjusted p value lower than 0.05, the last gene which has been selected (for example in the coefficient of treatment) has an adjusted p value of 3.27 E-04, which is the same value than the raw p value (the other coefficients have a similar behavior), but I would be waiting an adjusted p value close to 0.05. Am I doing something wrong? Many thanks for your time and suggestions. Yours faithfully, Jordi Altirriba Hospital Clinic, Barcelona, Spain ##Code ##Read and transform data library(affy) library(limma) Data<-ReadAffy() eset<-rma(Data) ##Phenodata sink(?pData.txt?) data.frame(DIABETES= c(rep("TRUE",6), rep("FALSE",6)), TREATMENT = c(rep("FALSE",3), rep("TRUE",3),rep("FALSE",3), rep("TRUE",3)), row.names= sampleNames(Data)) sink() pd1<-read.table("pData.txt") pData(eset)<-pd1 ##Limma analysis design <- model.matrix( ~ DIABETES*TREATMENT, data=pData(eset)) contrast.matrix<- cbind(DIABETESTRUE=c(0,1,0,0),TREATMENTTRUE=c(0,0,1,0), DIABETESTRUE.TREATMENTTRUE=c(0,0,0,1)) fit<-lmFit(eset,design) fit2<-contrasts.fit(fit,contrast.matrix) fit3<-eBayes(fit2) a) Analysis considering multiple testing with the method separate > topTable(fit3,number=1,coef="TREATMENTTRUE",adjust="BH") ID: 1388229_a_at logFC: -1.691638 AveExpr: 4.741326 t: -8.354672 P.Value: 1.194499e-06 adj.P.Val: 0.0190200 B: -0.8381138 b) Analysis considering multiple testing with the method global > results <- decideTests(fit3,method="global") > a <- as.logical(results[,"TREATMENTTRUE"]) > topTable(fit3[a,],n=1) ID: 1 1370027_a_at TREATMENTTRUE: -3.111467 DIABETESTRUE: -1.432111 DIABETESTRUE.TREATMENTTRUE: 1.548180 AveExpr: 7.242719 F: 94.78192 P.Value: 3.243243e-09 adj.P.Val: 1.621622e-08 > sessionInfo() R version 2.8.1 RC (2008-12-14 r47200) 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] tools stats graphics grDevices utils datasets methods base other attached packages: [1] rae230acdf_2.3.0 limma_2.16.4 affy_1.20.2 Biobase_2.2.2 loaded via a namespace (and not attached): [1] affyio_1.10.1 preprocessCore_1.4.0