I have an RNA-seq data set (8 treatments, 2 replicates) that I would like to process by voom transformation for analysis in limma. I also would like to include sample-specific quality weights, so I know using the function voomWithQualityWeights() is the way to go.

However, since I have only 2 replicates per group/treatment, I wonder whether it would be best to provide as argument for design  the actual design matrix of the experiment, or rather (still) the default (intercept only) design matrix.

The reason I am asking is based on this statement in the limma user guide (for the related function arrayWeights()) [section 14.4 "When to Use Array Weights", page 68]:

For example, in a two-group comparison with just 2 replicates in each group, the array weights should be estimated with the default (intercept) design matrix, otherwise each array is compared only to its partner rather than to the other 3 arrays.

Does this also apply for other designs with only 2 replicates per group? I noticed that when providing the design matrix of the experiment the sample-specific weights for each duplo sample are identical, whereas when providing the intercept design matrix the sample-weights are (indeed) different for all samples.

What would would be the best approach? Any insights on this would be appreciated!

> design
   Grp1 Grp2 Grp3 Grp4 Grp5 Grp6 Grp7 Grp8
1     0    0    0    0    1    0    0    0
2     0    0    0    0    1    0    0    0
3     0    0    0    0    0    0    1    0
4     0    0    0    0    0    0    1    0
5     0    0    0    0    0    0    0    1
6     0    0    0    0    0    0    0    1
7     0    1    0    0    0    0    0    0
8     0    1    0    0    0    0    0    0
9     0    0    0    1    0    0    0    0
10    0    0    0    1    0    0    0    0
11    1    0    0    0    0    0    0    0
12    1    0    0    0    0    0    0    0
13    0    0    1    0    0    0    0    0
14    0    0    1    0    0    0    0    0
15    0    0    0    0    0    1    0    0
16    0    0    0    0    0    1    0    0
attr(,"assign")
[1] 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1] "contr.treatment"

>

Yes, this advice also applies when using voomWithQualityWeights on RNA-seq data.

Using the actual design matrix will work best if you have 3 or more replicates per group with this model.

If you have fewer than this, an intercept model is the best you can do.

Just to refine Matt's answer: It isn't wrong to apply voomWithQualityWeights() or arrayWeights() with just n=2 in each group, but the result will be to apply a weight to each group rather than to each individual sample. It's easy to understand why it has to be like this -- when you have just two samples in a group, it isn't possible to tell which of them is the outlier (if any)!