Dear all,

I'm trying to use limma in order to analyze a large number of Affymetrix microarrays. The experimental setup includes 2-3 replicates, 2-4 different time points and 1-5 different dosage levels.

So a typical targets file could look like the following where the time and dose can be non-integer numbers in reality:

> my.targets = data.frame(group = c(rep("untreated", times = 6),
                                  rep("treated", times = 18)),
                        time = rep(1:3, each = 2, times = 4),
                        dose = rep(0:3, each = 6),
                        replicate = rep(c("A", "B"), times = 12))

> my.targets
       group time dose replicate
1  untreated    1    0         A
2  untreated    1    0         B
3  untreated    2    0         A
4  untreated    2    0         B
5  untreated    3    0         A
6  untreated    3    0         B
7    treated    1    1         A
8    treated    1    1         B
9    treated    2    1         A
10   treated    2    1         B
11   treated    3    1         A
12   treated    3    1         B
13   treated    1    2         A
14   treated    1    2         B
15   treated    2    2         A
16   treated    2    2         B
17   treated    3    2         A
18   treated    3    2         B
19   treated    1    3         A
20   treated    1    3         B
21   treated    2    3         A
22   treated    2    3         B
23   treated    3    3         A
24   treated    3    3         B

as you can see the "control" group is synonymous with no dosage.

I would like to evaluate the overall genetic response to treatment and time. So usually I would accomplish this with:

> my.design = model.matrix(~ time * dose, data = my.targets)

since I'm not sure whether the response is to dosage, time or their interaction.

However, I don't know what the baseline expression would be with this design and whether limma is aware of the within replicates variability. So is it enough to use the following design:

> my.design = model.matrix(~ group + time * dose, data = my.targets)

with which limma will then use the untreated samples to establish a baseline expression? Or do I need to explicitly include the replicates as well?

Thank you for your insights.

How you specify a design is dependent on what you are trying to show. If you want to be able to say that the expression of gene X increases some amount as the dosage increases by some amount, then you have to fit a linear model with dosage as a continuous variable. As a related point, if you have many different dosages and times, then collapsing to groups is probably not what you want to do.

If you condense your dose and time values into 'buckets' then you are losing information as well as degrees of freedom. And if you only have two replicates per time/dose combination, then you almost surely don't want to collapse to factors. But this question has much more to do with 'how should I analyze my data' rather than 'how do I get limma to analyze my data in the way I have chosen?' This support site is primarily dedicated to the latter, and we as a rule at least attempt to refrain from answering the former, as that requires much more knowledge of the experiment and its goals. And if you are asking random people on the intertubes how you should be analyzing your data, you will likely get what you pay for.