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.,:
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?