limma fold change problem

Question

[Learner]

Hello,

I know this has been asked many times but I did not figure it out programmatically. So I tried to make a small reproducible example to see if someone knows about it

Here is a data

df <- structure(list(control1 = c(28.2494785829057, 28.1476353572629, 
27.1929113349158, 25.6702823986081, 26.8927288390394, 0, 26.7207761577158, 
26.9122794096471, 31.7603001024038, 29.6522061602385, 28.6278754368064, 
28.6577526656622, 24.9082071983177, 26.7735417008271, 27.0352760202938
), control2 = c(28.2584923577297, 28.1021446453534, 26.2233747158687, 
26.0274297329297, 26.9801099459161, 25.358504222508, 26.7976110787887, 
27.2493429511711, 31.7579581611434, 29.7188834878171, 28.5234274780504, 
28.6846850151125, 0, 27.0461438308603, 27.2389971036243), treatment1 = c(28.8868697499915, 
28.3221674791755, 26.1941520812899, 25.1512405315526, 26.8643148321299, 
25.2739455625142, 26.3634923544076, 27.3338132331068, 31.4077578366191, 
29.6662104574918, 28.3396435112577, 28.5763261654449, 24.9068319368666, 
26.8338833938774, 27.3473948352458), treatment2 = c(28.6182441737351, 
28.3364546698625, 26.5573504729984, 25.4954893749752, 26.9793468747426, 
25.7376090355345, 26.6356105379829, 27.4102977563763, 31.5699592291352, 
29.8796437292758, 28.1989170942089, 28.5062377808235, 24.6275067774911, 
26.8883222436102, 26.8426609128988)), .Names = c("control1", 
"control2", "treatment1", "treatment2"), row.names = c("guk-1", 
"dab-1", "swsn-9", "row-1", "ZC434.7", "rad-50", "sin-3", "hmp-1", 
"dhs-3", "nra-2", "saps-1", "D1081.7", "aph-2", "F15B9.10", "C27A7.5"
), class = "data.frame")

This data has 2 control and 2 treatment conditions. I want to get fold change and p-values based on Limma This is continuous data which is suitable for limma, so please don't ask redundant question like what type of data etc .

Here is what I do

design <- model.matrix(~c(rep(0,2),rep(1,2)))
fit <- lmFit(df, design)
fit2 <- eBayes(fit)
t <- topTable(fit2, coef=2, n=Inf)

which gives me this

               logFC  AveExpr          t    P.Value adj.P.Val         B
dab-1     0.20442107 28.22710  3.9514485 0.03685753 0.2271467 -3.905281
guk-1     0.49857149 28.50327  3.9394325 0.03711921 0.2271467 -3.906530
saps-1   -0.30637115 28.42247 -3.4043286 0.05179096 0.2271467 -3.971773
dhs-3    -0.27027060 31.62399 -3.1719049 0.06057245 0.2271467 -4.007078
row-1    -0.52549111 25.58611 -2.3700569 0.11096701 0.3051381 -4.174433
D1081.7  -0.12993687 28.60625 -2.2558390 0.12205525 0.3051381 -4.205428
sin-3    -0.25964217 26.62937 -1.9634112 0.15745464 0.3298828 -4.294491
hmp-1     0.29124431 27.22643  1.8412188 0.17593747 0.3298828 -4.336009
rad-50   12.82652519 19.09251  1.1581479 0.34133736 0.5256822 -4.609835
aph-2    12.31306576 18.61064  1.1320028 0.35045482 0.5256822 -4.621237
nra-2     0.08738227 29.72924  0.8052242 0.48723341 0.6555833 -4.761318
swsn-9   -0.33239175 26.54195 -0.7312325 0.52446668 0.6555833 -4.791007
F15B9.10 -0.04873995 26.88547 -0.3740694 0.73651132 0.8498207 -4.908211
ZC434.7  -0.01458854 26.92913 -0.1849504 0.86662644 0.8745856 -4.943976
C27A7.5  -0.04210869 27.11608 -0.1737547 0.87458562 0.8745856 -4.945376

if I look at this post, R: How to convert log2FC (Fold Change) obtained by limma's topTable() function to FC I should be able to get the answer by one of the following example but I don't

example1 <- log2(rowMeans(df[,c(3,4)])/(rowMeans(df[,c(1,2)])))

example2 <- log2(rowMeans(df[,c(3,4)])-(rowMeans(df[,c(1,2)])))

example3 <- log(rowMeans(df[,c(3,4)])/(rowMeans(df[,c(1,2)])))

So I have two question,

1- the way I am using the limma is right? (especially the design I make).

2- how really the fold change is calculated. Please avoid words , just show it by R

Answer 1

[russhh]

limma will assume your data is already log2 transformed. So, re 2:

example4 <- rowMeans(df[, 3:4]) - rowMeans(df[, 1:2])

gives the log2-fold change. If your data has not been transformed into log-space, you should do this yourself prior to calling limma.

The way you are using limma is ok, it's just a bit ugly: df and t are internal R functions, you shouldn't really be overwriting their defintiions; and you model.matrix will be a bit more self-explanatory if you used a factor to distinguish treatment from control, ie:

treat <- factor(rep(c('control', 'treatment'), each = 2), levels = c('control', 'treatment'))
design <- model.matrix(~ treat)

all the rest is fine

[Learner]

@russhh thanks for the nice and clear answer right to the point which i liked don't you think that we should get treatment/control for fold change ? also do you know how they calculate p-value in limma?

[russhh]

re treatment / control values; these aren't provided by default but then, they're superfluous: just exponentiate the log2 FC: 2 ^ logFC

... It's kind of assumed that you have street-fighting levels of maths.

The p-values are a different kettle of fish. You'd be better looking into how p-values are calculated in basic t-tests / ANOVA, before learning how the limma/eBayes-regularised p-values (and subsequently the multiple-comparison adjusted p-values) are determined.

[fr]

(~4 years after.... sorry!) Where does it say that limma assumes log2 transformation?