LIMMA RNA-seq difference between coefficient and contrasts approach

Question

[gorenjakm]

Dear all, I wanted to adjust for continuous covariates using limma. I tried 2 approaches which so far work for 2 group designs as explained in limma user guide on page 40. However, if I add a single covariate/blocker into the equation, I get slightly different statistics results in spite of the same coefficient values.

First example:

design <- model.matrix(~ group+covar) 
design
  (Intercept) group1 covar
  1           1      0                 0.2352934
  2           1      0                 0.2497340
  3           1      1                 0.3405839
  4           1      0                 0.3807273
  5           1      1                 0.3847266
  6           1      1                 0.3096032

          (Intercept)      group1        covar
653635      2.0700705 -0.04815855  2.031578
729737      4.7529897  1.02896914 -6.946158
100133331   2.2457020  0.73567908 -3.811348
100288069   1.8392895  0.19191507 -2.054743
100287934  -0.8755875  0.19874871  3.724989
400728      1.7266959  0.38078881 -1.048634

y <- voom(dgeObj,design) 
prefit <- lmFit(y,design)
fit <- eBayes(prefit)

            logFC  AveExpr         t      P.Value adj.P.Val         B
55007  -1.5740320 5.953331 -6.612608 0.0002165173 0.9999677 -4.014277

Second example:

design1 <- model.matrix(~ 0 + covar)
design1
  group0 group1 covar
1      1      0                 0.2352934
2      1      0                 0.2497340
3      0      1                 0.3405839
4      1      0                 0.3807273
5      0      1                 0.3847266
6      0      1                 0.3096032

                 CvsCON
  653635    -0.04815855
  729737     1.02896914
  100133331  0.73567908
  100288069  0.19191507
  100287934  0.19874871
  400728     0.38078881

contrasts <- makeContrasts(CvsCON = group1 - group0, levels = design1)
y1 <- voom(dgeObj,design1) 
prefit1 <- lmFit(y1,design1)
fit.cont <- contrasts.fit(prefit1, contrasts)
fit1 <- eBayes(fit.cont)

            logFC  AveExpr         t      P.Value adj.P.Val         B
55007  -1.5740320 5.953331 -6.756469 0.0001879744 0.9999673 -3.986743

Can anyone help me to figure out where the difference in topTable stat calculation is? Thank you very much in advance!

Answer 1

[gorenjakm]

Already figured it out. Thank you!

If the design matrix is not orthogonal, contrasts.fit() uses an approximation (The approximation only arises when there is an additive term). In the case that design matrix is orthogonal, contrasts.fit() is going to produce exact values as they are with the included intercept (~group).