Adjusted p values in limma
0
0
Entering edit mode
@jordi-altirriba-gutierrez-682
Last seen 3.5 years ago
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
limma limma • 2.1k views

Login before adding your answer.

Traffic: 369 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6