[altirriba@hotmail.com: [BioC] Design in factDesign] (fwd)

[Denise Scholtens]

Hello Jordi,

My comments to your questions are below.  I hope this helps.  -Denise

__________________________________________________________________________
Denise Scholtens
Department of Biostatistics
Harvard School of Public Health
dscholte at hsph.harvard.edu

Hi all!
I’ve been using RMA and LIMMA to analyse my data and I am currently trying
to analyse it with the package factDesign. My design is a 2x2 factorial
design with 4 groups: diabetic treated, diabetic untreated, health treated
and health untreated with 3 biological replicates in each group. I want to
know what genes are differentially expressed due to diabetes,  to the
treatment and to the combination of both (diabetes + treatment).
My phenoData is:
>pData(eset)
         DIABETES    TREATMENT
DNT1     TRUE       FALSE
DNT2     TRUE       FALSE
DNT3     TRUE       FALSE
DT1        TRUE       TRUE
DT2        TRUE       TRUE
DT3        TRUE       TRUE
SNT1      FALSE     FALSE
SNT2      FALSE     FALSE
SNT3      FALSE     FALSE
ST1         FALSE     TRUE
ST2         FALSE     TRUE
ST3         FALSE     TRUE

Are these commands correct to get the results listed below? Where are the
errors?
>lm.full<-function(y) lm(y ~ DIABETES + TREATMENT + DIABETES * TREATMENT)
>lm.diabetes<-function(y) lm(y ~ DIABETES)
>lm.treatment<-function(y) lm(y ~ TREATMENT)
>lm.diabetestreatment<-function(y) lm(y ~ DIABETES + TREATMENT)
>lm.f<-esApply(eset, 1, lm.full)
>lm.d<-esApply(eset, 1, lm.diabetes)
>lm.t<-esApply(eset, 1, lm.treatment)
>lm.dt<-esApply(eset, 1, lm.diabetestreatment)

#####
# Yes, these commands look correct for making the linear models and
# running them for the exprSet called eset.
######

## To get the genes characteristics of the treatment:
>Fpvals<-rep(0, length(lm.f))
>for (i in 1:length(lm.f)) {Fpvals[i]<-anova(lm.d[[i]], lm.f[[i]])$P[2]}
>Fsub<-which(Fpvals<0.01)
>eset.Fsub<-eset[Fsub]
>lm.f.Fsub<-lm.f[Fsub]
>betaNames<-names(lm.f[[1]] [["coefficients"]])
>lambda<-par2lambda(betaNames, c("TREATMENTTRUE"), c(1)) ## I get the same
>genes if I write : > lambda2 <- par2lambda (betaNames,
>list(c("TREATMENTTRUE" , "DIABETESTRUE:TREATMENTTRUE")),list( c(1,1)))
>mainTR<-function(x) contrastTest(x,lambda,p=0.1) [[1]]
>mainTRgenes<-sapply(lm.f.Fsub, FUN=mainES)

#####
# I think the problem is the use of mainES rather than mainTR in the last
# sapply. mainES is a function that is defined in the factDesign vignette
# - your own function should be used here instead.  If you define the
# function differently for different contrasts, my guess is you will see
# different gene lists for the lambda and lambda2 defined above.
#####


## To get the genes characteristics of the diabetes:
>for (i in 1:length(lm.f)) {Fpvals[i]<-anova(lm.t[[i]], lm.f[[i]])$P[2]}
>Fsub<-which(Fpvals<0.01)
>eset.Fsub<-eset[Fsub]
>lm.f.Fsub<-lm.f[Fsub]
>betaNames<-names(lm.f[[1]] [["coefficients"]])
>lambda<-par2lambda(betaNames, c("DIABETESTRUE"), c(1)) ## I get also the
>same genes if I consider the intersection DIABETESTRUE:TREATMENTTRUE.
>mainDI<-function(x) contrastTest(x,lambda,p=0.1) [[1]]
>mainDIgenes<-sapply(lm.f.Fsub, FUN=mainES)

#####
# See above comments.
#####

## To get the genes characteristics of the diabetes+treatment:
>for (i in 1:length(lm.f)) {Fpvals[i]<-anova(lm.dt[[i]], lm.f[[i]])$P[2]}
>Fsub<-which(Fpvals<0.01)
>eset.Fsub<-eset[Fsub]
>lm.f.Fsub<-lm.f[Fsub]
>  betaNames<-names(lm.f[[1]] [["coefficients"]])
>lambda<-par2lambda(betaNames, c("DIABETESTRUE:TREATMENTTRUE"), c(1))
>mainDT<-function(x) contrastTest(x,lambda,p=0.1) [[1]]
>mainDTgenes<-sapply(lm.f.Fsub, FUN=mainES) ## I don’t get any “fail to
>reject” gene.

#####
# Again, I think changing mainES to mainDT will do the trick.
#####


When I get the “rejected” and the “failed to reject” genes, can I classify
them by their Fvalues? How?

#####
# Currently, the contrastTest function only returns the contrast estimate
# (cEst), the pvalue from the F-test (pvalue), and a statement of either
# "REJECT" or "FAIL TO REJECT" based on the p-value cutoff you specify.
# This can be changed to return the F-value as well, and I'm happy to put
# this change into the package.  Then you can use the Fvalues for whatever
# you would like.
#
# One thing to consider if you are going to use p-values from the F tests
# to select genes - you will want to corrent for multiple testing.  The
# multtest package is very useful for this.
######


Thank you very much for your comments and help.
Yours sincerely,

Jordi Altirriba
IDIBAPS-Hospital Clinic (Barcelona, Spain)

_________________________________________________________________
Déjanos tu CV y recibe ofertas de trabajo en tu buzón. Multiplica tus
oportunidades con MSN Empleo. http://www.msn.es/Empleo/

_______________________________________________
Bioconductor mailing list
Bioconductor at stat.math.ethz.ch
https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor