Hi all,

I have an intercept issue that I do not understand with limma.

Let's say that I have a simple model with one group 'cancer sub-type' which contains 4 factors (TNBC, NonTNBC, HER2, Control).

I wrote a model which has either an intercept or not. If I add an intercept, the Control factor is used as reference.

In term of contrast, it means that testing 'subtypeTNBC' in the model with intercept is the same than testing 'subtypeTNBC - subtypeControl' in the model without intercept. Both tests compare the mean expression of TNBC vs Control groups.

However, if I want now to compare two other groups, let's say NonTNBC and TNBC ... using the same contrasts (as the Intercept is expected to correspond to the Control group), I have different results ...

Many thanks for your help

N

design <- model.matrix(~ 0 + subtype, data=splan)
v <- voom(y, design, plot=FALSE)
fit <- lmFit(v, design)
contrast <-makeContrasts(subtypeTNBC - subtypeNonTNBC, levels=design)
fit <- contrasts.fit(fit, contrast)
eb.wI <- eBayes(fit)
res.wI <- topTable(eb.wI, number=1e6, adjust.method="BH")

design.noI <- model.matrix(~ subtype, data=splan)
v <- voom(y, design.noI, plot=FALSE)
fit <- lmFit(v, design.noI)
contrast <-makeContrasts(subtypeTNBC - subtypeNonTNBC, levels=design.noI)
fit <- contrasts.fit(fit, contrast)
eb.noI <- eBayes(fit)
res.noI <- topTable(eb.noI, number=1e6, adjust.method="BH")

Your code looks fine, so I would suspect any difference is due to an approximation in contrasts.fit - see the warning in ?contrasts.fit. Any difference should be small - the analysis with design.noI should be exact.

Edit: The warning mentions that both versions should be exact when using a one-way model, though this does not seem to be the case:

y <- matrix(rnorm(100000), ncol=10)
w <- 2^matrix(rnorm(100000), ncol=10)
g <- gl(5, 2)

X.in <- model.matrix(~g)
fit.in <- lmFit(y, X.in, weights=w)
con.in <- makeContrasts(g5 - g4, levels=X.in)
fit.in2 <- contrasts.fit(fit.in, con.in)
fit.in2 <- eBayes(fit.in2)
topTable(fit.in2)

X.out <- model.matrix(~0 + g)
fit.out <- lmFit(y, X.out, weights=w)
con.out <- makeContrasts(g5 - g4, levels=X.out)
fit.out2 <- contrasts.fit(fit.out, con.out)
fit.out2 <- eBayes(fit.out2)
topTable(fit.out2)

Thanks for your feedbacks.

I agree that  small differences could be observed but here, it's huge ... twice mode differentially expressed genes !

> colnames(fit$coefficients)
[1] "subtypeTNBC - subtypeNonTNBC"
> res.noI <- topTable(eb.noI, number=1e6, adjust.method="BH", coef=1)
> head(res.noI)
                       logFC  AveExpr          t      P.Value    adj.P.Val        B
ENSG00000091831.17 -6.217176 6.003674 -12.166565 1.662744e-12 8.190013e-08 18.26739
ENSG00000003989.12 -5.562645 5.743366 -11.666386 4.356657e-12 1.072958e-07 17.46547
ENSG00000173467.4  -7.040726 2.197144 -10.470702 4.900484e-11 8.045941e-07 14.64939
ENSG00000109436.7  -3.587246 7.342151  -9.876878 1.739701e-10 2.142268e-06 14.04120
ENSG00000151892.10 -4.852984 4.816535  -9.682996 2.656709e-10 2.617177e-06 13.52225
ENSG00000118513.14 -3.623776 3.977708  -9.382388 5.171364e-10 4.245345e-06 12.78780
> head(res.wI)
                       logFC    AveExpr          t      P.Value    adj.P.Val        B
ENSG00000091831.17 -6.217176  6.0036736 -11.497290 6.072332e-12 2.990988e-07 16.94490
ENSG00000003989.12 -5.562645  5.7433655  -9.621528 3.041428e-10 7.490428e-06 13.26293
ENSG00000109436.7  -3.587246  7.3421515  -9.101191 9.746033e-10 1.215850e-05 12.34451
ENSG00000151892.10 -4.852984  4.8165346  -8.917986 1.481058e-09 1.215850e-05 11.81884
ENSG00000171428.9  -4.148381  2.9173688  -8.935683 1.422108e-09 1.215850e-05 11.61055
ENSG00000233823.1  -5.547515 -0.9911218  -9.030907 1.143732e-09 1.215850e-05 11.32991
> length(which(res.noI$adj.P.Val<0.05))
[1] 3507
> length(which(res.wI$adj.P.Val<0.05))
[1] 1853