I have posted this question to BioStar. http://biostar.stackexchange.com/questions/9662/ Perhaps few limma users monitor BioStar yet, so I'm cross-posting here to the Bioconductor mailing list in hopes that a few lurk here. I will re-post the contents of the question below: The study [GSE13465][1] contains experiments of toxic effects of acetaminophen on two organisms: ?- human ?- rat I only care about data from the rat cultures. (The human data is present due to submitter error.) The rat cultures took place in two media types: ?- rat hepatocytes in normal media ?- rat hepatocytes in modified media Each media condition was paired with an acetaminophen (APAP) treatment condition for each biological replicate: ?- no exposure (0 mM APAP), the control ?- exposure to 5 mM APAP ?- exposure to 10 mM APAP Expression levels were measured for each biological replicate-medium combination using the two-channel Agilent 011868 Rat Oligo Microarray G4130A ([GPL890][4]). All gene expression measurements were performed with the labels as such: ?- Channel 1 (Cy5): 0 mM APAP (control) or 5 or 10 mM APAP (treatment) ?- Channel 2 (Cy3): 0 mM APAP (control) This experiment design was repeated for three biological replicates (rats 2, 3, and 4). I would like to perform the following contrasts: ?- 5 mM APAP vs. control, standard media ?- 10 mM APAP vs. control, standard media ?- 5 mM APAP vs. control, modified media ?- 10 mM APAP vs. control, modified media **My challenge is the following: there was no single mRNA sample common across all the microarray experiments within a medium.** >From what I understand, ideally the researchers should have pooled the mRNA from each 0 mM APAP standard medium culture, then labeled this pool with Cy3 and compared the pooled mRNA to each individual culture's mRNA labeled with Cy5. (They would have repeated this design for the modified medium, as well.) Unfortunately, this was not what the researchers did. Instead, within each biological replicate, the researchers used the mRNA from the untreated sample from the same rat as a reference on the Cy3 channel for that rat's cultures in that medium, and *only* that rat's cultures. For example, in experiment GSM339438, Rat 2's cells exposed to 5 mM APAP in standard medium were compared to Rat 2's cells exposed to no APAP in standard medium. A look at a subset of my `targets` may be helpful here: ??? > targets[targets$Medium == "standard",] ????????????????????????????????????????? FileNameCy3???????????? FileNameCy5 ??? Std_APAP0.vs.Std_APAP0.2? GSM339437_Cy3_36569.txt GSM339437_Cy5_36569.txt ??? Std_APAP5.vs.Std_APAP0.2? GSM339438_Cy3_36570.txt GSM339438_Cy5_36570.txt ??? Std_APAP10.vs.Std_APAP0.2 GSM339439_Cy3_36571.txt GSM339439_Cy5_36571.txt ??? Std_APAP0.vs.Std_APAP0.3? GSM339443_Cy3_36630.txt GSM339443_Cy5_36630.txt ??? Std_APAP5.vs.Std_APAP0.3? GSM339444_Cy3_36631.txt GSM339444_Cy5_36631.txt ??? Std_APAP10.vs.Std_APAP0.3 GSM339445_Cy3_36632.txt GSM339445_Cy5_36632.txt ??? Std_APAP0.vs.Std_APAP0.4? GSM339449_Cy3_37372.txt GSM339449_Cy5_37372.txt ??? Std_APAP5.vs.Std_APAP0.4? GSM339450_Cy3_37373.txt GSM339450_Cy5_37373.txt ??? Std_APAP10.vs.Std_APAP0.4 GSM339451_Cy3_37374.txt GSM339451_Cy5_37374.txt ????????????????????????????? GEOSample????????? Organism Donor?? Medium Cy3APAP ??? Std_APAP0.vs.Std_APAP0.2? GSM339437 Rattus norvegicus???? 2 standard?????? 0 ??? Std_APAP5.vs.Std_APAP0.2? GSM339438 Rattus norvegicus???? 2 standard?????? 0 ??? Std_APAP10.vs.Std_APAP0.2 GSM339439 Rattus norvegicus???? 2 standard?????? 0 ??? Std_APAP0.vs.Std_APAP0.3? GSM339443 Rattus norvegicus???? 3 standard?????? 0 ??? Std_APAP5.vs.Std_APAP0.3? GSM339444 Rattus norvegicus???? 3 standard?????? 0 ??? Std_APAP10.vs.Std_APAP0.3 GSM339445 Rattus norvegicus???? 3 standard?????? 0 ??? Std_APAP0.vs.Std_APAP0.4? GSM339449 Rattus norvegicus???? 4 standard?????? 0 ??? Std_APAP5.vs.Std_APAP0.4? GSM339450 Rattus norvegicus???? 4 standard?????? 0 ??? Std_APAP10.vs.Std_APAP0.4 GSM339451 Rattus norvegicus???? 4 standard?????? 0 ????????????????????????????? Cy5APAP?????? Cy3??????? Cy5 ??? Std_APAP0.vs.Std_APAP0.2??????? 0 Std_APAP0? Std_APAP0 ??? Std_APAP5.vs.Std_APAP0.2??????? 5 Std_APAP0? Std_APAP5 ??? Std_APAP10.vs.Std_APAP0.2????? 10 Std_APAP0 Std_APAP10 ??? Std_APAP0.vs.Std_APAP0.3??????? 0 Std_APAP0? Std_APAP0 ??? Std_APAP5.vs.Std_APAP0.3??????? 5 Std_APAP0? Std_APAP5 ??? Std_APAP10.vs.Std_APAP0.3????? 10 Std_APAP0 Std_APAP10 ??? Std_APAP0.vs.Std_APAP0.4??????? 0 Std_APAP0? Std_APAP0 ??? Std_APAP5.vs.Std_APAP0.4??????? 5 Std_APAP0? Std_APAP5 ??? Std_APAP10.vs.Std_APAP0.4????? 10 Std_APAP0 Std_APAP10 Can I still consider the 0 mM APAP conditions as one common reference, even though each 0 mM APAP sample is actually only a common mRNA reference within the same rat ("`Donor`")? In other words, does the code below build a limma design model that is valid, or is the model invalid? ??? > std.targets <- targets[targets$Medium == "standard",] ??? > std.design <- modelMatrix(std.targets, ref="Std_APAP0") ??? > std.design ????????????????????????????? Std_APAP10 Std_APAP5 ??? Std_APAP0.vs.Std_APAP0.2?????????? 0???????? 0 ??? Std_APAP5.vs.Std_APAP0.2?????????? 0???????? 1 ??? Std_APAP10.vs.Std_APAP0.2????????? 1???????? 0 ??? Std_APAP0.vs.Std_APAP0.3?????????? 0???????? 0 ??? Std_APAP5.vs.Std_APAP0.3?????????? 0???????? 1 ??? Std_APAP10.vs.Std_APAP0.3????????? 1???????? 0 ??? Std_APAP0.vs.Std_APAP0.4?????????? 0???????? 0 ??? Std_APAP5.vs.Std_APAP0.4?????????? 0???????? 1 ??? Std_APAP10.vs.Std_APAP0.4????????? 1???????? 0 ??? > # Line below includes a dye effect in the model. ??? > # NOTE: Apparently we can't check for a dye effect, even though the ??? > # self-self hybridization of APAP0s should tell this AFAIK. Uncomment ??? > # to see this code break. ??? > #std.design <- cbind(DyeEffect=1, std.design) ??? > std.fit <- lmFit(MA.std.genes.only, std.design) ??? > std.fit <- eBayes(std.fit) ??? > # I'm going to assume the Std_APAP5 and Std_APAP10 are my contrasts ??? > # against Std_APAP0. The control mRNA samples are from *biological* replicates, but not identical across arrays. Is this sufficient? If not, what sort of design and contrast matrices should I develop? For example, should I consider this to be like a paired-samples type study (Section 8.3 in the limma User's Guide, page 43). In this case, I would create a design matrix using the factors of rat ("`Donor`") and treatment ("`Cy5APAP`"), such as follows: ??? > std.targets <- targets[targets$Medium == "standard",] ??? > std.rats <- factor(std.targets$Donor) ??? > std.rats ??? [1] 2 2 2 3 3 3 4 4 4 ??? Levels: 2 3 4 ??? > std.treat <- factor(std.targets$Cy5APAP) ??? > std.treat ??? [1] 0? 5? 10 0? 5? 10 0? 5? 10 ??? Levels: 0 5 10 ??? > std.design <- model.matrix(~std.rats+std.treat) ??? > rownames(std.design) <- rownames(std.targets) ??? > std.design ????????????????????????????? (Intercept) std.rats3 std.rats4 std.treat5 ??? Std_APAP0.vs.Std_APAP0.2??????????? 1???????? 0???????? 0????????? 0 ??? Std_APAP5.vs.Std_APAP0.2??????????? 1???????? 0???????? 0????????? 1 ??? Std_APAP10.vs.Std_APAP0.2?????????? 1???????? 0???????? 0????????? 0 ??? Std_APAP0.vs.Std_APAP0.3??????????? 1???????? 1???????? 0????????? 0 ??? Std_APAP5.vs.Std_APAP0.3??????????? 1???????? 1???????? 0????????? 1 ??? Std_APAP10.vs.Std_APAP0.3?????????? 1???????? 1???????? 0????????? 0 ??? Std_APAP0.vs.Std_APAP0.4??????????? 1???????? 0???????? 1????????? 0 ??? Std_APAP5.vs.Std_APAP0.4??????????? 1???????? 0???????? 1????????? 1 ??? Std_APAP10.vs.Std_APAP0.4?????????? 1???????? 0???????? 1????????? 0 ????????????????????????????? std.treat10 ??? Std_APAP0.vs.Std_APAP0.2??????????? 0 ??? Std_APAP5.vs.Std_APAP0.2??????????? 0 ??? Std_APAP10.vs.Std_APAP0.2?????????? 1 ??? Std_APAP0.vs.Std_APAP0.3??????????? 0 ??? Std_APAP5.vs.Std_APAP0.3??????????? 0 ??? Std_APAP10.vs.Std_APAP0.3?????????? 1 ??? Std_APAP0.vs.Std_APAP0.4??????????? 0 ??? Std_APAP5.vs.Std_APAP0.4??????????? 0 ??? Std_APAP10.vs.Std_APAP0.4?????????? 1 ??? attr(,"assign") ??? [1] 0 1 1 2 2 ??? attr(,"contrasts") ??? attr(,"contrasts")$std.rats ??? [1] "contr.treatment" ??? attr(,"contrasts")$std.treat ??? [1] "contr.treatment" ??? > std.fit <- lmFit(MA.std.genes.only, std.design) ??? > std.fit <- eBayes(std.fit) ??? > # Is this correct that I don't need a contrasts matrix; that ??? > # std.treat5 is representing the contrast ??? > # "Std_APAP5 vs. Std_APAP0?" I take it in this design, the particular rat is being considered part of the model, with Rat 2 as the baseline, and then the treatment with APAP is considered an additional factor on top of the particular rat. It's not obvious to me that this is a correct approach either, though. Any thoughts and suggestions on the appropriate design and contrast matrices would be extremely helpful. ? [1]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13465