I have a question regarding limma-driven microarray analysis.
I have a total of 24 arrays. My experimental design consists in two fixed factors, i.e. "matrix" (F or G) and "dim" (2 or 3), and what I consider a Random Factor, i.e. "donor". There are six donors, so that each donor has provided a matrix * dim treatment. That is: 4 possible combinations matrix * dim * 6 donors = 24 arrays.
I run the following code:
eset <-new ("ExpressionSet", exprs = as.matrix(data.norm)) treat <- factor (substr(colnames(data.norm),1,2)) design <- model.matrix (~ 0 + treat) colnames(design) <- c("F2","F3","G2","G3") corfit <- duplicateCorrelation (eset, design, block = factor(substr(colnames(data.norm),4,6))) corfit$consensus fit <- lmFit(eset, design, block = factor (substr(colnames(data.norm),4,6)), correlation = corfit$consensus) cont.matrix <- makeContrasts( F3.F2 = F3 - F2, G3.G2 = G3 - G2, G2.F2 = G2 - F2, G3.F3 = G3 - F3, inter = (G3 - G2) - (F3 - F2), levels = design ) fit2 <- contrasts.fit (fit, cont.matrix) eb <- eBayes(fit2)
My main question is: is it right to run such code? I say this because on page 50 of the limma user's guide, it says regarding multi level designs that "Since we need to make comparaisons BOTH within AND between subjects, it is necessary to treat Patient as a Random Effect" (capital letters are mine). This is because in the example given in the document, patients were either "Normal" or "Diseased", but not both. In my design, however, a single donor covers all fixed effect combinations. I still consider my donors as random, but I wonder if the limma code I used is adapted.
A subsidiary question would deal with the fact that when running lmer models on every gene some of them present singular fits, that is over fitting. This is not visible when running the code shown above. I am more used to run lme or lmer models out of the transcriptomics studies, where thousands of models, i.e. one per gene, are fitted.
Thanks for your help.