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
1 answer
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
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.
Log in to answer this question.