Dear all,

I have a question about differential gene expression calculation using limma. I have two disease groups - D1 (n=6) and D2 (n=7) that I have first compared to each other and obtained about 30 genes that have adj.P-value < 0.05. Please see the code that I've used below.

> design <- model.matrix(~0 + Treat)
> fit <- lmFit(eset,design)
> cm <- makeContrasts(coef1 = D1 - D2, levels=design)
> fit2 <- contrasts.fit(fit, cm)
> fit2 <- eBayes(fit2)

After I had calculated that, I obtained some healthy samples - H (n=23) to compare each of the disease types to. However, now when I calculate D1 - D2, I get a lot more significant genes (~200) with adj.P-value < 0.05 as opposed to the 30 genes I was getting before with the same comparison. The sample numbers or anything else for my D1 and D2 groups haven't changed at all, I've only added the H samples. Please see code below: 

> design <- model.matrix(~0 + Treat)
> fit <- lmFit(eset,design)
> cm <- makeContrasts(coef1 = D1 - D2, coef2 = D1 - H, coef3= D2 - H, levels=design)
> fit2 <- contrasts.fit(fit, cm)
> fit2 <- eBayes(fit2)

I was wondering if anybody knew the reason for this discrepancy in the numbers of significant genes (with adj.P-value < 0.05) for 'coef1' from both of these codes. 

Thank you,

Akul

The answer by svlachavas is not correct in several ways.

First, the reason you get more significant genes when you include the healthy samples is because you have increased your degrees of freedom substantially by adding in 23 more samples. In other words, when you compute a contrast, the numerator is the difference in means between the two groups, and the denominator is a measure of the intra-group variability. The measure of intra-group variability is in essence an average of the intra-group variability over all the groups in your model. When you add a group with 23 samples, you have added another group that you can use to estimate intra-group variability. Sample variance estimators are not very efficient, which means that you need relatively more observations for the estimate to converge to the population parameter that you are trying to estimate, so adding in 23 more samples is a big deal.

If you increase your degrees of freedom substantially, then the variance estimate will be smaller, which will reduce the denominator for all of your contrasts. If you reduce the denominator while keeping the numerator constant, your statistic gets bigger, so comparisons that were not significant with just the two groups are now significant with three because you have made the denominator smaller by adding in the extra group. This has nothing to do with healthy being a 'baseline'. It is only due to the increase in observations.

Second, if you fit the model that svlachavas recommends:

design <- model,matrix(~Treat)

and then fit the contrast D1 - D2, you are NOT comparing anything to the 'healthy baseline'. This is simple algebra. The contrast you are fitting is

(D1 - H) - (D2 - H)

Which is easy to see will reduce to

D1 - D2

And the healthy samples play no part in the comparison.

Dear Akulsinghania,

in my opinion the most possible reason for the discrepancy in the number of significant genes, is the fact that in your first design matrix you didnt included the Healthy samples, so now your factor has one more level that it has to be taken into account. Also ,regarding your code you have to "refit" in way with contrasts.fit each disease group with the control group, in order for each of your coefficients to represent the log-fold change of your two disease groups versous your control group. In other words,

when you used cm<- makeContrasts(coef=D1-D2, levels=design) in the first code, you dont compute the difference between two two disease groups while your baseline group are the healthy samples. Thus, you could use:

Treat <- factor("your factor assignment", levels=c("H","D1","D2 )) # so the baseline level is the Healthy samples-group

design <- model,matrix(~Treat)

fit <- lmFit(eset, design) 
fit2 <- eBayes(fit, trend=TRUE) # trend argument optional if you have performed any initial non-intensity filtering prior testing

where you can see with head(fit2$coefficients) that the two coefficients represent each of your two disease groups vs Healthy group

Furthermore, if you intend to compare the two disease groups, in respect to the Healthy samples(while the "baseline level" is "H"),

you can move after the same design matrix with:

con <- makeContrasts("D1-D2", levels=design) # or D2-D1, as the coefficient output of design matrix
fit2 <- contrasts.fit(fit, contrasts=con)
ebayes <- eBayes(fit2, trend=TRUE)

The only think that i dont know for sure is if the difference in your number of samples(unequal number), hampers the analysis