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