Entering edit mode
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}}