Hello,

I have a design like below, and I would like to compare (1) cont1 vs. treat1, and (2) cont2 vs. treat2. It is interesting to see that the results are different by analyzing the data in the following 2 ways. For example, the topTables of cont1 vs. treat1 are different using the 2 ways. Could you please let me know what causes the difference??  Which is the correct way to go??

1. Use the design matrix with all 12 samples, i.e. model.matrix(~0+type), and set up contrasts like diff1=treat1-cont1, and diff2=treat2-cont2, and look at the result of topTable(fit, coef=1), and topTable(fit, coef=2)

samples               type       Experiment

1                            cont1       1

2                            cont1       1

3                            cont1       1

4                           treat1       1

5                           treat1       1

6                           treat1       1

7                           cont2       2

8                           cont2        2

9                           cont2        2

10                       treat2        2

11                       treat2        2

12                       treat2        2

           2. Subset the data to the first 6 samples, and use the design matrix based on the 6 samples, and look at the result of topTable(fit). For the other 6 samples in experiment 2, I do them separately.

samples               type       Experiment

1                            cont1       1

2                            cont1       1

3                            cont1       1

4                           treat1       1

5                           treat1       1

6                           treat1       1

Thank you very much !

Best regards,

Xiayu

The first approach uses the entire dataset and will give you more degrees of freedom for variance estimation. This will result in a more precise variance estimate for each gene. Ultimately, this should result in greater power for detecting DE genes across each of your contrasts.

The flipside is that the first approach is only sensible when the variance of each gene is similar between the first and second experiments. This is because limma will only estimate a single variance for each gene. If the variances of the two experiments are systematically different, the single estimated value will not accurately model the data from either experiment.

It's generally recommended to use limma on the entire dataset (i.e., your first approach). The second approach should only be considered if there is substantial evidence for different variances between experiments.