[BioC] Limma Dual Color Agilent : Mixed model or not?

[Gordon K Smyth]

Dear Guillaume,

Not all questions have brief answers!  However see the short paper

   http://www.statsci.org/smyth/pubs/ISI2005-116.pdf

which is referenced on the help page for lmscFit.  The paragraph on page 2 
starting "Replicated arrays" explains why I would recommend a standard 
log-ratio analysis for the model you have fitted.

In general, the separate channel analysis gives more DE genes than a 
standard log-ratio analysis because it makes more model assumptions (of 
common intra-spot correlation) and recovers information from the A-values 
which is otherwise ignored.  An additional reason in your case is that you 
are imputing data for the separate channel analysis, so there are more non 
NA data values being used than for the log-ratio analysis.  For the model 
you have fitted, there is actually no information in the A-values to 
recover, so I would not do it.

Your targets frame seems to include other variables (RepNum and Polarity) 
not used in the analysis, but I assume you have good reasons to ignore 
them.

Best wishes
Gordon

> Date: Tue, 29 Nov 2011 16:10:02 +0100
> From: Guillaume Meurice <guillaume.meurice at igr.fr>
> To: Bioconductor mailing list <bioconductor at r-project.org>
> Subject: [BioC] Limma Dual Color Agilent : Mixed model or not ?
> Content-Type: text/plain
>
> Dear all,
>
> I have the project with the following design (dual color, with dye swap):
>
>                        Cy5     Cy3 Status RepNum Polarity
> R1 SI31_vs_R1 CTRL+ R1_SI31 R1_CTRL      R      1        +
> R1 SI31_vs_R1 CTRL- R1_CTRL R1_SI31      R      1        -
> R2 SI31_vs_R2 CTRL+ R2_SI31 R2_CTRL      R      2        +
> R2 SI31_vs_R2 CTRL- R2_CTRL R2_SI31      R      2        -
> S1 SI31_vs_S1 CTRL+ S1_SI31 S1_CTRL      S      1        +
> S1 SI31_vs_S1 CTRL- S1_CTRL S1_SI31      S      1        -
> S3 SI31_vs_S3 CTRL+ S3_SI31 S3_CTRL      S      3        +
> S3 SI31_vs_S3 CTRL- S3_CTRL S3_SI31      S      3        -
> R3 SI31_vs_R3 CTRL+ R3_SI31 R3_CTRL      R      3        +
> R3 SI31_vs_R3 CTRL- R3_CTRL R3_SI31      R      3        -
> S2 SI31_vs_S2 CTRL+ S2_SI31 S2_CTRL      S      2        +
> S2 SI31_vs_S2 CTRL- S2_CTRL S2_SI31      S      2        -
>
>
> I want to compare R versus S.
>
> The first thing I've done is the following (where X is my processed log2(ratio) matrix):
> ==========================================================
> status = factor(target$Status)
> dye = rep(c(1,-1),times = 6)
> design = model.matrix(~0+dye+status)
> colnames(design) = c("Dye","R","S")
> fit <- lmFit(X, design)
> cont.matrix <- makeContrasts("R-S",levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
> res1 = topTable(fit2, adjust="BH", number = Inf, coef="R-S")
> ==========================================================
>
> which raise 93 transcript significantly differentially expressed.
>
>
>
>
> Then, I've tried to use the mixed model approach which raise about 4000 transcripts :
>
> ==========================================================
> MAmixed = normalizeBetweenArrays(RG,method = "Aquantile")
> # replace NA using KNN
> MnoNa = gestionNA_knn(MAmixed$M)
> AnoNa = gestionNA_knn(MAmixed$A)
>
> MAmixed$M = MnoNa
> MAmixed$A = AnoNa
>
> t2 = targetsA2C(target)
> u <- unique(t2$Status)
> f <- factor(t2$Status, levels=u)
> design <- model.matrix(~0+f)
> colnames(design) <- u
> corfit <- intraspotCorrelation(MAmixed, design)
> fit <- lmscFit(MAmixed, design, correlation=corfit$consensus)
> cont.matrix <- makeContrasts("R-S",levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
> res2 = topTable(fit2, adjust="BH", number = Inf)
> ==========================================================
>
>
>
>
> I was wondering which approach better fit the design of this project.
>
> I would also be very grateful if you could (briefly) explain why there is such a difference ?
>
> Many thanks by advance.
> --
> Guillaume

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}