[BioC] Getting different results with 2 models for factorial designs with LIMMA

[Garcia Orellana,Miriam]

Dear ALL (special request to Dr. G. Smith):
Just to give an additional information to my first email. NOW I AM SURE both of the models are considering as references those that I was expecting. Also both models for each of my 5 contrasts are giving me the top table with the same numerical values for: AveExp, t, Pvalue, adjPvalue and B. However the logFC for contrasts 1, 2 and 3 (detail of contrasts is lines below) in model B is exactly half of that in model A, while the logFC for contrasts 4 and 5 in model B is exactly one fourth of that one in model A.
How that can be possible if all other values are the same? and so what should I follow?

Thanks so much.
Miriam
________________________________________
From: bioconductor-bounces at r-project.org [bioconductor-bounces at r-project.org] on behalf of Garcia Orellana,Miriam [mgarciao at ufl.edu]
Sent: Saturday, June 23, 2012 9:29 PM
To: bioconductor at r-project.org
Subject: [BioC] Getting different results with 2 models for factorial designs with LIMMA

Dear Users:
I hope someone can help me to understand why the two models I analyzed for my data are giving me different outputs regarding differentially expressed gene. For example the MODEL A for the MR effect give me 47 up- and 84 down-regulated genes (adjPvalue <0.05 rawFC=1.5), while the model B give me only 18 up- and only 2 down-regulated genes under same cut-offs. So additional to the big difference in number of genes also the proportion of UP /DOWN in model A is lower than in MODEL B. So also I am wonder how I can be sure the program is using the right factors as the reference to calculate the logFC. My goal is to have as references the factors as indicated below
Briefly me data is a factorial design of 3 dam diets (DD: CTL, SFA, EFA) and 2 milk replacers (MR: LLA, HLA), I have three replicates for each of the interaction factors, then a total of 18 arrays. The data was filtered for informative/noninformative probes and plotted for array quality. So from a initial of 24118 bovine probes I endup with 8026 probes. My interest is to compare:

1.       Feeding FAT at prepartum= (SFA +EFA) vs CTL, with CTL as ref

2.       Feeding EFA prepartum = EFA vs SFA, with SFA as ref

3.       Feeding MR to calves= HLA vs LLA, with LLA as reference

4.       Interaction of feeding FAT by MR: (SFA +EFA) vs CTL by MR, with (SFA+EFA) vs CTL by LLA as ref

5.       Interaction of feeding EFA by MR:  EFA vs SFA by MR, with EFA vs SFA by LLA as ref

MODEL A (I created that with the guide of the LIMMA user guide for a factorial design:
TS <- paste(phenoDiet$DD, phenoDiet$MR, sep=".")
TS
TS <- factor(TS, levels=c("Ctl.LLA", "Ctl.HLA","SFA.LLA","SFA.HLA","EFA.LLA", "EFA.HLA"))
design <- model.matrix(~0+TS)
colnames(design) <- levels(TS)
fit <- lmFit(eset2, design, method="robust", maxit=1000)
efit <- eBayes(fit)
#Contrast results
MatContrast=makeContrasts(FAT=(SFA.LLA + SFA.HLA + EFA.LLA + EFA.HLA)/4 - (Ctl.LLA + Ctl.HLA)/2,
                                FA=(EFA.LLA + EFA.HLA)/2 - (SFA.LLA + SFA.HLA)/2,
                                MR=(EFA.HLA+SFA.HLA+Ctl.HLA)/3 - (EFA.LLA+SFA.LLA+Ctl.LLA)/3,
                                FATbyMR=((EFA.HLA+SFA.HLA)/2 - Ctl.HLA) - ((EFA.LLA+SFA.LLA)/2-Ctl.LLA),
                                FAbyMR=( EFA.HLA-SFA.HLA)-(EFA.LLA - SFA.LLA),
                                levels=design)
fitMat<-contrasts.fit(fit,MatContrast)
Contrast.eBayes=eBayes(fitMat)

MODEL B (this model was kindly provided by Dr G. Smith):
DD <-factor(phenoDie$DD, levels = c("Ctl", "SFA", "EFA"))
MR <-factor(phenoDie$MR, levels = c("LLA", "HLA"))
contrasts (DD) <- cbind (SFAEFAvsCtl=c(-2,1,1),EFAvsSFA=c(0,-1,1))
contrasts (MR) <- c(-1,1)
design <-model.matrix (~DD*MR)
design
fit <- lmFit (eset2, design, method="robust",maxit=1000)
efit <- eBayes(fit)

Thanks  in advance,
Miriam


        [[alternative HTML version deleted]]

_______________________________________________
Bioconductor mailing list
Bioconductor at r-project.org
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor