[BioC] (no subject)

[Jordi Altirriba Gutiérrez]

Dear Bioconductor users,

I have found part of the answer to my question in a previous post:
https://stat.ethz.ch/pipermail/bioconductor/2008-July/023593.html
Sorry…

I would like to suggest that the adjusted p values which appear after this code:
 
> 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

Would be the one obtained with this other code:
 
> p.adjust(fit3$p.value, method="BH")
 
So adj.P.val: 0.0049664 instead of adj.P.Val: 1.621622e-08.
 
Many thanks,
Yours faithfully,
 
Jordi Altirriba
Hospital Clinic, Barcelona, Spain
 

Jordi Altirriba Gutiérrez wrote:

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