[limma] LogFC interpretation

Question

[flappix]

Hi @ all,
I have cel files from 2 groups, 3 replicates each group. I want to use limma for the differential gene expression analysis. My script looks like this

prefix = "101053-0"
cels <- ReadAffy ( filenames = paste ( paste (prefix, c(4,5,6,1,2,3), sep=""), ".CEL", sep="" ),
                   sampleNames = c("N402_rep_1", "N402_rep_2", "N402_rep_3", "B36_rep_1", "B36_rep_2", "B36_rep_3") )



results <- expresso (cels, bgcorrect.method = "rma", normalize.method = "quantiles", pmcorrect.method = "pmonly", summary.method="medianpolish")
design <- model.matrix(~factor(c(1,1,1,2,2,2)))
fit <- lmFit (results, design=design)
ebayes <- eBayes (fit)
topTable (ebayes)

The first gene of the list is this

                  logFC  AveExpr          t      P.Value    adj.P.Val        B
An00g11929_at -7.347305 8.465857 -110.96206 6.249263e-18 9.095178e-14 27.44659

So logFC is reported as -7.34 But when I look at the absolute expression values

> exprs(results)["An00g11929_at", ]
N402_rep_1 N402_rep_2 N402_rep_3  B36_rep_1  B36_rep_2  B36_rep_3 
 12.127017  12.148645  12.142866   4.683787   4.842077   4.850748

The logFC is more about -1.5 than -7.34
How can this be?

I read that limma expect log2-transformed values. I also read that using median polish to summarize already produces log2 transformed expression values. Anyway, I transformed the value by my self

exprs(results) <- log2(exprs(results))

And now the logFC is what I would expect from the un-transformed values

                   logFC  AveExpr         t      P.Value    adj.P.Val        B
An00g11929_at -1.3411349 2.931070 -80.01571 2.751338e-12 4.004298e-08 16.27317

But now the log2-transformed values do not fit to the logFC

N402_rep_1 N402_rep_2 N402_rep_3  B36_rep_1  B36_rep_2  B36_rep_3 
  3.600153   3.602724   3.602037   2.227675   2.275626   2.278207

So I'm very confused. Why do the absolute expression values do not match the logFC? Do I have to apply the log2-transformation on my own?

Would be very nice if someone could help :)

[cpad0112]

> n=c(12.127017,12.148645,12.142866)
> t=c(4.683787,4.842077,4.850748)
> mean(t)-mean(n)
[1] -7.347305

Answer 1

[russhh]

The data stored in exprs is already log-transformed. The intensities on the array will be in the thousands.

A _back-of-the envelope_ calculation is as follows:

estimate of fold-change from the de-transformed array data:

geometric_mean(array_data(group_2)) / geometric_mean(array_data(group_1))
~ typical_value(group_2) / typical_value(group_1)
~ (2 ^ log_transformed_typical_value (group_2)) / (2 ^ log_transformed_typical_value(group_1))
~ (2^5) / (2^12)
~ 0.008
~ 1 / 125

So you get (an inaccurately calculated) ~ 125-fold change

estimate of log-fold-change from the transformed array data (as stored in exprs)

arithmetic_mean(exprs(group_2)) - arithmetic_mean(exprs(group_1))
~ typical_value(group_2) - typcial_value(group_1)
~ 5 - 12
~ -7

So your estimate of fold-change is 2^-7 ~ 0.008; again 1/125

[flappix]

Thank you very much, I see my mistake.

Just to be sure: If I want to retrieve the absolute expression values, I have to use 2^exprs(results) ?

[russhh]

That would give you an approximation to the intensity values from the array. However, some of your steps (eg, quantile normalisation) will mean that they aren't exactly the same as the intensity values that were present on the original array.

Answer 2

[Devon Ryan]

The logFC of log2 transformed values are (approximately) the difference in the group means, which is indeed ~1.35 in your last example and ~7.3 the former example.