I am analyzing a microarray experiment. The design of my study is as follows:
Each mouse has 3 measurements (baseline, 14 days, 30 days) There are NO replicates (technical or biological). I have 28 samples in total that corresponds to 10 mice and there are 3 treatments. So per treatment I have 3 mice.
Group Treatment MiceID Status GR1 TrtA 1 baseline GR1 TrtA 1 Day14 GR1 TrtA 1 Day30 GR3 Combo 2 baseline GR3 Combo 2 Day14 GR3 Combo 2 Day30 GR3 Combo 3 baseline GR3 Combo 3 Day30 GR2 TrtB 4 baseline GR2 TrtB 4 Day14 GR2 TrtB 4 Day30 GR2 TrtB 5 baseline GR2 TrtB 5 Day14 GR2 TrtB 6 baseline GR2 TrtB 6 Day14 GR2 TrtB 6 Day30 GR1 TrtA 7 baseline GR1 TrtA 7 Day14 GR1 TrtA 7 Day30 GR3 Combo 8 baseline GR3 Combo 8 Day14 GR3 Combo 8 Day30 GR3 Combo 9 baseline GR3 Combo 9 Day14 GR3 Combo 9 Day30 GR1 TrtA 10 baseline GR1 TrtA 10 Day14 GR1 TrtA 10 Day30
Goal: The question is what are the differentially expressed genes for day 14 vs baseline and day30 vs day14 per treatment?
To account for the correlation within each mice, I used the following Limma model and my corfit$consensus = -0.0068
1) Given this is a negative and nearly zero corfit correlation, should I be including correlation in my lmFit model? If not, then how do I account for the correlation per mice?
2) Does it make sense to run two Limma, one comparing baseline vs day14 and another comparing day30 vs day 14 samples
LIMMA R Code: annot : sampleInfo shown above, data.log= log2(gene expression signals)
Treat <- factor(paste(annot$Treatment,annot$Status,sep="."))
design <- model.matrix(~0+Treat)
colnames(design) <- c("TrtA.baseline", "TrtA.Day30", "TrtA.Day14", "Combo.baseline", "Combo.Day30", "Combo.Day14", "TrtB.baseline", "TrtB.Day30", "TrtB.Day14")
corfit <- duplicateCorrelation(data.log, design, block=annot$MiceID)corfit$consensus
fit1 <- lmFit(data.log,design,block=annot$miceID,correlation=corfit$consensus)
Differences <- makeContrasts(TrtALatevsTrtAwk2= TrtA.Day30-TrtA.Day14, ComboLatevsCombowk2=Combo.Day30-Combo.Day14, TrtBLate=TrtB.Day30-TrtB.Day14, TrtAwk2vsTrtA0= TrtA.Day14-TrtA.baseline, Combowk2vsCombo.0=Combo.Day14-Combo.baseline, TrtB=TrtB.Day14-TrtB.baseline,levels=design)
fit2 <- contrasts.fit(fit1,Differences) fit <- eBayes(fit2)
Thanks Aaron. Let me run this model and update you on the outcome.