Hi,

I have RNAseq for some liver metastases which have been classified into three independent groups based upon their expression profiles.

I did pairwise differential expression analyses on these groups using both DESeq and limma, e.g., grp1 - grp2, grp1 - grp3, and grp2 - grp3. There were no significant DE following multiple testing corrections for either of the methods (not necessarily unexpected).

I then tried a GLM analysis with limma's lmFit() using the following design, e.g.,:

```
groups <- factor(c(1,1,1,2,2,2,3,3,3))
design <- model.matrix(~ 0 + groups)
contrasts <- makeContrasts(group1, group2, group3, group1 - group2, group1 - group3, group2 - group3, levels=design)
fit <- lmFit(data, design)
fit2 <- lmFit(fit, contrasts)
fit2 < eBayes(fit2)
```

When looking at results, all of the group by group contrasts (e.g., group1 - group2) are consistent with the pairwise tests showing no significantly differential expressed genes. However, when I look at the individual group coefficients, e.g.,:

```
topTable(fit2, coef="group1")
```

Nearly every gene is strongly differentially expressed. For example, for group1 nearly 90% of gene annotations showed adjusted pvalues less than 0.05. So I've come to the conclusion that I'm misinterpreting what these individual group coefficients mean, that my design is incorrect, and that I'm not reporting the right statistics. After all, how can group1 not be statistically different from group2 or group3, but strongly different from the combination of group2 + group3?

Ah, that makes sense. So essentially I'm just looking at an intercept for each group.

One more follow up if you don't mind.

In order to capture the combined group contrasts (e.g., between g1 - (g2 + g3)) I've tried creating multiple additional factors which have the appropriate group memberships:

However, I get an error from lmFit saying the coefficients for the new factors are not calculable. I assume this is because they are dependent. I'm not sure how to set up the contrasts without them though. The only time I got reasonable results was when I used one binary factor for a single comparison that included an intercept:

I'm not sure his is a correct design either. What am I overlooking with respect to model.matrix and makeContrasts?

Thanks for your help, BTW.

A very nice and clear answer indeed

Yes, that is correct. I tried:

But that gave incorrect results. As you noted before it doesn't sum to 0.

I figured it out. To do the comparison you would do something like this:

makeContrasts(group1 - (group2 + group3) / 2, ....)

Thanks again for the help.